source: trunk/LMDZ.COMMON/libf/evolution/fast_subs_mars.f90 @ 3470

Last change on this file since 3470 was 3470, checked in by evos, 4 weeks ago

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