source: LMDZ6/trunk/libf/phylmd/soil.f90 @ 5409

Last change on this file since 5409 was 5298, checked in by abarral, 7 weeks ago

Turn planete.h comsoil.h into module
Remove obsolete message related to /scratch/ common

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.8 KB
RevLine 
[996]1!
2! $Header$
3!
[3974]4SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, qsol, &
5     lon, lat, ptsoil, pcapcal, pfluxgrd)
[996]6 
[5285]7  USE yomcst_mod_h
[996]8  USE dimphy
9  USE mod_phys_lmdz_para
[1785]10  USE indice_sol_mod
[2311]11  USE print_control_mod, ONLY: lunout
[5273]12  USE dimsoil_mod_h, ONLY: nsoilmx
[5298]13  USE comsoil_mod_h
[996]14  IMPLICIT NONE
15
16!=======================================================================
17!
18!   Auteur:  Frederic Hourdin     30/01/92
19!   -------
20!
21!   Object:  Computation of : the soil temperature evolution
22!   -------                   the surfacic heat capacity "Capcal"
23!                            the surface conduction flux pcapcal
24!
[3974]25!   Update: 2021/07 : soil thermal inertia, formerly a constant value,
26!   ------   can also be now a function of soil moisture (F Cheruy's idea)
27!            depending on iflag_inertie, read from physiq.def via conf_phys_m.F90
28!            ("Stage L3" Eve Rebouillat, with E Vignon, A Sima, F Cheruy)
[996]29!
30!   Method: Implicit time integration
31!   -------
32!   Consecutive ground temperatures are related by:
33!           T(k+1) = C(k) + D(k)*T(k)  (*)
34!   The coefficients C and D are computed at the t-dt time-step.
35!   Routine structure:
36!   1) C and D coefficients are computed from the old temperature
37!   2) new temperatures are computed using (*)
38!   3) C and D coefficients are computed from the new temperature
39!      profile for the t+dt time-step
40!   4) the coefficients A and B are computed where the diffusive
41!      fluxes at the t+dt time-step is given by
42!             Fdiff = A + B Ts(t+dt)
43!      or     Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
44!             with F0 = A + B (Ts(t))
45!                 Capcal = B*dt
[5274]46!
[996]47!   Interface:
48!   ----------
49!
50!   Arguments:
51!   ----------
52!   ptimestep            physical timestep (s)
53!   indice               sub-surface index
54!   snow(klon)           snow
55!   ptsrf(klon)          surface temperature at time-step t (K)
[3974]56!   qsol(klon)           soil moisture (kg/m2 or mm)
57!   lon(klon)            longitude in radian
58!   lat(klon)            latitude in radian
[996]59!   ptsoil(klon,nsoilmx) temperature inside the ground (K)
60!   pcapcal(klon)        surfacic specific heat (W*m-2*s*K-1)
61!   pfluxgrd(klon)       surface diffusive flux from ground (Wm-2)
[5274]62!
[996]63!=======================================================================
64! Arguments
65! ---------
66  REAL, INTENT(IN)                     :: ptimestep
[3974]67  INTEGER, INTENT(IN)                  :: indice, knon !, knindex
[996]68  REAL, DIMENSION(klon), INTENT(IN)    :: snow
69  REAL, DIMENSION(klon), INTENT(IN)    :: ptsrf
[3974]70  REAL, DIMENSION(klon), INTENT(IN)    :: qsol
71  REAL, DIMENSION(klon), INTENT(IN)    :: lon
72  REAL, DIMENSION(klon), INTENT(IN)    :: lat
73
[996]74  REAL, DIMENSION(klon,nsoilmx), INTENT(INOUT) :: ptsoil
75  REAL, DIMENSION(klon), INTENT(OUT)           :: pcapcal
76  REAL, DIMENSION(klon), INTENT(OUT)           :: pfluxgrd
77
78!-----------------------------------------------------------------------
79! Local variables
80! ---------------
81  INTEGER                             :: ig, jk, ierr
82  REAL                                :: min_period,dalph_soil
83  REAL, DIMENSION(nsoilmx)            :: zdz2
84  REAL                                :: z1s
85  REAL, DIMENSION(klon)               :: ztherm_i
86  REAL, DIMENSION(klon,nsoilmx,nbsrf) :: C_coef, D_coef
87
88! Local saved variables
89! ---------------------
90  REAL, SAVE                     :: lambda
91!$OMP THREADPRIVATE(lambda)
92  REAL, DIMENSION(nsoilmx), SAVE :: dz1, dz2
93!$OMP THREADPRIVATE(dz1,dz2)
94  LOGICAL, SAVE                  :: firstcall=.TRUE.
95!$OMP THREADPRIVATE(firstcall)
96   
97!-----------------------------------------------------------------------
98!   Depthts:
99!   --------
100  REAL fz,rk,fz1,rk1,rk2
101  fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
102
103
104!-----------------------------------------------------------------------
105! Calculation of some constants
106! NB! These constants do not depend on the sub-surfaces
107!-----------------------------------------------------------------------
108
109  IF (firstcall) THEN
110!-----------------------------------------------------------------------
111!   ground levels
112!   grnd=z/l where l is the skin depth of the diurnal cycle:
113!-----------------------------------------------------------------------
114
115     min_period=1800. ! en secondes
116     dalph_soil=2.    ! rapport entre les epaisseurs de 2 couches succ.
117!$OMP MASTER
118     IF (is_mpi_root) THEN
119        OPEN(99,file='soil.def',status='old',form='formatted',iostat=ierr)
120        IF (ierr == 0) THEN ! Read file only if it exists
121           READ(99,*) min_period
122           READ(99,*) dalph_soil
[1575]123           WRITE(lunout,*)'Discretization for the soil model'
124           WRITE(lunout,*)'First level e-folding depth',min_period, &
[996]125                '   dalph',dalph_soil
126           CLOSE(99)
127        END IF
128     ENDIF
129!$OMP END MASTER
130     CALL bcast(min_period)
131     CALL bcast(dalph_soil)
132
133!   la premiere couche represente un dixieme de cycle diurne
134     fz1=SQRT(min_period/3.14)
135     
136     DO jk=1,nsoilmx
137        rk1=jk
138        rk2=jk-1
139        dz2(jk)=fz(rk1)-fz(rk2)
140     ENDDO
141     DO jk=1,nsoilmx-1
142        rk1=jk+.5
143        rk2=jk-.5
144        dz1(jk)=1./(fz(rk1)-fz(rk2))
145     ENDDO
146     lambda=fz(.5)*dz1(1)
[1575]147     WRITE(lunout,*)'full layers, intermediate layers (seconds)'
[996]148     DO jk=1,nsoilmx
149        rk=jk
150        rk1=jk+.5
151        rk2=jk-.5
[1575]152        WRITE(lunout,*)'fz=', &
[996]153             fz(rk1)*fz(rk2)*3.14,fz(rk)*fz(rk)*3.14
154     ENDDO
155
156     firstcall =.FALSE.
157  END IF
158
159
160!-----------------------------------------------------------------------
161!   Calcul de l'inertie thermique a partir de la variable rnat.
[2915]162!   on initialise a inertie_sic meme au-dessus d'un point de mer au cas
[996]163!   ou le point de mer devienne point de glace au pas suivant
164!   on corrige si on a un point de terre avec ou sans glace
165!
[2915]166!   iophys can be used to write the ztherm_i variable in a phys.nc file
167!   and check the results; to do so, add "CALL iophys_ini" in physiq_mod
168!   and add knindex to the list of inputs in all the calls to soil.F90
169!   (and to soil.F90 itself !)
[996]170!-----------------------------------------------------------------------
[2915]171
[996]172  IF (indice == is_sic) THEN
173     DO ig = 1, knon
[2915]174        ztherm_i(ig)   = inertie_sic
[996]175     ENDDO
[2915]176     IF (iflag_sic == 0) THEN
177       DO ig = 1, knon
178         IF (snow(ig) > 0.0) ztherm_i(ig)   = inertie_sno
179       ENDDO
180!      Otherwise sea-ice keeps the same inertia, even when covered by snow
181     ENDIF
182!    CALL iophys_ecrit_index('ztherm_sic', 1, 'ztherm_sic', 'USI', &
183!      knon, knindex, ztherm_i)
[996]184  ELSE IF (indice == is_lic) THEN
185     DO ig = 1, knon
[2915]186        ztherm_i(ig)   = inertie_lic
[996]187        IF (snow(ig) > 0.0) ztherm_i(ig)   = inertie_sno
188     ENDDO
[2915]189!    CALL iophys_ecrit_index('ztherm_lic', 1, 'ztherm_lic', 'USI', &
190!      knon, knindex, ztherm_i)
[996]191  ELSE IF (indice == is_ter) THEN
[3974]192     !
193     ! La relation entre l'inertie thermique du sol et qsol change d'apres
194     !   iflag_inertie, defini dans physiq.def, et appele via comsoil.h
195     !
[996]196     DO ig = 1, knon
[3974]197        ! iflag_inertie=0 correspond au cas inertie=constant, comme avant
198        IF (iflag_inertie==0) THEN         
199           ztherm_i(ig)   = inertie_sol
200        ELSE IF (iflag_inertie == 1) THEN
201          ! I = a_qsol * qsol + b  modele lineaire deduit d'une
202          ! regression lineaire I = a_mrsos * mrsos + b obtenue sur
203          ! sorties MO d'une simulation LMDZOR(CMIP6) sur l'annee 2000
204          ! sur tous les points avec frac_snow=0
205          ! Difference entre qsol et mrsos prise en compte par un
206          ! facteur d'echelle sur le coefficient directeur de regression:
207          ! fact = 35./150. = mrsos_max/qsol_max
208          ! et a_qsol = a_mrsos * fact (car a = dI/dHumidite)
209            ztherm_i(ig) = 30.0 *35.0/150.0 *qsol(ig) +770.0
210          ! AS : pour qsol entre 0 - 150, on a I entre 770 - 1820
211        ELSE IF (iflag_inertie == 2) THEN
212          ! deux regressions lineaires, sur les memes sorties, 
213          ! distinguant le type de sol : sable ou autre (limons/argile)
214          ! Implementation simple : regression type "sable" seulement pour
215          ! Sahara, defini par une "boite" lat/lon (NB : en radians !! )
216          IF (lon(ig)>-0.35 .AND. lon(ig)<0.70 .AND. lat(ig)>0.17 .AND. lat(ig)<0.52) THEN
217              ! Valeurs theoriquement entre 728 et 2373 ; qsol valeurs basses
218              ztherm_i(ig) = 47. *35.0/150.0 *qsol(ig) +728.  ! boite type "sable" pour Sahara
219          ELSE
220              ! Valeurs theoriquement entre 550 et 1940 ; qsol valeurs moyennes et hautes
221              ztherm_i(ig) = 41. *35.0/150.0 *qsol(ig) +505.
222          ENDIF
223        ELSE IF (iflag_inertie == 3) THEN
224          ! AS : idee a tester :
225          ! si la relation doit etre une droite,
226          ! definissons-la en fonction des valeurs min et max de qsol (0:150),
227          ! et de l'inertie (900 : 2000 ou 2400 ; choix ici: 2000)
228          ! I = I_min + qsol * (I_max - I_min)/(qsol_max - qsol_min)
229              ztherm_i(ig) = 900. + qsol(ig) * (2000. - 900.)/150.
230        ELSE         
231          WRITE (lunout,*) "Le choix iflag_inertie = ",iflag_inertie," n'est pas defini. Veuillez choisir un entier entre 0 et 3"
232        ENDIF
233     !
234     ! Fin de l'introduction de la relation entre l'inertie thermique du sol et qsol
235     !-------------------------------------------
236        !AS : donc le moindre flocon de neige sur un point de grid
237        ! fait que l'inertie du point passe a la valeur pour neige !
[996]238        IF (snow(ig) > 0.0) ztherm_i(ig)   = inertie_sno
[3974]239       
[996]240     ENDDO
[2915]241!    CALL iophys_ecrit_index('ztherm_ter', 1, 'ztherm_ter', 'USI', &
242!      knon, knindex, ztherm_i)
[996]243  ELSE IF (indice == is_oce) THEN
244     DO ig = 1, knon
[2915]245!       This is just in case, but SST should be used by the model anyway
246        ztherm_i(ig)   = inertie_sic
[996]247     ENDDO
[2915]248!    CALL iophys_ecrit_index('ztherm_oce', 1, 'ztherm_oce', 'USI', &
249!      knon, knindex, ztherm_i)
[996]250  ELSE
[1575]251     WRITE(lunout,*) "valeur d indice non prevue", indice
[2311]252     call abort_physic("soil", "", 1)
[996]253  ENDIF
254
255
256!-----------------------------------------------------------------------
257! 1)
258! Calculation of Cgrf and Dgrd coefficients using soil temperature from
259! previous time step.
260!
261! These variables are recalculated on the local compressed grid instead
262! of saved in restart file.
263!-----------------------------------------------------------------------
264  DO jk=1,nsoilmx
265     zdz2(jk)=dz2(jk)/ptimestep
266  ENDDO
267 
268  DO ig=1,knon
269     z1s = zdz2(nsoilmx)+dz1(nsoilmx-1)
270     C_coef(ig,nsoilmx-1,indice)= &
271          zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1s
272     D_coef(ig,nsoilmx-1,indice)=dz1(nsoilmx-1)/z1s
273  ENDDO
274 
275  DO jk=nsoilmx-1,2,-1
276     DO ig=1,knon
277        z1s = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk) &
278             *(1.-D_coef(ig,jk,indice)))
279        C_coef(ig,jk-1,indice)= &
280             (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*C_coef(ig,jk,indice)) * z1s
281        D_coef(ig,jk-1,indice)=dz1(jk-1)*z1s
282     ENDDO
283  ENDDO
284
285!-----------------------------------------------------------------------
286! 2)
287! Computation of the soil temperatures using the Cgrd and Dgrd
288! coefficient computed above
289!
290!-----------------------------------------------------------------------
291
292!    Surface temperature
293  DO ig=1,knon
294     ptsoil(ig,1)=(lambda*C_coef(ig,1,indice)+ptsrf(ig))/  &
295          (lambda*(1.-D_coef(ig,1,indice))+1.)
296  ENDDO
297 
298!   Other temperatures
299  DO jk=1,nsoilmx-1
300     DO ig=1,knon
301        ptsoil(ig,jk+1)=C_coef(ig,jk,indice)+D_coef(ig,jk,indice) &
302             *ptsoil(ig,jk)
303     ENDDO
304  ENDDO
305
306  IF (indice == is_sic) THEN
307     DO ig = 1 , knon
308        ptsoil(ig,nsoilmx) = RTT - 1.8
309     END DO
310  ENDIF
311
312!-----------------------------------------------------------------------
313! 3)
314! Calculate the Cgrd and Dgrd coefficient corresponding to actual soil
315! temperature
316!-----------------------------------------------------------------------
317  DO ig=1,knon
318     z1s = zdz2(nsoilmx)+dz1(nsoilmx-1)
319     C_coef(ig,nsoilmx-1,indice) = zdz2(nsoilmx)*ptsoil(ig,nsoilmx)/z1s
320     D_coef(ig,nsoilmx-1,indice) = dz1(nsoilmx-1)/z1s
321  ENDDO
322 
323  DO jk=nsoilmx-1,2,-1
324     DO ig=1,knon
325        z1s = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk) &
326             *(1.-D_coef(ig,jk,indice)))
327        C_coef(ig,jk-1,indice) = &
328             (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*C_coef(ig,jk,indice)) * z1s
329        D_coef(ig,jk-1,indice) = dz1(jk-1)*z1s
330     ENDDO
331  ENDDO
332
333!-----------------------------------------------------------------------
334! 4)
335! Computation of the surface diffusive flux from ground and
336! calorific capacity of the ground
337!-----------------------------------------------------------------------
338  DO ig=1,knon
339     pfluxgrd(ig) = ztherm_i(ig)*dz1(1)* &
340          (C_coef(ig,1,indice)+(D_coef(ig,1,indice)-1.)*ptsoil(ig,1))
341     pcapcal(ig)  = ztherm_i(ig)* &
342          (dz2(1)+ptimestep*(1.-D_coef(ig,1,indice))*dz1(1))
343     z1s = lambda*(1.-D_coef(ig,1,indice))+1.
344     pcapcal(ig)  = pcapcal(ig)/z1s
345     pfluxgrd(ig) = pfluxgrd(ig) &
346          + pcapcal(ig) * (ptsoil(ig,1) * z1s &
347          - lambda * C_coef(ig,1,indice) &
348          - ptsrf(ig)) &
349          /ptimestep
350  ENDDO
351   
352END SUBROUTINE soil
Note: See TracBrowser for help on using the repository browser.