Changeset 4245 for dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90

Timestamp:

Jun 24, 2020, 11:27:16 AM (4 years ago)

Author:

dubos

Message:

simple_physics : turn zc, zd, capcal, fluxgrd into local temporaries

File:

• dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90 (modified) (12 diffs)

dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90

@@ -3 +3 @@
   USE callkeys
   USE comgeomfi
-  USE surface
   IMPLICIT NONE
   PRIVATE
@@ -25 +24 @@
        &            pdu,pdv,pdt,pdpsrf) BIND(C, name='phyparam_phyparam')
     USE phys_const, ONLY : g, rcp, r, unjours
+    USE surface, ONLY : soil_forward, soil_backward
+    USE surface, ONLY : z0, inertie, emissiv, albedo  ! precomputed
+    USE surface, ONLY : tsurf, tsoil ! state variables
+!    USE surface, ONLY : capcal, fluxgrd, zc, zd ! should be temporaries
     USE turbulence, ONLY : vdif
     USE convection, ONLY : convadj
@@ -66 +69 @@
 
     !    Local variables :
-    REAL, DIMENSION(ngrid) :: mu0
-    INTEGER :: j,l,ig,nlevel,igout
-    LOGICAL :: lwrite
-    !
-    REAL :: zday, zdtime
-    REAL zh(ngrid,nlayer),z1,z2
-    REAL zzlev(ngrid,nlayer+1),zzlay(ngrid,nlayer)
-    REAL zdvfr(ngrid,nlayer),zdufr(ngrid,nlayer)
-    REAL zdhfr(ngrid,nlayer),zdtsrf(ngrid),zdtsrfr(ngrid)
-    REAL zflubid(ngrid),zpmer(ngrid)
-    REAL zpopsk(ngrid,nlayer)
+    REAL :: zh(ngrid,nlayer),      & ! potential temperature
+         &  zpopsk(ngrid,nlayer),  & ! Exner function
+         &  zzlev(ngrid,nlayer+1), & ! altitude of interfaces
+         &  zzlay(ngrid,nlayer),   & ! altitude of full levels
+         &  fluxrad(ngrid),        & ! radiative flux at surface
+         &  zc(ngrid, nsoilmx),    & ! LU coefficients for soil implicit solve
+         &  zd(ngrid, nsoilmx),    &
+         &  fluxgrd(ngrid),        & ! heat flux from deep soil
+         &  capcal(ngrid)          & ! effective heat capacity of soil
+         &  zdufr(ngrid,nlayer),   & ! partial tendencies for zonal velocity,
+         &  zdvfr(ngrid,nlayer),   & !   meridional velocity,
+         &  zdhfr(ngrid,nlayer),   & !   potential temperature,
+         &  zdtsrfr(ngrid),        & !   surface temperature
+         &  zdtsrf(ngrid),         & ! total tendency of surface temperature
+         &  zflubid(ngrid),        & ! radiative + deep soil fluxes
+         &  zpmer(ngrid),          & ! sea-level pressure
     REAL zdum1(ngrid,nlayer)
     REAL zdum2(ngrid,nlayer)
     REAL zdum3(ngrid,nlayer)
-    REAL fluxrad(ngrid)
+
+    INTEGER :: j,l,ig,nlevel,igout
+    LOGICAL :: lwrite
+    REAL    :: zday, zdtime, z1, z2
 
     WRITELOG(*,*) 'latitude0', ngrid, lati(1:2), lati(ngrid-1:ngrid)
@@ -88 +99 @@
     IF (ngrid.NE.ngridmax) THEN
        PRINT*,'STOP in inifis'
-       PRINT*,'Probleme de dimenesions :'
+       PRINT*,'Probleme de dimensions :'
        PRINT*,'ngrid     = ',ngrid
        PRINT*,'ngridmax   = ',ngridmax
@@ -159 +170 @@
     ENDDO
 
+    !-------------------------------------------------------------
+    !  soil temperatures : 1st half of implicit time integration
+    !  forward sweep from deep ground to surface
+    !  yields LU coefficients zc,zd and capcal, fluxgrd
+    !   ----------------------------------------------------------
+
+    IF (callsoil) THEN
+       CALL soil_forward(ngrid,nsoilmx, ptimestep, inertie, tsurf, tsoil, &
+         &            zc, zd, capcal, fluxgrd)
+
+!       CALL soil_new(ngrid,nsoilmx,ptimestep,inertie, &
+!            tsurf, tsoil, capcal,fluxgrd)
+!       CALL soil(ngrid,nsoilmx,.false.,inertie, &
+!            &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
+    ELSE
+       capcal(:)  = capcal_nosoil
+       fluxgrd(:) = 0.
+    ENDIF
+
+    IF(lverbose) THEN
+       WRITELOG(*,*) 'Surface Heat capacity, conduction Flux, Ts'
+       WRITELOG(*,*) capcal(igout), fluxgrd(igout), tsurf(igout)
+       LOG_DBG('phyparam')
+    ENDIF
+
     !-----------------------------------------------------------------------
     !    2. Compute radiative tendencies :
@@ -169 +205 @@
     !-----------------------------------------------------------------------
     !    3. Vertical diffusion (turbulent mixing):
-    !    -----------------------------------------
+    ! Kz is computed then vertical diffusion is integrated in time implicitly
+    ! using a linear relationship between surface heat flux and air temperature
+    ! in lowest level (Robin-type BC)
+    !    -------------------------------------------------------------------
     !
     IF(calldifv) THEN
@@ -212 +251 @@
        ENDDO
     ENDIF
+
+    !-------------------------------------------------------------
+    !   soil temperatures : 2nd half of implicit time integration
+    !   using updated tsurf as input
+    !   ----------------------------------------------------------
+
+    DO ig=1,ngrid
+       tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig)
+    ENDDO
+
+    WRITE(55,'(2e15.5)') zday,tsurf(ngrid/2+1)
+
+    IF (callsoil) THEN
+       CALL soil_backward(ngrid,nsoilmx, zc,zd, tsurf,tsoil)
+       IF(lverbose) THEN
+          WRITELOG(*,*) 'Surface Ts, dTs, dt'
+          WRITELOG(*,*) tsurf(igout), zdtsrf(igout), ptimestep
+          LOG_DBG('phyparam')
+       ENDIF
+    END IF
+
+
     !
     !-----------------------------------------------------------------------
