soil_temp_arp.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     #########
      SUBROUTINE soil_temp_arp(PTSTEP,PA,PB,PC,PGAMMAT,PTDEEP,PSODELX,PTG)
!     ############################################################
!
!!****  *SOIL_TEMP_ARP*  
!
!!    PURPOSE
!!    -------
!     This subroutine solves the ARPEGE 1-d surface and deep force-restore 
!     'PTG' using the backward-difference scheme (implicit) as in soil_heatdiff. 
!     The eqs are solved rapidly by taking advantage of the
!     fact that the matrix is tridiagonal. Soln to the eqs:
!
!                   dTi       S1   
!                   --- =  Ct --  (Gi - Gi+1)
!                   dt        Si
!
!     with        |  G1 = Rn-H-LE
!                 |
!                 |         2Pi         1
!                 |  Gi = ------  -------------- (Ti-1 - Ti) 
!                 |       Ct*Day  S1 (Si-1 + Si)
!                 
!
!     where Si = pulsation depth, i=1 is the surface
!     
!!**  METHOD
!!    ------
!
!     Direct calculation
!
!!    EXTERNAL
!!    --------
!
!     None
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    USE MODD_SURF_PAR
!!    USE MODI_TRIDIAG_GROUND
!!      
!!    REFERENCE
!!    ---------
!!      
!!    AUTHOR
!!    ------
!!      B. Decharme          * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    21/01/09   B. Decharme
!
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_csts,       ONLY : xpi, xday
!
USE modd_surf_par,   ONLY : xundef
!
USE modi_tridiag_ground
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!
REAL, INTENT(IN)                   :: ptstep ! Model time step (s)
!
REAL, DIMENSION(:), INTENT(IN)     :: pa,pb,pc ! terms for the linearization of Ts(t)
!
REAL, DIMENSION(:), INTENT(IN)     :: ptdeep, pgammat
!                                     PTDEEP   = Deep soil temperature (prescribed)
!                                                which models heating/cooling from
!                                                below the diurnal wave penetration
!                                               (surface temperature) depth.
!                                      PGAMMAT  = Deep soil heat transfer coefficient:
!                                                assuming homogeneous soil so that
!                                                this can be prescribed in units of 
!                                                (1/days): associated time scale with
!                                                PTDEEP.
REAL, DIMENSION(:), INTENT (IN)     ::  psodelx   ! Pulsation for each layer (Only used if LTEMP_ARP=True)
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: ptg
!                                      PTG    = soil temperature (K)
!
!
!*      0.2    declarations of local variables
!
INTEGER                             :: jj ! Loop counter
!
INTEGER                             :: inlvld ! Number of points and grid layers
!
REAL, DIMENSION(SIZE(PTG,1),SIZE(PTG,2))  :: ztgm, zfrcv, zamtrx, zbmtrx,     &
                                             zcmtrx, zlambda, zalpha  
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------
!
! Initialize local variables:
!
IF (lhook) CALL dr_hook('SOIL_TEMP_ARP',0,zhook_handle)
inlvld = SIZE(ptg(:,:),2)
!
ztgm(:,:)      = ptg(:,:)
!
zfrcv(:,:)   = 0.0
zamtrx(:,:)   = 0.0
zbmtrx(:,:)   = 0.0
zcmtrx(:,:)   = 0.0
zalpha(:,:)   = 0.0
zlambda(:,:)   = 0.0
!
!-------------------------------------------------------------------------------
!
! Calculate tri-diagonal matrix coefficients:
!
!
DO jj=1,inlvld-1
   zalpha(:,jj) = psodelx(jj)*xday/(ptstep*psodelx(1)*2.0*xpi)
   zlambda(:,jj) = 1.0/(psodelx(1)*(psodelx(jj)+psodelx(jj+1)))
ENDDO
zalpha(:,inlvld) = psodelx(inlvld)*xday/(ptstep*psodelx(1)*2.0*xpi)
!
!-------------------------------------------------------------------------------
!
! Upper BC
!
zamtrx(:,1) =  0.0
zbmtrx(:,1) =  pa(:)+(2.0*xpi/xday)*(1.0/(psodelx(1)*(psodelx(1)+psodelx(2)))-1.0)
zcmtrx(:,1) = -2.0*xpi/(xday*psodelx(1)*(psodelx(1)+psodelx(2)))
zfrcv(:,1) =  pb(:)*ztgm(:,1)+pc(:)-2.0*xpi*ztgm(:,2)/xday
!
!
! Interior Grid
!
DO jj=2,inlvld-1
   zamtrx(:,jj) = -zlambda(:,jj-1) 
   zbmtrx(:,jj) =  zalpha(:,jj) + zlambda(:,jj-1) + zlambda(:,jj)
   zcmtrx(:,jj) = -zlambda(:,jj)
   zfrcv(:,jj)  =  zalpha(:,jj)*ztgm(:,jj) 
ENDDO
!
! Lower BC: 2 currently accounted for, Either zero flux
! or a fixed temperature 'TDEEP' 
!
zamtrx(:,inlvld) = -zlambda(:,inlvld-1) 
zcmtrx(:,inlvld) =  0.0                
!
WHERE(ptdeep(:) /= xundef .AND. pgammat(:) /= xundef)
   zbmtrx(:,inlvld) =  zalpha(:,inlvld) + zlambda(:,inlvld-1) + ptstep*pgammat(:)/xday
   zfrcv(:,inlvld) =  zalpha(:,inlvld)*ztgm(:,inlvld) + ptstep*pgammat(:)*ptdeep(:)/xday
ELSEWHERE
   zbmtrx(:,inlvld) =  zalpha(:,inlvld) + zlambda(:,inlvld-1) 
   zfrcv(:,inlvld) =  zalpha(:,inlvld)*ztgm(:,inlvld) 
END WHERE
!
!-------------------------------------------------------------------------------
!
! Compute ZTGM (solution vector) 
! used for systems of equations involving tridiagonal 
! matricies.
!
 CALL tridiag_ground(zamtrx,zbmtrx,zcmtrx,zfrcv,ztgm)
!
! Update values in time:
!
ptg(:,:) = ztgm(:,:)
!
IF (lhook) CALL dr_hook('SOIL_TEMP_ARP',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE soil_temp_arp