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

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

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