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

Last change on this file since 3026 was 2919, checked in by llange, 22 months ago

PCM
Flux_geo is now correctled initialized (wasnot correctly read in the startfi) for the 3D and 1D
Minor fix in the pem (wrong name for a variable)
LL
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,islope)/
164     &      (coefq(nsoil-1) + coefd(ig,nsoil-1,islope))
165      enddo
166   
167! Other layers
168      do ik=nsoil-2,1,-1
169        do ig=1,ngrid
170          beta(ig,ik,islope)=(coefq(ik)*tsoil(ig,ik+1,islope)+
171     &                 coefd(ig,ik+1,islope)*beta(ig,ik+1,islope))/
172     &                 (coefq(ik)+coefd(ig,ik+1,islope)*
173     &                 (1.0-alph(ig,ik+1,islope))
174     &                  +coefd(ig,ik,islope))
175        enddo
176      enddo
177
178!  3. Compute surface diffusive flux & calorific capacity
179      do ig=1,ngrid
180! Cstar
181!        capcal(ig)=volcapa(ig,1)*layer(ig,1)+
182!     &              (thermdiff(ig,1)/(mlayer(ig,1)-mlayer(ig,0)))*
183!     &              (timestep*(1.-alph(ig,1)))
184! Fstar
185        fluxgrd(ig,islope)=(thermdiff(ig,1,islope)/
186     &              (mlayer(1)-mlayer(0)))*
187     &              (beta(ig,1,islope)+(alph(ig,1,islope)-1.0)
188     &              *tsoil(ig,1,islope))
189
190!        mu=mlayer(ig,0)/(mlayer(ig,1)-mlayer(ig,0))
191!        capcal(ig)=capcal(ig)/(1.+mu*(1.0-alph(ig,1))*
192!     &                         thermdiff(ig,1)/mthermdiff(ig,0))
193! Fs
194        fluxgrd(ig,islope)=fluxgrd(ig,islope)+(capcal(ig,islope)
195     &              /timestep)*
196     &              (tsoil(ig,1,islope)*
197     &              (1.+mu*(1.0-alph(ig,1,islope))*
198     &               thermdiff(ig,1,islope)/mthermdiff(ig,0,islope))
199     &               -tsurf(ig,islope)-mu*beta(ig,1,islope)*
200     &                thermdiff(ig,1,islope)/mthermdiff(ig,0,islope))
201      enddo
202      enddo ! islope
203      end
204
Note: See TracBrowser for help on using the repository browser.