source: trunk/LMDZ.MARS/libf/phymars/soil.F @ 2917

Last change on this file since 2917 was 2900, checked in by romain.vande, 3 years ago

Mars PCM:
Following r2896, further implementation of subslope parametrisation.
Carefull ! This is a devolpment revision and it still need improvements and tests.
However, this commit should not change anything for nslope=1.
Only nslope=1 is possible for now!
Utilitaries are not adapted yet.
Dimension of some variables (30 listed below) is changed to add the nslope dimension.
Outputs (diagfi and restartfi) are not changed yet.
nslope is now read and written in the startfi.nc

Changed variables:
_ In surfdat_h:

  • qsurf
  • tsurf
  • watercap
  • emis
  • capcal
  • fluxgrd

_ In comsoil_h

  • tsoil
  • mthermdiff
  • thermdiff
  • coefd
  • alph
  • beta

_ In dimradmars_mod

  • albedo
  • fluxrad_sky
  • fluxrad

_ In physiq_mod

  • inertiesoil
  • fluxsurf_lw
  • fluxsurf_dn_sw
  • dqsurf
  • zdtsurf
  • zdtsdif
  • zdtsurfc
  • zdqsdif
  • zdqsc
  • dwatercap
  • dwatercap_dif
  • zflubid
  • fluxsurf_dn_sw_tot
  • inertiedat_tmp