@@ -243 +304 @@
        ENDDO
 
-    ENDIF
-
-    !-----------------------------------------------------------------------
-    !   On ajoute les tendances physiques a la temperature du sol:
-    !   ---------------------------------------------------------------
-
-    DO ig=1,ngrid
-       tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig)
-    ENDDO
-
-    WRITE(55,'(2e15.5)') zday,tsurf(ngrid/2+1)
-
-    !-----------------------------------------------------------------------
-    !   soil temperatures:
-    !   --------------------
-
-    IF (callsoil) THEN
-       CALL soil(ngrid,nsoilmx,.false.,inertie, &
-            &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
-       IF(lverbose) THEN
-          WRITELOG(*,*) 'Surface Heat capacity,conduction Flux, Ts, dTs, dt'
-          WRITELOG(*,*) capcal(igout), fluxgrd(igout), tsurf(igout), &
-               &        zdtsrf(igout), ptimestep
-          LOG_DBG('phyparam')
-       ENDIF
     ENDIF
 
@@ -309 +345 @@
     !$cython wrapper def alloc(ngrid, nlayer) : phy.phyparam_alloc(ngrid, nlayer)
     USE astronomy, ONLY : iniorbit
+    USE surface
     INTEGER, INTENT(IN), VALUE :: ngrid, nlayer
+    ! allocate precomputed arrays
+    ALLOCATE(rnatur(ngrid), albedo(ngrid), emissiv(ngrid))
+    ALLOCATE(z0(ngrid),inertie(ngrid))
     ! allocate arrays for internal state
     ALLOCATE(tsurf(ngrid))
     ALLOCATE(tsoil(ngrid,nsoilmx))
-    ! we could avoid the arrays below with a different implementation of surface / radiation / turbulence coupling
-    ALLOCATE(capcal(ngrid),fluxgrd(ngrid))
-    ALLOCATE(zc(ngrid,nsoilmx),zd(ngrid,nsoilmx))
-    ! allocate precomputed arrays
-    ALLOCATE(rnatur(ngrid), albedo(ngrid), emissiv(ngrid))
-    ALLOCATE(z0(ngrid),inertie(ngrid))
+    IF(.FALSE.) THEN ! arrays below are now local temporaries in phyparam
+       ALLOCATE(capcal(ngrid),fluxgrd(ngrid))
+       ALLOCATE(zc(ngrid,nsoilmx),zd(ngrid,nsoilmx))
+    END IF
     CALL iniorbit
   END SUBROUTINE alloc
@@ -326 +364 @@
     !$cython wrapper def precompute() : phy.phyparam_precompute()
     ! precompute time-independent arrays
+    USE surface
     rnatur(:)  = 1.
     inertie(:) = (1.-rnatur(:))*I_mer+rnatur(:)*I_ter
@@ -337 +376 @@
     !$cython wrapper def coldstart (ngrid, timestep): phy.phyparam_coldstart(ngrid, timestep)
     ! create internal state to start a run without a restart file
+    USE surface
     INTEGER, INTENT(IN), VALUE :: ngrid
     REAL, INTENT(IN),    VALUE :: ptimestep
@@ -343 +383 @@
     icount=0
     IF(callsoil) THEN
-       ! initializes zc, zd, capcal, fluxgrd
-       CALL soil(ngrid,nsoilmx,.TRUE.,inertie, &
-            &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
+       IF(.FALSE.) THEN
+          ! init_soil is now called by iniphyparam
+          ! initializes zc, zd, capcal, fluxgrd
+          CALL soil(ngrid,nsoilmx,.TRUE.,inertie, &
+               &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
+       END IF
+       IF(.FALSE.) THEN ! soil_forward is now called by phyparam
+          CALL soil_forward(ngrid, nsoilmx, ptimestep, inertie, tsurf, tsoil, &
+               &            zc, zd, capcal, fluxgrd)
+       END IF
     ELSE
        WRITELOG(*,*) 'WARNING!!! Thermal conduction in the soil turned off'