source: trunk/LMDZ.GENERIC/libf/phystd/soil.F @ 2537

Last change on this file since 2537 was 1543, checked in by emillour, 9 years ago

All models: Further adaptations to keep up with changes in LMDZ5 concerning
physics/dynamics separation:

  • dyn3d:
  • adapted gcm.F so that all physics initializations are now done in iniphysiq.
  • dyn3dpar:
  • adapted gcm.F so that all physics initializations are now done in iniphysiq.
  • updated calfis_p.F to follow up with changes.
  • copied over updated "bands.F90" from LMDZ5.
  • dynphy_lonlat:
  • calfis_p.F90, mod_interface_dyn_phys.F90, follow up of changes in phy_common/mod_* routines
  • phy_common:
  • added "geometry_mod.F90" to store information about the grid (replaces phy*/comgeomphy.F90) and give variables friendlier names: rlond => longitude , rlatd => latitude, airephy => cell_area, cuphy => dx , cvphy => dy
  • added "physics_distribution_mod.F90"
  • updated "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_mpi_data.F90", "mod_phys_lmdz_para.F90", "mod_phys_lmdz_mpi_transfert.F90", "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_omp_data.F90", "mod_phys_lmdz_omp_transfert.F90", "write_field_phy.F90" and "ioipsl_getin_p_mod.F90" to LMDZ5 versions.
  • phy[venus/titan/mars/std]:
  • removed "init_phys_lmdz.F90", "comgeomphy.F90"; adapted routines to use geometry_mod (longitude, latitude, cell_area, etc.)

EM