New variables called VAR_meshavg correspond to the mesh average over all the subslope distribution of the variable
File size: 6.7 KB
Line 
1      subroutine soil(ngrid,nsoil,firstcall,
2     &          therm_i,
3     &          timestep,tsurf,tsoil,
4     &          capcal,fluxgrd)
5
6      use comsoil_h, only: layer, mlayer, volcapa,
7     &                     mthermdiff, thermdiff, coefq,
8     &                     coefd, alph, beta, mu,flux_geo
9      use surfdat_h, only: watercaptag, inert_h2o_ice
10      use comslope_mod, ONLY: nslope
11      implicit none
12
13!-----------------------------------------------------------------------
14!  Author: Ehouarn Millour
15!
16!  Purpose: Compute soil temperature using an implict 1st order scheme
17! 
18!  Note: depths of layers and mid-layers, soil thermal inertia and
19!        heat capacity are commons in comsoil_h
20!-----------------------------------------------------------------------
21
22#include "callkeys.h"
23
24c-----------------------------------------------------------------------
25!  arguments
26!  ---------
27!  inputs:
28      integer ngrid     ! number of (horizontal) grid-points
29      integer nsoil     ! number of soil layers
30      logical firstcall ! identifier for initialization call
31      real therm_i(ngrid,nsoil,nslope) ! thermal inertia
32      real timestep         ! time step
33      real tsurf(ngrid,nslope)   ! surface temperature
34! outputs:
35      real tsoil(ngrid,nsoil,nslope) ! soil (mid-layer) temperature
36      real capcal(ngrid,nslope) ! surface specific heat
37      real fluxgrd(ngrid,nslope) ! surface diffusive heat flux
38
39! local variables:
40      integer ig,ik,islope
41
42! 0. Initialisations and preprocessing step
43      if (firstcall.or.tifeedback) then
44      ! note: firstcall is set to .true. or .false. by the caller
45      !       and not changed by soil.F
46! 0.1 Build mthermdiff(:), the mid-layer thermal diffusivities
47      do ig=1,ngrid
48       do islope = 1,nslope
49        if (watercaptag(ig)) then
50          do ik=0,nsoil-1
51! If we have permanent ice, we use the water ice thermal inertia from ground to last layer.
52              mthermdiff(ig,ik,islope)=inert_h2o_ice*
53     &            inert_h2o_ice/volcapa
54          enddo
55        else
56          do ik=0,nsoil-1
57           mthermdiff(ig,ik,islope)=therm_i(ig,ik+1,islope)*
58     &         therm_i(ig,ik+1,islope)/volcapa
59          enddo
60        endif
61      enddo
62      enddo
63
64#ifdef MESOSCALE
65      do ig=1,ngrid
66       do islope = 1,nslope
67        if ( therm_i(ig,1,islope) .ge. inert_h2o_ice ) then
68          print *, "limit max TI ", therm_i(ig,1,islope), inert_h2o_ice
69          do ik=0,nsoil-1
70               mthermdiff(ig,ik,islope)=inert_h2o_ice*
71     &    inert_h2o_ice/volcapa
72          enddo
73        endif
74       enddo
75      enddo
76#endif
77
78! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
79      do ig=1,ngrid
80       do islope = 1,nslope
81        do ik=1,nsoil-1
82      thermdiff(ig,ik,islope)=((layer(ik)-mlayer(ik-1))
83     &                        *mthermdiff(ig,ik,islope)
84     &                +(mlayer(ik)-layer(ik))
85     &                *mthermdiff(ig,ik-1,islope))
86     &                    /(mlayer(ik)-mlayer(ik-1))
87!       write(*,*),'soil: ik: ',ik,' thermdiff:',thermdiff(ig,ik)
88        enddo
89       enddo
90      enddo
91
92! 0.3 Build coefficients mu, q_{k+1/2}, d_k, alpha_k and capcal
93      ! mu
94      mu=mlayer(0)/(mlayer(1)-mlayer(0))
95
96      ! q_{1/2}
97      coefq(0)=volcapa*layer(1)/timestep
98        ! q_{k+1/2}
99        do ik=1,nsoil-1
100          coefq(ik)=volcapa*(layer(ik+1)-layer(ik))
101     &                 /timestep
102        enddo
103
104      do ig=1,ngrid
105       do islope = 1,nslope
106        ! d_k
107        do ik=1,nsoil-1
108          coefd(ig,ik,islope)=thermdiff(ig,ik,islope)
109     &        /(mlayer(ik)-mlayer(ik-1))
110        enddo
111       
112        ! alph_{N-1}
113        alph(ig,nsoil-1,islope)=coefd(ig,nsoil-1,islope)/
114     &                  (coefq(nsoil-1)+coefd(ig,nsoil-1,islope))
115        ! alph_k
116        do ik=nsoil-2,1,-1
117          alph(ig,ik,islope)=coefd(ig,ik,islope)/
118     &         (coefq(ik)+coefd(ig,ik+1,islope)*
119     &          (1.-alph(ig,ik+1,islope))+coefd(ig,ik,islope))
120        enddo
121
122        ! capcal
123! Cstar
124        capcal(ig,islope)=volcapa*layer(1)+
125     &         (thermdiff(ig,1,islope)/(mlayer(1)-mlayer(0)))*
126     &         (timestep*(1.-alph(ig,1,islope)))
127! Cs
128        capcal(ig,islope)=capcal(ig,islope)/
129     &            (1.+mu*(1.0-alph(ig,1,islope))*
130     &            thermdiff(ig,1,islope)/mthermdiff(ig,0,islope))
131!      write(*,*)'soil: ig=',ig,' capcal(ig)=',capcal(ig)
132      enddo ! islope
133      enddo ! of do ig=1,ngrid
134           
135      endif ! of if (firstcall.or.tifeedback)
136
137!  1. Compute soil temperatures
138      IF (.not.firstcall) THEN
139! First layer:
140      do islope = 1,nslope
141      do ig=1,ngrid
142        tsoil(ig,1,islope)=(tsurf(ig,islope)+mu*beta(ig,1,islope)*
143     &      thermdiff(ig,1,islope)/mthermdiff(ig,0,islope))/
144     &      (1.+mu*(1.0-alph(ig,1,islope))*
145     &      thermdiff(ig,1,islope)/mthermdiff(ig,0,islope))
146      enddo
147! Other layers:
148      do ik=1,nsoil-1
149        do ig=1,ngrid
150          tsoil(ig,ik+1,islope)=alph(ig,ik,islope)*
151     &    tsoil(ig,ik,islope)+beta(ig,ik,islope)
152        enddo
153      enddo
154       enddo ! islope
155      ENDIF! of if (.not.firstcall)
156
157!  2. Compute beta coefficients (preprocessing for next time step)
158! Bottom layer, beta_{N-1}
159       do islope = 1,nslope
160      do ig=1,ngrid
161        beta(ig,nsoil-1,islope)=coefq(nsoil-1)*tsoil(ig,nsoil,islope)
162     &      /(coefq(nsoil-1)+coefd(ig,nsoil-1,islope))
163     &      +flux_geo(ig)/(coefq(nsoil-1) + coefd(ig,nsoil-1,islope))
164      enddo
165   
166! Other layers
167      do ik=nsoil-2,1,-1
168        do ig=1,ngrid
169          beta(ig,ik,islope)=(coefq(ik)*tsoil(ig,ik+1,islope)+
170     &                 coefd(ig,ik+1,islope)*beta(ig,ik+1,islope))/
171     &                 (coefq(ik)+coefd(ig,ik+1,islope)*
172     &                 (1.0-alph(ig,ik+1,islope))
173     &                  +coefd(ig,ik,islope))
174        enddo
175      enddo
176
177!  3. Compute surface diffusive flux & calorific capacity
178      do ig=1,ngrid
179! Cstar
180!        capcal(ig)=volcapa(ig,1)*layer(ig,1)+
181!     &              (thermdiff(ig,1)/(mlayer(ig,1)-mlayer(ig,0)))*
182!     &              (timestep*(1.-alph(ig,1)))
183! Fstar
184        fluxgrd(ig,islope)=(thermdiff(ig,1,islope)/
185     &              (mlayer(1)-mlayer(0)))*
186     &              (beta(ig,1,islope)+(alph(ig,1,islope)-1.0)
187     &              *tsoil(ig,1,islope))
188
189!        mu=mlayer(ig,0)/(mlayer(ig,1)-mlayer(ig,0))
190!        capcal(ig)=capcal(ig)/(1.+mu*(1.0-alph(ig,1))*
191!     &                         thermdiff(ig,1)/mthermdiff(ig,0))
192! Fs
193        fluxgrd(ig,islope)=fluxgrd(ig,islope)+(capcal(ig,islope)
194     &              /timestep)*
195     &              (tsoil(ig,1,islope)*
196     &              (1.+mu*(1.0-alph(ig,1,islope))*
197     &               thermdiff(ig,1,islope)/mthermdiff(ig,0,islope))
198     &               -tsurf(ig,islope)-mu*beta(ig,1,islope)*
199     &                thermdiff(ig,1,islope)/mthermdiff(ig,0,islope))
200      enddo
201      enddo ! islope
202      end
203
Note: See TracBrowser for help on using the repository browser.