source: trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_mars.F90 @ 3538

Last change on this file since 3538 was 3538, checked in by jbclement, 13 days ago

PEM:
Fixing Norbert Schorghofer's subroutine which modified soil temperatures and was unwanted in the PEM loop.
JBC
File size: 11.3 KB
Line 
1MODULE fast_subs_mars
2
3implicit none
4
5  ! parameters for thermal model
6  ! they are only used in the subroutines below
7  real(8), parameter :: dt = 0.02  ! in units of Mars solar days
8  !real(8), parameter :: Fgeotherm = 0.
9  real(8), parameter :: Fgeotherm = 0 !0.028  ! [W/m^2]
10  real(8), parameter :: Lco2frost=6.0e5, co2albedo=0.60, co2emiss=1.
11  real(8), parameter :: emiss0 = 1.  ! emissivity of dry surface
12  integer, parameter :: EQUILTIME =15 ! [Mars years]
13
14  integer, parameter :: NMAX = 1000
15
16CONTAINS
17!*****************************************************
18! Subroutines for fast method
19! written by Norbert Schorghofer 2007-2011
20!*****************************************************
21
22
23subroutine icelayer_mars(bigstep,nz,NP,thIn,rhoc,z,porosity,pfrost, &
24     & Tb,zdepthF,zdepthE,porefill,Tmean1,Tmean3,zdepthG, &
25     & albedo,p0,icefrac,zdepthT,avrho1, &
26     & avrho1prescribed)
27!*************************************************************************
28! bigstep = time step [Earth years]
29! latitude  [degree]
30!*************************************************************************
31  use ice_table_mod, only: rho_ice
32  use fast_subs_univ, only: icechanges
33  !use omp_lib
34  implicit none
35  integer, intent(IN) :: nz, NP
36  real(8), intent(IN) :: bigstep
37  real(8), intent(IN) :: thIn(NP), rhoc(NP), z(NMAX), porosity, pfrost(NP)
38  real(8), intent(INOUT) :: porefill(nz,NP), zdepthF(NP), zdepthT(NP)
39  real(8), intent(INOUT) :: Tb(nz)
40  real(8), intent(OUT), dimension(NP) :: zdepthE, Tmean1, Tmean3, zdepthG
41  real(8), intent(IN), dimension(NP) ::  albedo, p0
42  real(8), intent(IN) ::  icefrac
43  real(8), intent(OUT) :: avrho1(NP)
44  real(8), intent(IN), optional :: avrho1prescribed(NP)  ! <0 means absent
45
46  integer k, typeF, typeG, typeT, j, jump, typeP
47  real(8) fracIR, fracDust, ti(NMAX), rhocv(NMAX)
48  real(8) Diff, ypp(nz), avdrho(NP), avdrhoP(NP), B, deltaz
49  real(8), SAVE :: avdrho_old(100), zdepth_old(100)  ! NP<=100
50  logical mode2
51
52  !$omp parallel &
53  !$omp    private(Diff,fracIR,fracDust,B,typeT,j,ti,rhocv,typeF,jump,typeG)
54  !$omp do
55  do k=1,NP   ! big loop
56
57     Diff = 4e-4*600./p0(k)
58     fracIR = 0.04*p0(k)/600.; fracDust = 0.02*p0(k)/600.
59     B = Diff*bigstep*86400.*365.24/(porosity*927.)
60     !B = Diff*bigstep*86400.*365.24/(porosity*rho_ice(Tb(),'h2o'))
61
62     typeT = -9
63     if (zdepthT(k)>=0. .and. zdepthT(k)<z(nz)) then
64        do j=1,nz
65           if (z(j)>zdepthT(k)) then ! ice
66              typeT = j  ! first point in ice
67              exit
68           endif
69        enddo
70     endif
71
72   !  call assignthermalproperties(nz,thIn(k),rhoc(k),ti,rhocv, &
73   !       & typeT,icefrac,porosity,porefill(:,k))
74
75     !----run thermal model
76     call ajsub_mars(typeT, albedo(k), pfrost(k), nz, z, &
77          &     ti, rhocv, fracIR, fracDust, p0(k), &
78          &     avdrho(k), avdrhoP(k), avrho1(k), Tb(k), zdepthE(k), typeF, &
79          &     zdepthF(k), ypp, porefill(:,k), Tmean1(k), Tmean3(k), B, &
80          &     typeG, zdepthG(k), avrho1prescribed(k))
81
82     if (typeF*zdepthF(k)<0.) stop 'error in icelayer_mars'
83     ! diagnose
84     if (zdepthT(k)>=0.) then
85        jump = 0
86        do j=1,nz
87           if (zdepth_old(k)<z(j) .and. zdepthT(k)>z(j)) jump = jump+1
88        enddo
89     else
90        jump = -9
91     endif
92     if (zdepthT(k)>=0. .and. avdrho(k)*avdrho_old(k)<0.) then
93        write(34,*) '# zdepth arrested'
94        if (jump>1) write(34,*) '# previous step was too large',jump
95     endif
96!     write(34,'(f8.3,1x,f6.2,1x,f11.5,2(1x,g11.4),1x,i3)') &
97!          &        bigstep,latitude(k),zdepthT(k),avdrho(k),avdrhoP(k),jump
98     zdepth_old(k) = zdepthT(k)
99     avdrho_old(k) = avdrho(k)
100
101!----mode 2 growth
102     typeP = -9;  mode2 = .FALSE.
103     do j=1,nz
104        if (porefill(j,k)>0.) then
105           typeP = j  ! first point with ice
106           exit
107        endif
108     enddo
109     if (typeT>0 .and. typeP>2 .and. zdepthE(k)>0.) then
110        if (porefill(typeP,k)>=1. .and. porefill(typeP-1,k)==0. .and. &
111             & zdepthE(k)<z(typeP) .and. &
112             & z(typeP)-zdepthE(k)>2*(z(typeP)-z(typeP-1))) then  ! trick that avoids oscillations
113           deltaz = -avdrhoP(k)/z(typeP)*18./8314.*B  ! conservation of mass
114           if (deltaz>z(typeP)-z(typeP-1)) then  ! also implies avdrhoP<0.
115              mode2 = .TRUE.
116           endif
117        endif
118     endif
119
120     !call icechanges_poreonly(nz,z,typeF,typeG,avdrhoP(k),ypp,B,porefill(:,k))
121     call icechanges(nz,z(:),typeF,avdrho(k),avdrhoP(k),ypp(:), &
122          & Diff,porosity,icefrac,bigstep,zdepthT(k),porefill(:,k),typeG)
123     if (typeP>2) then
124     if (mode2 .and. porefill(typeP,k)>=1. .and. porefill(typeP-1,k)==0.) then  ! nothing changed
125        porefill(typeP-1,k)=1.  ! paste a layer
126   !     write(34,*) '# mode 2 growth occurred',typeP,typeF,typeT
127     endif
128     endif
129
130  enddo  ! end of big loop
131  !$omp end do
132  !$omp end parallel
133end subroutine icelayer_mars
134
135
136
137subroutine ajsub_mars(typeT, albedo0, pfrost, nz, z, ti, rhocv, &
138     &     fracIR, fracDust, patm, avdrho, avdrhoP, avrho1, &
139     &     Tb, zdepthE, typeF, zdepthF, ypp, porefill, Tmean1, Tmean3, &
140     &     B, typeG, zdepthG, avrho1prescribed)
141!***********************************************************************
142!  A 1D thermal model that returns various averaged quantities
143!
144!  Tb<0 initializes temperatures
145!  Tb>0 initializes temperatures with Tb
146!***********************************************************************
147  use fast_subs_univ, only: depths_avmeth, equildepth
148  use constants_marspem_mod, only: sols_per_my, sec_per_sol
149  implicit none
150  integer, intent(IN) :: nz, typeT
151!  real(8), intent(IN) :: latitude  ! in radians
152  real(8), intent(IN) :: albedo0, pfrost, z(NMAX)
153  real(8), intent(IN) :: ti(NMAX), rhocv(NMAX), fracIR, fracDust, patm
154  real(8), intent(IN) ::  porefill(nz)
155  real(8), intent(OUT) :: avdrho, avdrhoP  ! difference in annual mean vapor density
156  real(8), intent(OUT) :: avrho1  ! mean vapor density on surface
157  real(8), intent(INOUT) :: Tmean1
158  real(8), intent(INOUT) :: Tb(nz)
159  integer, intent(OUT) :: typeF  ! index of depth below which filling occurs
160  real(8), intent(OUT) :: zdepthE, zdepthF
161  real(8), intent(OUT) :: ypp(nz) ! (d rho/dt)/Diff
162  real(8), intent(OUT) :: Tmean3, zdepthG
163  real(8), intent(IN) :: B  ! just passing through
164  integer, intent(OUT) :: typeG
165  real(8), intent(IN), optional :: avrho1prescribed
166  real(8), parameter :: a = 1.52366 ! Mars semimajor axis in A.U.
167  integer nsteps, n, i, nm, typeP
168  real(8) tmax, time, Qn, Qnp1, tdays
169  real(8) marsR, marsLs, marsDec, HA
170  real(8) Tsurf, Tco2frost, albedo, Fsurf, m, dE, emiss, T(NMAX)
171  real(8) Told(nz), Fsurfold, Tsurfold, Tmean0, avrho2
172  real(8) rhosatav0, rhosatav(nz), rlow
173
174  tmax = EQUILTIME*sols_per_my
175  nsteps = int(tmax/dt)     ! calculate total number of timesteps
176
177  Tco2frost = tfrostco2(patm)
178
179  !if (Tb<=0.) then  ! initialize
180     !Tmean0 = 210.15       ! black-body temperature of planet
181     !Tmean0 = (589.*(1.-albedo0)*cos(latitude)/(pi*emiss0*sigSB))**0.25 ! better estimate
182     !Tmean0 = Tmean0-5.
183     !write(34,*) '# initialized with temperature estimate at',latitude/d2r,'of',Tmean0,'K'
184     !T(1:nz) = Tmean0
185  !else
186     T(1:nz) = Tb
187     ! not so good when Fgeotherm is on
188  !endif
189
190  albedo = albedo0
191  emiss = emiss0
192  do i=1,nz
193     if (T(i)<Tco2frost) T(i)=Tco2frost
194  enddo
195  Tsurf = T(1)
196  m=0.; Fsurf=0.
197
198  nm=0
199  avrho1=0.; avrho2=0.
200  Tmean1=0.; Tmean3=0.
201  rhosatav0 = 0.
202  rhosatav(:) = 0.
203
204  time=0.
205  !call generalorbit(0.d0,a,ecc,omega,eps,marsLs,marsDec,marsR)
206  !HA = 2.*pi*time            ! hour angle
207!  Qn=flux(marsR,marsDec,latitude,HA,albedo,fracir,fracdust,0.d0,0.d0)
208  !Qn = flux_mars77(marsR,marsDec,latitude,HA,albedo,fracir,fracdust)
209  !----loop over time steps
210  do n=0,nsteps-1
211     time = (n+1)*dt         !   time at n+1
212   !  tdays = time*(sec_per_sol/earthDay) ! parenthesis may improve roundoff
213   !  call generalorbit(tdays,a,ecc,omega,eps,marsLs,marsDec,marsR)
214   !  HA = 2.*pi*mod(time,1.d0)  ! hour angle
215!     Qnp1=flux(marsR,marsDec,latitude,HA,albedo,fracir,fracdust,0.d0,0.d0)
216     !Qnp1 = flux_mars77(marsR,marsDec,latitude,HA,albedo,fracir,fracdust)
217
218     Tsurfold = Tsurf
219     Fsurfold = Fsurf
220     Told(1:nz) = T(1:nz)
221     if (m<=0. .or. Tsurf>Tco2frost) then
222       ! call conductionQ(nz,z,dt*sec_per_sol,Qn,Qnp1,T,ti,rhocv,emiss, &
223       !      &           Tsurf,Fgeotherm,Fsurf)
224     endif
225     if (Tsurf<Tco2frost .or. m>0.) then ! CO2 condensation
226        T(1:nz) = Told(1:nz)
227        !call conductionT(nz,z,dt*sec_per_sol,T,Tsurfold,Tco2frost,ti, &
228             !&              rhocv,Fgeotherm,Fsurf)
229        Tsurf = Tco2frost
230    !    dE = (- Qn - Qnp1 + Fsurfold + Fsurf + &
231    !         &      emiss*sigSB*(Tsurfold**4+Tsurf**4)/2.
232        m = m + dt*sec_per_sol*dE/Lco2frost
233     endif
234     if (Tsurf>Tco2frost .or. m<=0.) then
235        albedo = albedo0
236        emiss = emiss0
237     else
238        albedo = co2albedo
239        emiss = co2emiss
240     endif
241     !Qn=Qnp1
242
243     if (time>=tmax-sols_per_my) then
244        Tmean1 = Tmean1 + Tsurf
245        Tmean3 = Tmean3 + T(nz)
246        avrho1 = avrho1 + min(psv(Tsurf),pfrost)/Tsurf
247        rhosatav0 = rhosatav0 + psv(Tsurf)/Tsurf
248        do i=1,nz
249           rhosatav(i) = rhosatav(i) + psv(T(i))/T(i)
250        enddo
251        nm = nm+1
252     endif
253
254  enddo  ! end of time loop
255
256  avrho1 = avrho1/nm
257  if (present(avrho1prescribed)) then
258     if (avrho1prescribed>=0.) avrho1=avrho1prescribed
259  endif
260  Tmean1 = Tmean1/nm; Tmean3 = Tmean3/nm
261  rhosatav0 = rhosatav0/nm; rhosatav(:)=rhosatav(:)/nm
262  if (typeT>0) then
263     avrho2 = rhosatav(typeT)
264  else
265     avrho2 = rhosatav(nz) ! no ice
266  endif
267  avdrho = avrho2-avrho1
268  typeP = -9
269  do i=1,nz
270     if (porefill(i)>0.) then
271        typeP = i  ! first point with ice
272        exit
273     endif
274  enddo
275  if (typeP>0) then
276     avdrhoP = rhosatav(typeP) - avrho1
277  else
278     avdrhoP = -9999.
279  end if
280
281  zdepthE = equildepth(nz, z, rhosatav, rhosatav0, avrho1)
282
283  if (Fgeotherm>0.) then
284     !Tb = Tmean1
285     typeG = 1   ! will be overwritten by depths_avmeth
286     rlow = 2*rhosatav(nz)-rhosatav(nz-1)
287  else
288     !Tb = T(nz)
289     typeG = -9
290     rlow = rhosatav(nz-1)
291  endif
292  call depths_avmeth(typeT, nz, z, rhosatav(:), rhosatav0, rlow, avrho1,  &
293       & porefill(:), typeF, zdepthF, B, ypp(:), typeG, zdepthG)
294
295end subroutine ajsub_mars
296
297
298
299real*8 function tfrostco2(p)
300!     the inverse of function psvco2
301!     input is pressure in Pascal, output is temperature in Kelvin
302implicit none
303real*8 p
304tfrostco2 = 3182.48/(23.3494+log(100./p))
305end function
306
307!=======================================================================
308
309      real*8 function psv(T)
310!     saturation vapor pressure of H2O ice [Pascal]
311!     input is temperature [Kelvin]
312      implicit none
313      real*8 T
314
315!-----parametrization 1
316!      real*8 DHmelt,DHvap,DHsub,R,pt,Tt,C
317!      parameter (DHmelt=6008.,DHvap=45050.)
318!      parameter (DHsub=DHmelt+DHvap) ! sublimation enthalpy [J/mol]
319!      parameter (R=8.314,pt=6.11e2,Tt=273.16)
320!      C = (DHsub/R)*(1./T - 1./Tt)
321!      psv = pt*exp(-C)
322
323!-----parametrization 2
324!     eq. (2) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)
325!     differs from parametrization 1 by only 0.1%
326      real*8 A,B
327      parameter (A=-6143.7, B=28.9074)
328      psv = exp(A/T+B)  ! Clapeyron
329
330!-----parametrization 3
331!     eq. (7) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005)
332!     psv = exp(9.550426 - 5723.265/T + 3.53068*log(T) - 0.00728332*T)
333
334      end function psv
335
336END MODULE fast_subs_mars
Note: See TracBrowser for help on using the repository browser.