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

Last change on this file since 3495 was 3493, checked in by jbclement, 3 weeks ago

PEM:

  • Renaming of Norbert Schorghofer's subroutines with the prefix 'NS_';
  • Making the extension of all NS's subroutines as '.F90';
  • Deletion of the wrapper subroutine;
  • Making the initialization, variables management and arguments of the main subroutine for the dynamic computation of ice table to be more suitable.

JBC
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) :: 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
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
124  enddo  ! end of big loop
125  !$omp end do
126  !$omp end parallel
127end subroutine icelayer_mars
128
129
130
131subroutine ajsub_mars(typeT, albedo0, pfrost, nz, z, ti, rhocv, &
132     &     fracIR, fracDust, patm, avdrho, avdrhoP, avrho1, &
133     &     Tb, zdepthE, typeF, zdepthF, ypp, porefill, Tmean1, Tmean3, &
134     &     B, typeG, zdepthG, avrho1prescribed)
135!***********************************************************************
136!  A 1D thermal model that returns various averaged quantities
137!
138!  Tb<0 initializes temperatures
139!  Tb>0 initializes temperatures with Tb
140!***********************************************************************
141  use miscparameters
142  use thermalmodelparam_mars
143  use allinterfaces, except_this_one => ajsub_mars
144  implicit none
145  integer, intent(IN) :: nz, typeT
146!  real(8), intent(IN) :: latitude  ! in radians
147  real(8), intent(IN) :: albedo0, pfrost, z(NMAX)
148  real(8), intent(IN) :: ti(NMAX), rhocv(NMAX), fracIR, fracDust, patm
149  real(8), intent(IN) ::  porefill(nz)
150  real(8), intent(OUT) :: avdrho, avdrhoP  ! difference in annual mean vapor density
151  real(8), intent(OUT) :: avrho1  ! mean vapor density on surface
152  real(8), intent(INOUT) :: Tmean1
153  real(8), intent(INOUT) :: Tb(nz)
154  integer, intent(OUT) :: typeF  ! index of depth below which filling occurs
155  real(8), intent(OUT) :: zdepthE, zdepthF
156  real(8), intent(OUT) :: ypp(nz) ! (d rho/dt)/Diff
157  real(8), intent(OUT) :: Tmean3, zdepthG
158  real(8), intent(IN) :: B  ! just passing through
159  integer, intent(OUT) :: typeG
160  real(8), intent(IN), optional :: avrho1prescribed
161  real(8), parameter :: a = 1.52366 ! Mars semimajor axis in A.U.
162  integer nsteps, n, i, nm, typeP
163  real(8) tmax, time, Qn, Qnp1, tdays
164  real(8) marsR, marsLs, marsDec, HA
165  real(8) Tsurf, Tco2frost, albedo, Fsurf, m, dE, emiss, T(NMAX)
166  real(8) Told(nz), Fsurfold, Tsurfold, Tmean0, avrho2
167  real(8) rhosatav0, rhosatav(nz), rlow
168  real(8), external :: psv, tfrostco2
169 
170  tmax = EQUILTIME*solsperyear
171  nsteps = int(tmax/dt)     ! calculate total number of timesteps
172
173  Tco2frost = tfrostco2(patm)
174
175  !if (Tb<=0.) then  ! initialize
176     !Tmean0 = 210.15       ! black-body temperature of planet
177     !Tmean0 = (589.*(1.-albedo0)*cos(latitude)/(pi*emiss0*sigSB))**0.25 ! better estimate
178     !Tmean0 = Tmean0-5.
179     !write(34,*) '# initialized with temperature estimate at',latitude/d2r,'of',Tmean0,'K'
180     !T(1:nz) = Tmean0
181  !else
182     T(1:nz) = Tb
183     ! not so good when Fgeotherm is on
184  !endif
185 
186  albedo = albedo0
187  emiss = emiss0
188  do i=1,nz
189     if (T(i)<Tco2frost) T(i)=Tco2frost
190  enddo
191  Tsurf = T(1)
192  m=0.; Fsurf=0.
193
194  nm=0
195  avrho1=0.; avrho2=0.
196  Tmean1=0.; Tmean3=0.
197  rhosatav0 = 0.
198  rhosatav(:) = 0.
199
200  time=0.
201  !call generalorbit(0.d0,a,ecc,omega,eps,marsLs,marsDec,marsR)
202  !HA = 2.*pi*time            ! hour angle
203!  Qn=flux(marsR,marsDec,latitude,HA,albedo,fracir,fracdust,0.d0,0.d0)
204  !Qn = flux_mars77(marsR,marsDec,latitude,HA,albedo,fracir,fracdust)
205  !----loop over time steps
206  do n=0,nsteps-1
207     time = (n+1)*dt         !   time at n+1
208     tdays = time*(marsDay/earthDay) ! parenthesis may improve roundoff
209   !  call generalorbit(tdays,a,ecc,omega,eps,marsLs,marsDec,marsR)
210   !  HA = 2.*pi*mod(time,1.d0)  ! hour angle
211!     Qnp1=flux(marsR,marsDec,latitude,HA,albedo,fracir,fracdust,0.d0,0.d0)
212     !Qnp1 = flux_mars77(marsR,marsDec,latitude,HA,albedo,fracir,fracdust)
213     
214     Tsurfold = Tsurf
215     Fsurfold = Fsurf
216     Told(1:nz) = T(1:nz)
217     if (m<=0. .or. Tsurf>Tco2frost) then
218       ! call conductionQ(nz,z,dt*marsDay,Qn,Qnp1,T,ti,rhocv,emiss, &
219       !      &           Tsurf,Fgeotherm,Fsurf)
220     endif
221     if (Tsurf<Tco2frost .or. m>0.) then ! CO2 condensation
222        T(1:nz) = Told(1:nz)
223        !call conductionT(nz,z,dt*marsDay,T,Tsurfold,Tco2frost,ti, &
224             !&              rhocv,Fgeotherm,Fsurf)
225        Tsurf = Tco2frost
226    !    dE = (- Qn - Qnp1 + Fsurfold + Fsurf + &
227    !         &      emiss*sigSB*(Tsurfold**4+Tsurf**4)/2.
228        m = m + dt*marsDay*dE/Lco2frost
229     endif
230     if (Tsurf>Tco2frost .or. m<=0.) then
231        albedo = albedo0
232        emiss = emiss0
233     else
234        albedo = co2albedo
235        emiss = co2emiss
236     endif
237     !Qn=Qnp1
238     
239     if (time>=tmax-solsperyear) then
240        Tmean1 = Tmean1 + Tsurf
241        Tmean3 = Tmean3 + T(nz)
242        avrho1 = avrho1 + min(psv(Tsurf),pfrost)/Tsurf
243        rhosatav0 = rhosatav0 + psv(Tsurf)/Tsurf
244        do i=1,nz
245           rhosatav(i) = rhosatav(i) + psv(T(i))/T(i)
246        enddo
247        nm = nm+1
248     endif
249
250  enddo  ! end of time loop
251 
252  avrho1 = avrho1/nm
253  if (present(avrho1prescribed)) then
254     if (avrho1prescribed>=0.) avrho1=avrho1prescribed
255  endif
256  Tmean1 = Tmean1/nm; Tmean3 = Tmean3/nm
257  rhosatav0 = rhosatav0/nm; rhosatav(:)=rhosatav(:)/nm
258  if (typeT>0) then
259     avrho2 = rhosatav(typeT)
260  else
261     avrho2 = rhosatav(nz) ! no ice
262  endif
263  avdrho = avrho2-avrho1
264  typeP = -9
265  do i=1,nz
266     if (porefill(i)>0.) then
267        typeP = i  ! first point with ice
268        exit
269     endif
270  enddo
271  if (typeP>0) then
272     avdrhoP = rhosatav(typeP) - avrho1
273  else 
274     avdrhoP = -9999.
275  end if
276
277  zdepthE = equildepth(nz, z, rhosatav, rhosatav0, avrho1)
278
279  if (Fgeotherm>0.) then
280     Tb = Tmean1
281     typeG = 1   ! will be overwritten by depths_avmeth
282     rlow = 2*rhosatav(nz)-rhosatav(nz-1)
283  else
284     Tb = T(nz)
285     typeG = -9
286     rlow = rhosatav(nz-1)
287  endif
288  call depths_avmeth(typeT, nz, z, rhosatav(:), rhosatav0, rlow, avrho1,  &
289       & porefill(:), typeF, zdepthF, B, ypp(:), typeG, zdepthG)
290
291end subroutine ajsub_mars
292
293
294
295subroutine outputmoduleparameters
296  use miscparameters
297  use thermalmodelparam_mars
298  implicit none
299!  print *,'Parameters stored in modules'
300!  print *,'  Ice bulk density',icedensity,'kg/m^3'
301!  print *,'  dt=',dt,'Mars solar days'
302!  print *,'  Fgeotherm=',Fgeotherm,'W/m^2'
303!  write(6,'(2x,a27,1x,f5.3)') 'Emissivity of dry surface=',emiss0
304!  write(6,'(2x,a12,1x,f5.3,2x,a16,1x,f5.3)') 'CO2 albedo=',co2albedo,'CO2 emissivity=',co2emiss
305!  print *,'  Thermal model equilibration time',EQUILTIME,'Mars years'
306end subroutine outputmoduleparameters
307
308
309
Note: See TracBrowser for help on using the repository browser.