File size: 8.2 KB
Line 
1      subroutine soil(ngrid,nsoil,firstcall,lastcall,
2     &          therm_i,
3     &          timestep,tsurf,tsoil,
4     &          capcal,fluxgrd)
5
6      use comsoil_h, only: layer, mlayer, volcapa, inertiedat
7      use comcstfi_mod, only: pi
8      use time_phylmdz_mod, only: daysec
9      use planete_mod, only: year_day
10      use geometry_mod, only: longitude, latitude ! in radians
11
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,intent(in) :: ngrid       ! number of (horizontal) grid-points
28      integer,intent(in) :: nsoil       ! number of soil layers
29      logical,intent(in) :: firstcall ! identifier for initialization call
30      logical,intent(in) :: lastcall
31      real,intent(in) :: therm_i(ngrid,nsoil) ! thermal inertia
32      real,intent(in) :: timestep           ! time step
33      real,intent(in) :: tsurf(ngrid)   ! surface temperature
34! outputs:
35      real,intent(out) :: tsoil(ngrid,nsoil) ! soil (mid-layer) temperature
36      real,intent(out) :: capcal(ngrid) ! surface specific heat
37      real,intent(out) :: fluxgrd(ngrid) ! surface diffusive heat flux
38
39! local saved variables:
40      real,dimension(:,:),save,allocatable :: mthermdiff ! mid-layer thermal diffusivity
41      real,dimension(:,:),save,allocatable :: thermdiff ! inter-layer thermal diffusivity
42      real,dimension(:),save,allocatable :: coefq ! q_{k+1/2} coefficients
43      real,dimension(:,:),save,allocatable :: coefd ! d_k coefficients
44      real,dimension(:,:),save,allocatable :: alph ! alpha_k coefficients
45      real,dimension(:,:),save,allocatable :: beta ! beta_k coefficients
46      real,save :: mu
47!$OMP THREADPRIVATE(mthermdiff,thermdiff,coefq,coefd,alph,beta,mu)
48           
49! local variables:
50      integer ig,ik
51      real :: inertia_min,inertia_max
52      real :: diurnal_skin ! diurnal skin depth (m)
53      real :: annual_skin ! anuual skin depth (m)
54
55! 0. Initialisations and preprocessing step
56      if (firstcall) then
57      ! note: firstcall is set to .true. or .false. by the caller
58      !       and not changed by soil.F
59
60      ALLOCATE(mthermdiff(ngrid,0:nsoil-1)) ! mid-layer thermal diffusivity
61      ALLOCATE(thermdiff(ngrid,nsoil-1))    ! inter-layer thermal diffusivity
62      ALLOCATE(coefq(0:nsoil-1))              ! q_{k+1/2} coefficients
63      ALLOCATE(coefd(ngrid,nsoil-1))        ! d_k coefficients
64      ALLOCATE(alph(ngrid,nsoil-1))         ! alpha_k coefficients
65      ALLOCATE(beta(ngrid,nsoil-1))         ! beta_k coefficients
66
67! 0.1 Build mthermdiff(:), the mid-layer thermal diffusivities
68      do ig=1,ngrid
69        do ik=0,nsoil-1
70          mthermdiff(ig,ik)=therm_i(ig,ik+1)*therm_i(ig,ik+1)/volcapa
71!         write(*,*),'soil: ik: ',ik,' mthermdiff:',mthermdiff(ig,ik)
72        enddo
73      enddo
74
75! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
76      do ig=1,ngrid
77        do ik=1,nsoil-1
78      thermdiff(ig,ik)=((layer(ik)-mlayer(ik-1))*mthermdiff(ig,ik)
79     &                +(mlayer(ik)-layer(ik))*mthermdiff(ig,ik-1))
80     &                    /(mlayer(ik)-mlayer(ik-1))
81!       write(*,*),'soil: ik: ',ik,' thermdiff:',thermdiff(ig,ik)
82        enddo
83      enddo
84
85! 0.3 Build coefficients mu, q_{k+1/2}, d_k, alpha_k and capcal
86      ! mu
87      mu=mlayer(0)/(mlayer(1)-mlayer(0))
88
89      ! q_{1/2}
90      coefq(0)=volcapa*layer(1)/timestep
91        ! q_{k+1/2}
92        do ik=1,nsoil-1
93          coefq(ik)=volcapa*(layer(ik+1)-layer(ik))
94     &                 /timestep
95        enddo
96
97      do ig=1,ngrid
98        ! d_k
99        do ik=1,nsoil-1
100          coefd(ig,ik)=thermdiff(ig,ik)/(mlayer(ik)-mlayer(ik-1))
101        enddo
102       
103        ! alph_{N-1}
104        alph(ig,nsoil-1)=coefd(ig,nsoil-1)/
105     &                  (coefq(nsoil-1)+coefd(ig,nsoil-1))
106        ! alph_k
107        do ik=nsoil-2,1,-1
108          alph(ig,ik)=coefd(ig,ik)/(coefq(ik)+coefd(ig,ik+1)*
109     &                              (1.-alph(ig,ik+1))+coefd(ig,ik))
110        enddo
111
112        ! capcal
113! Cstar
114        capcal(ig)=volcapa*layer(1)+
115     &              (thermdiff(ig,1)/(mlayer(1)-mlayer(0)))*
116     &              (timestep*(1.-alph(ig,1)))
117! Cs
118        capcal(ig)=capcal(ig)/(1.+mu*(1.0-alph(ig,1))*
119     &                         thermdiff(ig,1)/mthermdiff(ig,0))
120      !write(*,*)'soil: ig=',ig,' capcal(ig)=',capcal(ig)
121      enddo ! of do ig=1,ngrid
122     
123      ! Additional checks: is the vertical discretization sufficient
124      ! to resolve diurnal and annual waves?
125      do ig=1,ngrid
126        ! extreme inertia for this column
127        inertia_min=minval(inertiedat(ig,:))
128        inertia_max=maxval(inertiedat(ig,:))
129        ! diurnal and annual skin depth
130        diurnal_skin=(inertia_min/volcapa)*sqrt(daysec/pi)
131        annual_skin=(inertia_max/volcapa)*sqrt(year_day*daysec/pi)
132        if (0.5*diurnal_skin<layer(1)) then
133        ! one should have the fist layer be at least half of diurnal skin depth
134          write(*,*) "soil Error: grid point ig=",ig
135          write(*,*) "            longitude=",longitude(ig)*(180./pi)
136          write(*,*) "             latitude=",latitude(ig)*(180./pi)
137          write(*,*) "  first soil layer depth ",layer(1)
138          write(*,*) "  not small enough for a diurnal skin depth of ",
139     &                diurnal_skin
140          write(*,*) " change soil layer distribution (comsoil_h.F90)"
141          stop
142        endif
143        if (2.*annual_skin>layer(nsoil)) then
144        ! one should have the full soil be at least twice the diurnal skin depth
145          write(*,*) "soil Error: grid point ig=",ig
146          write(*,*) "            longitude=",longitude(ig)*(180./pi)
147          write(*,*) "             latitude=",latitude(ig)*(180./pi)
148          write(*,*) "  total soil layer depth ",layer(nsoil)
149          write(*,*) "  not large enough for an annual skin depth of ",
150     &                annual_skin
151          write(*,*) " change soil layer distribution (comsoil_h.F90)"
152          stop
153        endif
154      enddo ! of do ig=1,ngrid
155     
156      else ! of if (firstcall)
157
158
159!  1. Compute soil temperatures
160! First layer:
161      do ig=1,ngrid
162        tsoil(ig,1)=(tsurf(ig)+mu*beta(ig,1)*
163     &                         thermdiff(ig,1)/mthermdiff(ig,0))/
164     &              (1.+mu*(1.0-alph(ig,1))*
165     &               thermdiff(ig,1)/mthermdiff(ig,0))
166      enddo
167! Other layers:
168      do ik=1,nsoil-1
169        do ig=1,ngrid
170          tsoil(ig,ik+1)=alph(ig,ik)*tsoil(ig,ik)+beta(ig,ik)
171        enddo
172      enddo
173     
174      endif! of if (firstcall)
175
176!  2. Compute beta coefficients (preprocessing for next time step)
177! Bottom layer, beta_{N-1}
178      do ig=1,ngrid
179        beta(ig,nsoil-1)=coefq(nsoil-1)*tsoil(ig,nsoil)
180     &                   /(coefq(nsoil-1)+coefd(ig,nsoil-1))
181      enddo
182! Other layers
183      do ik=nsoil-2,1,-1
184        do ig=1,ngrid
185          beta(ig,ik)=(coefq(ik)*tsoil(ig,ik+1)+
186     &                 coefd(ig,ik+1)*beta(ig,ik+1))/
187     &                 (coefq(ik)+coefd(ig,ik+1)*(1.0-alph(ig,ik+1))
188     &                  +coefd(ig,ik))
189        enddo
190      enddo
191
192
193!  3. Compute surface diffusive flux & calorific capacity
194      do ig=1,ngrid
195! Cstar
196!        capcal(ig)=volcapa(ig,1)*layer(ig,1)+
197!     &              (thermdiff(ig,1)/(mlayer(ig,1)-mlayer(ig,0)))*
198!     &              (timestep*(1.-alph(ig,1)))
199! Fstar
200
201!         print*,'this far in soil 1'
202!         print*,'thermdiff=',thermdiff(ig,1)
203!         print*,'mlayer=',mlayer
204!         print*,'beta=',beta(ig,1)
205!         print*,'alph=',alph(ig,1)
206!         print*,'tsoil=',tsoil(ig,1)
207
208        fluxgrd(ig)=(thermdiff(ig,1)/(mlayer(1)-mlayer(0)))*
209     &              (beta(ig,1)+(alph(ig,1)-1.0)*tsoil(ig,1))
210
211!        mu=mlayer(ig,0)/(mlayer(ig,1)-mlayer(ig,0))
212!        capcal(ig)=capcal(ig)/(1.+mu*(1.0-alph(ig,1))*
213!     &                         thermdiff(ig,1)/mthermdiff(ig,0))
214! Fs
215        fluxgrd(ig)=fluxgrd(ig)+(capcal(ig)/timestep)*
216     &              (tsoil(ig,1)*(1.+mu*(1.0-alph(ig,1))*
217     &                         thermdiff(ig,1)/mthermdiff(ig,0))
218     &               -tsurf(ig)-mu*beta(ig,1)*
219     &                          thermdiff(ig,1)/mthermdiff(ig,0))
220      enddo
221
222      end
223
Note: See TracBrowser for help on using the repository browser.