Changeset 1516 for trunk/LMDZ.GENERIC

Timestamp:

Mar 17, 2016, 9:45:36 AM (9 years ago)

Author:

emillour

Message:

Generic model:

• Added tests to ensure that soil model settings are adequate to resolve sub-surface diurnal and annual thermal waves.

EM

Location:

trunk/LMDZ.GENERIC

Files:

• README (modified) (1 diff)
• libf/phystd/comsoil_h.F90 (modified) (1 diff)
• libf/phystd/soil.F (modified) (3 diffs)
• libf/phystd/soil_settings.F (modified) (4 diffs)

trunk/LMDZ.GENERIC/README

@@ -1117 +1117 @@
 - Added a generic routine (in utilities directory) to extrapolate ice fields.
 - Added a generic version of streamfunction (in utilities directory).
+
+== 17/03/2016 == EM
+- Added tests to ensure that soil model settings are adequate to resolve
+  sub-surface diurnal and annual thermal waves.

trunk/LMDZ.GENERIC/libf/phystd/comsoil_h.F90

@@ -6 +6 @@
 !integer, parameter :: nsoilmx = 13 ! for z1=0.03 m, depth = 104.8 m => earth case
   integer, parameter :: nsoilmx = 18
+! Full soil layer depths are set as: layer(k)=lay1*alpha**(k-1) , k=1,nsoil
+! Mid soil layer depths are set as: mlayer(k)=lay1*alpha**(k-1/2) , k=0,nsoil-1 
+  real,parameter :: lay1=2.e-4 ! depth of first full soil layer
+  real,parameter :: alpha=2
 
   real,save,allocatable,dimension(:) :: layer      ! soil layer depths

trunk/LMDZ.GENERIC/libf/phystd/soil.F

@@ -4 +4 @@
      &          capcal,fluxgrd)
 
-      use comsoil_h, only: layer, mlayer, volcapa
+      use comsoil_h, only: layer, mlayer, volcapa, inertiedat
+      use comcstfi_mod, only: pi, daysec
+      use planete_mod, only: year_day
+      use comgeomfi_h, only: long, lati
 
       implicit none
@@ -49 +52 @@
 ! local variables:
       integer ig,ik
-
+      real :: inertia_min,inertia_max
+      real :: diurnal_skin ! diurnal skin depth (m)
+      real :: annual_skin ! anuual skin depth (m)
 
 ! 0. Initialisations and preprocessing step
@@ -118 +123 @@
       !write(*,*)'soil: ig=',ig,' capcal(ig)=',capcal(ig)
       enddo ! of do ig=1,ngrid
-            
+      
+      ! Additional checks: is the vertical discretization sufficient
+      ! to resolve diurnal and annual waves?
+      do ig=1,ngrid
+        ! extreme inertia for this column
+        inertia_min=minval(inertiedat(ig,:))
+        inertia_max=maxval(inertiedat(ig,:))
+        ! diurnal and annual skin depth
+        diurnal_skin=(inertia_min/volcapa)*sqrt(daysec/pi)
+        annual_skin=(inertia_max/volcapa)*sqrt(year_day*daysec/pi)
+        if (0.5*diurnal_skin<layer(1)) then
+        ! one should have the fist layer be at least half of diurnal skin depth
+          write(*,*) "soil Error: grid point ig=",ig
+          write(*,*) "            longitude=",long(ig)*(180./pi)
+          write(*,*) "             latitude=",lati(ig)*(180./pi)
+          write(*,*) "  first soil layer depth ",layer(1)
+          write(*,*) "  not small enough for a diurnal skin depth of ",
+     &                diurnal_skin
+          write(*,*) " change soil layer distribution (comsoil_h.F90)"
+          stop
+        endif
+        if (2.*annual_skin>layer(nsoil)) then
+        ! one should have the full soil be at least twice the diurnal skin depth
+          write(*,*) "soil Error: grid point ig=",ig
+          write(*,*) "            longitude=",long(ig)*(180./pi)
+          write(*,*) "             latitude=",lati(ig)*(180./pi)
+          write(*,*) "  total soil layer depth ",layer(nsoil)
+          write(*,*) "  not large enough for an annual skin depth of ",
+     &                annual_skin
+          write(*,*) " change soil layer distribution (comsoil_h.F90)"
+          stop
+        endif
+      enddo ! of do ig=1,ngrid
+      
       else ! of if (firstcall)
 

trunk/LMDZ.GENERIC/libf/phystd/soil_settings.F

@@ -2 +2 @@
 
 !      use netcdf
-      use comsoil_h, only: layer, mlayer, inertiedat, volcapa
+      use comsoil_h, only: layer, mlayer, inertiedat, volcapa,
+     &                     lay1, alpha
       use iostart, only: inquire_field_ndims, get_var, get_field,
      &                   inquire_field, inquire_dimension_length
@@ -70 +71 @@
       real,dimension(:),allocatable :: newval
 
-      real alpha,lay1 ! coefficients for building layers
+      real malpha,mlay1 ! coefficients for building layers
       real xmin,xmax ! to display min and max of a field
 
@@ -133 +134 @@
       ! default mlayer distribution, following a power law:
       !  mlayer(k)=lay1*alpha**(k-1/2)
-        lay1=2.e-4
-        alpha=2
         do iloop=0,nsoil-1
           mlayer(iloop)=lay1*(alpha**(iloop-0.5))
@@ -142 +141 @@
       ! Assuming layer distribution follows mid-layer law:
       ! layer(k)=lay1*alpha**(k-1)
-      lay1=sqrt(mlayer(0)*mlayer(1))
-      alpha=mlayer(1)/mlayer(0)
+      mlay1=sqrt(mlayer(0)*mlayer(1))
+      malpha=mlayer(1)/mlayer(0)
+      ! Check that these values are the same as those prescibed for mlayers
+      if ((abs(mlay1-lay1)/lay1).gt.1.e-6) then
+         write(*,*) "soil_settings error: mlay1=",mlay1
+         write(*,*) " does not match comsoil_h lay1=",lay1
+         stop
+      endif
+      if ((abs(malpha-alpha)/alpha).gt.1.e-6) then
+         write(*,*) "soil_settings error: malpha=",malpha
+         write(*,*) " does not match comsoil_h alpha=",alpha
+         stop
+      endif
       do iloop=1,nsoil
-        layer(iloop)=lay1*(alpha**(iloop-1))
+        layer(iloop)=mlay1*(malpha**(iloop-1))
       enddo
+
 ! 2. Volumetric heat capacity (note: it is declared in comsoil_h)
 ! ---------------------------