Changeset 1538 for trunk/LMDZ.GENERIC/libf/phystd/soil_settings.F

Timestamp:

Apr 8, 2016, 4:57:00 PM (9 years ago)

Author:

emillour

Message:

Generic GCM:

• Made nsoilmx be no longer a "parameter" and thus added the possibility to define the number of subsurface layers nsoilmx, along with first layer thickness "lay1_soil" and companion coefficient "alpha_soil", in callphys.def at run time. As before (when these were hard-coded), these are such that the depth of soil mid-layers are: mlayer(k)=lay1_soil*alpha_soil(k-1/2), for k=0,nsoil-1

EM

File:

• trunk/LMDZ.GENERIC/libf/phystd/soil_settings.F (modified) (4 diffs)

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

@@ -3 +3 @@
 !      use netcdf
       use comsoil_h, only: layer, mlayer, inertiedat, volcapa,
-     &                     lay1, alpha
+     &                     lay1_soil, alpha_soil
       use iostart, only: inquire_field_ndims, get_var, get_field,
      &                   inquire_field, inquire_dimension_length
@@ -71 +71 @@
       real,dimension(:),allocatable :: newval
 
-      real malpha,mlay1 ! coefficients for building layers
+      real malpha,mlay1_soil ! coefficients for building layers
       real xmin,xmax ! to display min and max of a field
 
@@ -84 +84 @@
 
         if (dimlen.ne.nsoil) then
-          write(*,*)'soil_settings: Interpolation of soil temperature ',
-     &              'and thermal inertia will be required!'
         ! if dimlen doesn't match nsoil, then interpolation of
         ! soil temperatures and thermal inertia will be requiered
           interpol=.true.
-        ! allocate olmlayer  
-          allocate(oldmlayer(dimlen),stat=ierr)
-          if (ierr.ne.0) then
-            write(*,*) 'soil_settings: failed allocation of oldmlayer!'
-            stop
-          endif
         endif
+
+        ! allocate oldmlayer  
+        allocate(oldmlayer(dimlen),stat=ierr)
+        if (ierr.ne.0) then
+          write(*,*) 'soil_settings: failed allocation of oldmlayer!'
+          stop
+        endif
+
+        ! check if olmlayer distribution matches current one
+        call get_var("soildepth",oldmlayer,found)
+        if (found) then
+          malpha=oldmlayer(2)/oldmlayer(1)
+          if ((abs(malpha-alpha_soil)/alpha_soil).gt.1.e-6) then
+            ! alpha values are too different, intepolation needed
+            interpol=.true.
+          endif
+          ! check if loaded mid-layer depth value differs from current one
+          if (abs((oldmlayer(1)-lay1_soil*alpha_soil**(-1./2.))/
+     &             (lay1_soil*alpha_soil**(-1./2.))).gt.1.e-6) then
+            interpol=.true.
+          endif
+        endif
+
+        if (interpol) then
+          write(*,*)'soil_settings: Interpolation of soil temperature ',
+     &              'and thermal inertia will be required!'
+        endif
+        
 ! 1.2 Find out the # of dimensions <inertiedat> was defined as using
       ndims=inquire_field_ndims("inertiedat")
@@ -133 +153 @@
       if (interpol) then
       ! default mlayer distribution, following a power law:
-      !  mlayer(k)=lay1*alpha**(k-1/2)
+      !  mlayer(k)=lay1_soil*alpha_soil**(k-1/2)
         do iloop=0,nsoil-1
-          mlayer(iloop)=lay1*(alpha**(iloop-0.5))
+          mlayer(iloop)=lay1_soil*(alpha_soil**(iloop-0.5))
         enddo
       endif
 ! 1.5 Build layer(); following the same law as mlayer()
       ! Assuming layer distribution follows mid-layer law:
-      ! layer(k)=lay1*alpha**(k-1)
-      mlay1=sqrt(mlayer(0)*mlayer(1))
+      ! layer(k)=lay1_soil*alpha_soil**(k-1)
+      mlay1_soil=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
+      if ((abs(mlay1_soil-lay1_soil)/lay1_soil).gt.1.e-6) then
+         write(*,*) "soil_settings error: mlay1_soil=",mlay1_soil
+         write(*,*) " does not match comsoil_h lay1_soil=",lay1_soil
          stop
       endif
-      if ((abs(malpha-alpha)/alpha).gt.1.e-6) then
+      if ((abs(malpha-alpha_soil)/alpha_soil).gt.1.e-6) then
          write(*,*) "soil_settings error: malpha=",malpha
-         write(*,*) " does not match comsoil_h alpha=",alpha
+         write(*,*) " does not match comsoil_h alpha_soil=",alpha_soil
          stop
       endif
       do iloop=1,nsoil
-        layer(iloop)=mlay1*(malpha**(iloop-1))
+        layer(iloop)=mlay1_soil*(malpha**(iloop-1))
       enddo