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

Last change on this file since 3567 was 3563, checked in by jbclement, 3 weeks ago

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