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

Last change on this file since 3726 was 3726, checked in by emillour, 3 months ago

Mars PCM:
Turn "callkeys.h" into module "callkeys_mod.F90"
EM

File size: 6.8 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      use callkeys_mod, only: surfaceice_tifeedback, poreice_tifeedback
12      implicit none
13
14!-----------------------------------------------------------------------
15!  Author: Ehouarn Millour
16!
17!  Purpose: Compute soil temperature using an implict 1st order scheme
18
19!  Note: depths of layers and mid-layers, soil thermal inertia and
20!        heat capacity are commons in comsoil_h
21!-----------------------------------------------------------------------
22
23c-----------------------------------------------------------------------
24!  arguments
25!  ---------
26!  inputs:
27      integer ngrid     ! number of (horizontal) grid-points
28      integer nsoil     ! number of soil layers
29      logical firstcall ! identifier for initialization call
30      real therm_i(ngrid,nsoil,nslope) ! thermal inertia
31      real timestep         ! time step
32      real tsurf(ngrid,nslope)   ! surface temperature
33! outputs:
34      real tsoil(ngrid,nsoil,nslope) ! soil (mid-layer) temperature
35      real capcal(ngrid,nslope) ! surface specific heat
36      real fluxgrd(ngrid,nslope) ! surface diffusive heat flux
37
38! local variables:
39      integer ig,ik,islope
40
41! 0. Initialisations and preprocessing step
42      if(firstcall.or.surfaceice_tifeedback.or.poreice_tifeedback) then
43      ! note: firstcall is set to .true. or .false. by the caller
44      !       and not changed by soil.F
45! 0.1 Build mthermdiff(:), the mid-layer thermal diffusivities
46      do ig=1,ngrid
47       do islope = 1,nslope
48        if (watercaptag(ig)) then
49          do ik=0,nsoil-1
50! If we have permanent ice, we use the water ice thermal inertia from ground to last layer.
51              mthermdiff(ig,ik,islope)=inert_h2o_ice*
52     &            inert_h2o_ice/volcapa
53          enddo
54        else
55          do ik=0,nsoil-1
56           mthermdiff(ig,ik,islope)=therm_i(ig,ik+1,islope)*
57     &         therm_i(ig,ik+1,islope)/volcapa
58          enddo
59        endif
60      enddo
61      enddo
62
63#ifdef MESOSCALE
64      do ig=1,ngrid
65       do islope = 1,nslope
66        if ( therm_i(ig,1,islope) .ge. inert_h2o_ice ) then
67          print *, "limit max TI ", therm_i(ig,1,islope), inert_h2o_ice
68          do ik=0,nsoil-1
69               mthermdiff(ig,ik,islope)=inert_h2o_ice*
70     &    inert_h2o_ice/volcapa
71          enddo
72        endif
73       enddo
74      enddo
75#endif
76
77! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
78      do ig=1,ngrid
79       do islope = 1,nslope
80        do ik=1,nsoil-1
81      thermdiff(ig,ik,islope)=((layer(ik)-mlayer(ik-1))
82     &                        *mthermdiff(ig,ik,islope)
83     &                +(mlayer(ik)-layer(ik))
84     &                *mthermdiff(ig,ik-1,islope))
85     &                    /(mlayer(ik)-mlayer(ik-1))
86!       write(*,*),'soil: ik: ',ik,' thermdiff:',thermdiff(ig,ik)
87        enddo
88       enddo
89      enddo
90
91! 0.3 Build coefficients mu, q_{k+1/2}, d_k, alpha_k and capcal
92      ! mu
93      mu=mlayer(0)/(mlayer(1)-mlayer(0))
94
95      ! q_{1/2}
96      coefq(0)=volcapa*layer(1)/timestep
97        ! q_{k+1/2}
98        do ik=1,nsoil-1
99          coefq(ik)=volcapa*(layer(ik+1)-layer(ik))
100     &                 /timestep
101        enddo
102
103      do ig=1,ngrid
104       do islope = 1,nslope
105        ! d_k
106        do ik=1,nsoil-1
107          coefd(ig,ik,islope)=thermdiff(ig,ik,islope)
108     &        /(mlayer(ik)-mlayer(ik-1))
109        enddo
110       
111        ! alph_{N-1}
112        alph(ig,nsoil-1,islope)=coefd(ig,nsoil-1,islope)/
113     &                  (coefq(nsoil-1)+coefd(ig,nsoil-1,islope))
114        ! alph_k
115        do ik=nsoil-2,1,-1
116          alph(ig,ik,islope)=coefd(ig,ik,islope)/
117     &         (coefq(ik)+coefd(ig,ik+1,islope)*
118     &          (1.-alph(ig,ik+1,islope))+coefd(ig,ik,islope))
119        enddo
120
121        ! capcal
122! Cstar
123        capcal(ig,islope)=volcapa*layer(1)+
124     &         (thermdiff(ig,1,islope)/(mlayer(1)-mlayer(0)))*
125     &         (timestep*(1.-alph(ig,1,islope)))
126! Cs
127        capcal(ig,islope)=capcal(ig,islope)/
128     &            (1.+mu*(1.0-alph(ig,1,islope))*
129     &            thermdiff(ig,1,islope)/mthermdiff(ig,0,islope))
130!      write(*,*)'soil: ig=',ig,' capcal(ig)=',capcal(ig)
131      enddo ! islope
132      enddo ! of do ig=1,ngrid
133           
134      endif ! of if (firstcall.or.tifeedback)
135
136!  1. Compute soil temperatures
137      IF (.not.firstcall) THEN
138! First layer:
139      do islope = 1,nslope
140      do ig=1,ngrid
141        tsoil(ig,1,islope)=(tsurf(ig,islope)+mu*beta(ig,1,islope)*
142     &      thermdiff(ig,1,islope)/mthermdiff(ig,0,islope))/
143     &      (1.+mu*(1.0-alph(ig,1,islope))*
144     &      thermdiff(ig,1,islope)/mthermdiff(ig,0,islope))
145      enddo
146! Other layers:
147      do ik=1,nsoil-1
148        do ig=1,ngrid
149          tsoil(ig,ik+1,islope)=alph(ig,ik,islope)*
150     &    tsoil(ig,ik,islope)+beta(ig,ik,islope)
151        enddo
152      enddo
153       enddo ! islope
154      ENDIF! of if (.not.firstcall)
155
156!  2. Compute beta coefficients (preprocessing for next time step)
157! Bottom layer, beta_{N-1}
158       do islope = 1,nslope
159      do ig=1,ngrid
160        beta(ig,nsoil-1,islope)=coefq(nsoil-1)*tsoil(ig,nsoil,islope)
161     &      /(coefq(nsoil-1)+coefd(ig,nsoil-1,islope))
162     &      +flux_geo(ig,islope)/
163     &      (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.