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

Last change on this file since 3807 was 3727, checked in by emillour, 2 months ago

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