Changeset 2539 for trunk/LMDZ.VENUS

Timestamp:

Jun 23, 2021, 5:55:48 PM (3 years ago)

Author:

emillour

Message:

Venus GCM:
Fix soil model so that it correctly computes soil temperatures at the first time step. And in the process make it 1+1=2 compliant.
VB+EM

File:

• trunk/LMDZ.VENUS/libf/phyvenus/soil.F (modified) (11 diffs)

trunk/LMDZ.VENUS/libf/phyvenus/soil.F

@@ -46 +46 @@
 c   -------------
 
-      use dimphy
+      use dimphy, only: klon
       IMPLICIT NONE
-#include "YOMCST.h"
-#include "dimsoil.h"
-#include "clesphys.h"
+      include "YOMCST.h"
+      include "dimsoil.h"
+      include "clesphys.h"
 
 c-----------------------------------------------------------------------
@@ -56 +56 @@
 c  ---------
 
-      REAL ptimestep
-      INTEGER knon
-      REAL ptsrf(klon),ptsoil(klon,nsoilmx)
-      REAL pcapcal(klon),pfluxgrd(klon)
+      REAL, intent(IN) :: ptimestep
+      INTEGER, intent(IN) :: knon
+      REAL, intent(IN) :: ptsrf(klon)
+      REAL, intent(OUT) :: ptsoil(klon,nsoilmx)
+      REAL, intent(OUT) :: pcapcal(klon),pfluxgrd(klon)
 
 c-----------------------------------------------------------------------
@@ -72 +73 @@
 c   local saved variables:
 c   ----------------------
-      REAL dz1(nsoilmx),dz2(nsoilmx)
-      REAL,allocatable :: zc(:,:),zd(:,:)
-      REAL lambda
-      SAVE dz1,dz2,zc,zd,lambda
-      LOGICAL firstcall
-      SAVE firstcall
-
-      DATA firstcall/.true./
-
-c-----------------------------------------------------------------------
-c   Depthts:
-c   --------
+      REAL,SAVE :: dz1(nsoilmx),dz2(nsoilmx)
+      REAL,allocatable,save :: zc(:,:),zd(:,:)
+      REAL,SAVE :: lambda
+      LOGICAL,SAVE :: firstcall=.true.
+
+c-----------------------------------------------------------------------
+c   Depths:
+c   -------
 
       REAL fz,rk,fz1,rk1,rk2
       fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
+
       pfluxgrd(:) = 0.
  
+      ! on Venus thermal inertia is assumed constant over the globe 
       DO ig = 1, knon
           ztherm_i(ig)   = inertie
@@ -95 +94 @@
       IF (firstcall) THEN
 
-      allocate(zc(klon,nsoilmx),zd(klon,nsoilmx))
+         allocate(zc(klon,nsoilmx),zd(klon,nsoilmx))
 
 c-----------------------------------------------------------------------
@@ -103 +102 @@
 
 c VENUS : A REVOIR !!!!
-         min_period=20000. ! en secondes
-         dalph_soil=2.    ! rapport entre les epaisseurs de 2 couches succ.
+         min_period=20000. ! in seconds
+         dalph_soil=2.    ! ratio between successive layer sizes
 
          OPEN(99,file='soil.def',status='old',form='formatted',err=9999)
@@ -115 +114 @@
 9999     CONTINUE
 
-c   la premiere couche represente un dixieme de cycle diurne
+c   The first soil layer depth, based on min_period
          fz1=sqrt(min_period/3.14)
 
@@ -137 +136 @@
      .               fz(rk1)*fz(rk2)*3.14,fz(rk)*fz(rk)*3.14
          ENDDO
+         
+c-----------------------------------------------------------------------
+c   Computation of the Cgrd and Dgrd coefficient for the next step:
+c   ---------------------------------------------------------------
+         DO jk=1,nsoilmx
+            zdz2(jk)=dz2(jk)/ptimestep
+         ENDDO
+
+         DO ig=1,knon
+            z1(ig)=zdz2(nsoilmx)+dz1(nsoilmx-1)
+            zc(ig,nsoilmx-1)=
+     $          zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1(ig)
+            zd(ig,nsoilmx-1)=dz1(nsoilmx-1)/z1(ig)
+         ENDDO
+
+         DO jk=nsoilmx-1,2,-1
+            DO ig=1,knon
+               z1(ig)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)
+     $            *(1.-zd(ig,jk)))
+               zc(ig,jk-1)=
+     s         (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))
+     $             *z1(ig)
+               zd(ig,jk-1)=dz1(jk-1)*z1(ig)
+            ENDDO
+         ENDDO
          firstcall =.false.
 
-      ELSE   !--not firstcall
+      ENDIF !--not firstcall
+
 c-----------------------------------------------------------------------
 c   Computation of the soil temperatures using the Cgrd and Dgrd
@@ -145 +170 @@
 c  -----------------------------------------------
 
-c    surface temperature
+c        temperature in the first soil layer
          DO ig=1,knon
             ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/
@@ -151 +176 @@
          ENDDO
 
-c   other temperatures
+c        temperatures in the other soil layers
          DO jk=1,nsoilmx-1
             DO ig=1,knon
@@ -158 +183 @@
          ENDDO
 
-      ENDIF !--not firstcall
 c-----------------------------------------------------------------------
 c   Computation of the Cgrd and Dgrd coefficient for the next step:
@@ -204 +228 @@
       ENDDO
 
-      RETURN
+      
       END