Index: trunk/LMDZ.COMMON/libf/evolution/NS_conductionQ.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/NS_conductionQ.F90 (revision 3493) +++ trunk/LMDZ.COMMON/libf/evolution/NS_conductionQ.F90 (revision 3493) @@ -0,0 +1,113 @@ +subroutine conductionQ(nz,z,dt,Qn,Qnp1,T,ti,rhoc,emiss,Tsurf, & + & Fgeotherm,Fsurf) +!*********************************************************************** +! conductionQ: program to calculate the diffusion of temperature +! into the ground and thermal emission at the surface +! with variable thermal properties on irregular grid +! Crank-Nicolson scheme, flux conservative +! uses Samar's radiation formula +! Eqn: rhoc*T_t = (k*T_z)_z +! BC (z=0): Q(t) + kT_z = em*sig*T^4 +! BC (z=L): heat flux = Fgeotherm +! +! nz = number of grid points +! dt = time step +! Qn,Qnp1 = net solar insolation at time steps n and n+1 [W/m^2] +! T = vertical temperature profile [K] (in- and output) +! ti = thermal inertia [J m^-2 K^-1 s^-1/2] VECTOR +! rhoc = rho*c VECTOR where rho=density [kg/m^3] and +! c=specific heat [J K^-1 kg^-1] +! ti and rhoc are not allowed to vary in the layers immediately +! adjacent to the surface or the bottom +! emiss = emissivity +! Tsurf = surface temperature [K] (in- and output) +! Fgeotherm = geothermal heat flux at bottom boundary [W/m^2] +! Fsurf = heat flux at surface [W/m^2] (output) +! +! Grid: surface is at z=0 +! z(0)=0, z(2)=3*z(1), i.e., the top layer has half the width +! T(1) is at z(1); ...; T(i) is at z(i) +! k(i), rhoc(i), ti(i) are midway between z(i-1) and z(i) +! +! originally written by Samar Khatiwala, 2001 +! extended to variable thermal properties +! and irregular grid by Norbert Schorghofer +! added predictor-corrector 9/2019 +!*********************************************************************** + + implicit none + real*8, parameter :: sigSB=5.6704d-8 + + integer, intent(IN) :: nz + real*8, intent(IN) :: z(nz), dt, Qn, Qnp1, ti(nz),rhoc(nz) + real*8, intent(IN) :: emiss, Fgeotherm + real*8, intent(INOUT) :: T(nz), Tsurf + real*8, intent(OUT) :: Fsurf + integer i, iter + real*8 k(nz), k1, alpha(nz), gamma(nz), Tr + real*8 a(nz), b(nz), c(nz), r(nz), Told(nz) + real*8 arad, brad, ann, annp1, bn, buf, dz, beta + + ! set some constants + k(:) = ti(:)**2/rhoc(:) ! thermal conductivity + dz = 2.*z(1) + beta = dt/rhoc(1)/(2.*dz**2) ! assumes rhoc(0)=rhoc(1) + alpha(1) = beta*k(2) + gamma(1) = beta*k(1) + do i=2,nz-1 + buf = dt/(z(i+1)-z(i-1)) + alpha(i) = 2.*k(i+1)*buf/(rhoc(i)+rhoc(i+1))/(z(i+1)-z(i)) + gamma(i) = 2.*k(i)*buf/(rhoc(i)+rhoc(i+1))/(z(i)-z(i-1)) + enddo + buf = dt/(z(nz)-z(nz-1))**2 + gamma(nz) = k(nz)*buf/(2*rhoc(nz)) ! assumes rhoc(nz+1)=rhoc(nz) + + k1 = k(1)/dz + + ! elements of tridiagonal matrix + a(:) = -gamma(:) ! a(1) is not used + b(:) = 1. + alpha(:) + gamma(:) ! b(1) has to be reset at every timestep + c(:) = -alpha(:) ! c(nz) is not used + b(nz) = 1. + gamma(nz) + + Tr = Tsurf ! 'reference' temperature + iter = 0 + Told(:) = T(:) +30 continue ! in rare case of iterative correction + + ! Emission + arad = -3.*emiss*sigSB*Tr**4 + brad = 2.*emiss*sigSB*Tr**3 + ann = (Qn-arad)/(k1+brad) + annp1 = (Qnp1-arad)/(k1+brad) + bn = (k1-brad)/(k1+brad) + b(1) = 1. + alpha(1) + gamma(1) - gamma(1)*bn + + ! Set RHS + r(1) = gamma(1)*(annp1+ann) + & + & (1.-alpha(1)-gamma(1)+gamma(1)*bn)*T(1) + alpha(1)*T(2) + do concurrent (i=2:nz-1) + r(i) = gamma(i)*T(i-1) + (1.-alpha(i)-gamma(i))*T(i)+ alpha(i)*T(i+1) + enddo + r(nz) = gamma(nz)*T(nz-1) + (1.-gamma(nz))*T(nz) & + & + dt/rhoc(nz)*Fgeotherm/(z(nz)-z(nz-1)) ! assumes rhoc(nz+1)=rhoc(nz) + + ! Solve for T at n+1 + call tridag(a,b,c,r,T,nz) ! update by tridiagonal inversion + + Tsurf = 0.5*(annp1 + bn*T(1) + T(1)) ! (T0+T1)/2 + + ! iterative predictor-corrector + if ((Tsurf > 1.2*Tr .or. Tsurf < 0.8*Tr) .and. iter<10) then ! linearization error expected + ! redo until Tr is within 20% of new surface temperature + ! (under most circumstances, the 20% threshold is never exceeded) + iter = iter+1 + Tr = sqrt(Tr*Tsurf) ! linearize around an intermediate temperature + T(:) = Told(:) + goto 30 + endif + !if (iter>=5) print *,'consider taking shorter time steps',iter + !if (iter>=10) print *,'warning: too many iterations' + + Fsurf = - k(1)*(T(1)-Tsurf)/z(1) ! heat flux into surface +end subroutine conductionQ Index: trunk/LMDZ.COMMON/libf/evolution/NS_conductionT.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/NS_conductionT.F90 (revision 3493) +++ trunk/LMDZ.COMMON/libf/evolution/NS_conductionT.F90 (revision 3493) @@ -0,0 +1,70 @@ +subroutine conductionT(nz,z,dt,T,Tsurf,Tsurfp1,ti,rhoc,Fgeotherm,Fsurf) +!*********************************************************************** +! conductionT: program to calculate the diffusion of temperature +! into the ground with prescribed surface temperature +! and variable thermal properties on irregular grid +! Crank-Nicholson scheme, flux conservative +! +! Eqn: rhoc*T_t = (k*T_z)_z +! BC (z=0): T=T(t) +! BC (z=L): heat flux = Fgeotherm +! +! nz = number of grid points +! dt = time step +! T = vertical temperature profile [K] (in- and output) +! Tsurf, Tsurfp1 = surface temperatures at times n and n+1 +! ti = thermal inertia [J m^-2 K^-1 s^-1/2] VECTOR +! rhoc = rho*c heat capacity per volume [J m^-3 K^-1] VECTOR +! ti and rhoc are not allowed to vary in the layers immediately +! adjacent to the surface or the bottom +! Fgeotherm = geothermal heat flux at bottom boundary [W/m^2] +! Fsurf = heat flux at surface [W/m^2] (output) +! +! Grid: surface is at z=0 +! T(1) is at z(1); ...; T(i) is at z(i) +! k(i) is midway between z(i-1) and z(i) +! rhoc(i) is midway between z(i-1) and z(i) +!*********************************************************************** + implicit none + + integer, intent(IN) :: nz + real*8, intent(IN) :: z(nz), dt, Tsurf, Tsurfp1, ti(nz), rhoc(nz) + real*8, intent(IN) :: Fgeotherm + real*8, intent(INOUT) :: T(nz) + real*8, intent(OUT) :: Fsurf + integer i + real*8 alpha(nz), k(nz), gamma(nz), buf + real*8 a(nz), b(nz), c(nz), r(nz) + + ! set some constants + k(:) = ti(:)**2/rhoc(:) ! thermal conductivity + alpha(1) = k(2)*dt/rhoc(1)/(z(2)-z(1))/z(2) + gamma(1) = k(1)*dt/rhoc(1)/z(1)/z(2) + do i=2,nz-1 + buf = dt/(z(i+1)-z(i-1)) + alpha(i) = k(i+1)*buf*2./(rhoc(i)+rhoc(i+1))/(z(i+1)-z(i)) + gamma(i) = k(i)*buf*2./(rhoc(i)+rhoc(i+1))/(z(i)-z(i-1)) + enddo + buf = dt/(z(nz)-z(nz-1))**2 + gamma(nz) = k(nz)*buf/(rhoc(nz)+rhoc(nz)) ! assumes rhoc(nz+1)=rhoc(nz) + + ! elements of tridiagonal matrix + a(:) = -gamma(:) ! a(1) is not used + b(:) = 1. + alpha(:) + gamma(:) + c(:) = -alpha(:) ! c(nz) is not used + b(nz) = 1. + gamma(nz) + + ! Set RHS + r(1)= alpha(1)*T(2) + (1.-alpha(1)-gamma(1))*T(1) + gamma(1)*(Tsurf+Tsurfp1) + do concurrent (i=2:nz-1) + r(i) = gamma(i)*T(i-1) + (1.-alpha(i)-gamma(i))*T(i) + alpha(i)*T(i+1) + enddo + r(nz) = gamma(nz)*T(nz-1) + (1.-gamma(nz))*T(nz) + & + & dt/rhoc(nz)*Fgeotherm/(z(nz)-z(nz-1)) ! assumes rhoc(nz+1)=rhoc(nz) + + ! Solve for T at n+1 + call tridag(a,b,c,r,T,nz) ! update by tridiagonal inversion + + Fsurf = -k(1)*(T(1)-Tsurfp1)/z(1) ! heat flux into surface + +end subroutine conductionT Index: trunk/LMDZ.COMMON/libf/evolution/NS_derivs.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/NS_derivs.F90 (revision 3493) +++ trunk/LMDZ.COMMON/libf/evolution/NS_derivs.F90 (revision 3493) @@ -0,0 +1,115 @@ + +subroutine deriv1(z,nz,y,y0,yNp1,yp) + ! first derivative of a function y(z) on irregular grid + ! upper b.c.: y(0)=y0 + ! lower b.c.: yp =0 + implicit none + integer, intent(IN) :: nz + real(8), intent(IN) :: z(nz),y(nz),y0,yNp1 + real(8), intent(OUT) :: yp(nz) + integer j + real(8) hm,hp,c1,c2,c3 + + hp = z(2)-z(1) + !hm = z(1) + c1 = z(1)/(hp*z(2)) + c2 = 1/z(1) - 1/(z(2)-z(1)) + c3 = -hp/(z(1)*z(2)) + yp(1) = c1*y(2) + c2*y(1) + c3*y0 + do j=2,nz-1 + hp = z(j+1)-z(j) + hm = z(j)-z(j-1) + c1 = +hm/(hp*(z(j+1)-z(j-1))) + c2 = 1/hm - 1/hp + c3 = -hp/(hm*(z(j+1)-z(j-1))) + yp(j) = c1*y(j+1) + c2*y(j) + c3*y(j-1) + enddo + yp(nz) = (yNp1 - y(nz-1))/(2.*(z(nz)-z(nz-1))) +end subroutine deriv1 + + + +subroutine deriv2_full(z,nz,a,b,a0,b0,bNp1,yp2) + ! second derivative (a*b_z)_z on irregular grid + ! upper b.c.: a(0)=a0, b(0)=b0 + ! 2nd order, without cross-terms + implicit none + integer, intent(IN) :: nz + real(8), intent(IN) :: z(nz),a(nz),b(nz),a0,b0,bNp1 + real(8), intent(OUT) :: yp2(nz) + integer j + real(8) hm,hp,c1,c2,c3 + + c1 = 1./(z(1)*z(2)) + c2 = 1./((z(2)-z(1))*z(1)) + c3 = 1./((z(2)-z(1))*z(2)) + yp2(1) = -c2*a(1)*b(1) & + & -c1*a0*b(1) + c1*(a(1)+a0)*b0 & + & -c3*a(2)*b(1) + c3*(a(1)+a(2))*b(2) + do j=2,nz-1 + hp = z(j+1)-z(j) + hm = z(j)-z(j-1) + c1 = 1./(hm*(z(j+1)-z(j-1))) + c2 = 1./(hp*hm) + c3 = 1./(hp*(z(j+1)-z(j-1))) + yp2(j) = -c2*a(j)*b(j) & + & -c1*a(j-1)*b(j) + c1*(a(j)+a(j-1))*b(j-1) & + & -c3*a(j+1)*b(j) + c3*(a(j)+a(j+1))*b(j+1) + enddo + + ! lower bc: b_z = 0 + !yp(nz)= 2.*a(nz)*(b(nz-1)-b(nz))/(z(nz)-z(nz-1))**2 + !more general lower bc + yp2(nz) = a(nz)*(bNp1 - 2*b(nz) + b(nz-1))/(z(nz)-z(nz-1))**2 +end subroutine deriv2_full + + + +subroutine deriv2_simple(z,nz,y,y0,yNp1,yp2) + ! second derivative y_zz on irregular grid + ! boundary conditions: y(0)=y0, y(nz+1)=yNp1 + implicit none + integer, intent(IN) :: nz + real(8), intent(IN) :: z(nz),y(nz),y0,yNp1 + real(8), intent(OUT) :: yp2(nz) + integer j + real(8) hm,hp,c1,c2,c3 + + c1 = +2./((z(2)-z(1))*z(2)) + c2 = -2./((z(2)-z(1))*z(1)) + c3 = +2./(z(1)*z(2)) + yp2(1) = c1*y(2) + c2*y(1) + c3*y0 + + do j=2,nz-1 + hp = z(j+1)-z(j) + hm = z(j)-z(j-1) + c1 = +2./(hp*(z(j+1)-z(j-1))) + c2 = -2./(hp*hm) + c3 = +2./(hm*(z(j+1)-z(j-1))) + yp2(j) = c1*y(j+1) + c2*y(j) + c3*y(j-1) + enddo + + yp2(nz) = (yNp1 - 2*y(nz) + y(nz-1))/(z(nz)-z(nz-1))**2 +end subroutine deriv2_simple + + + +real(8) function deriv1_onesided(j,z,nz,y) + ! first derivative of function y(z) at z(j) + ! one-sided derivative on irregular grid + implicit none + integer, intent(IN) :: nz,j + real(8), intent(IN) :: z(nz),y(nz) + real(8) h1,h2,c1,c2,c3 + if (j<1 .or. j>nz-2) then + deriv1_onesided = -9999. + else + h1 = z(j+1)-z(j) + h2 = z(j+2)-z(j+1) + c1 = -(2*h1+h2)/(h1*(h1+h2)) + c2 = (h1+h2)/(h1*h2) + c3 = -h1/(h2*(h1+h2)) + deriv1_onesided = c1*y(j) + c2*y(j+1) + c3*y(j+2) + endif +end function deriv1_onesided + Index: trunk/LMDZ.COMMON/libf/evolution/NS_dyn_ss_ice_m.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/NS_dyn_ss_ice_m.F90 (revision 3493) +++ trunk/LMDZ.COMMON/libf/evolution/NS_dyn_ss_ice_m.F90 (revision 3493) @@ -0,0 +1,177 @@ + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! +!!! Purpose: Retreat and growth of subsurface ice on Mars +!!! orbital elements remain constant +!!! +!!! +!!! Author: EV, updated NS MSIM dynamical program for the PEM +!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + +SUBROUTINE dyn_ss_ice_m(ssi_depth_in,T1,Tb,nz,thIn,p0,pfrost,porefill_in,porefill,ssi_depth) + +!*********************************************************************** +! Retreat and growth of subsurface ice on Mars +! orbital elements remain constant +!*********************************************************************** + use miscparameters, only : pi, d2r, NMAX, marsDay, solsperyear + !use allinterfaces + implicit none + integer, parameter :: NP=1 ! # of sites + integer nz, i, k, iloop + real(8) zmax, delta, z(NMAX), icetime, porosity, icefrac + real(8), dimension(NP) :: albedo, thIn, rhoc + real(8), dimension(NP) :: pfrost, p0 + real(8) newti, stretch, newrhoc, ecc, omega, eps, timestep + real(8) ssi_depth_in, ssi_depth, T1 + real(8), dimension(NP) :: zdepthF, zdepthE, zdepthT, zdepthG + real(8), dimension(NMAX,NP) :: porefill, porefill_in + real(8), dimension(nz) :: Tb + real(8), dimension(NP) :: Tmean1, Tmean3, avrho1 + real(8) tmax, tlast, avrho1prescribed(NP), l1 + real(8), external :: smartzfac + + !if (iargc() /= 1) then + ! stop 'USAGE: icages ext' + !endif + !call getarg( 1, ext ) + + if (NP>100) stop 'subroutine icelayer_mars cannot handle this many sites' + + ! parameters that never ever change + porosity = 0.4d0 ! porosity of till + !rhoc(:) = 1500.*800. ! will be overwritten + icefrac = 0.98 + tmax = 1 + tlast = 0. + avrho1prescribed(:) = pfrost/T1 ! <0 means absent + albedo=0.23 + !avrho1prescribed(:) = 0.16/200. ! units are Pa/K + + !open(unit=21,file='lats.'//ext,action='read',status='old',iostat=ierr) + !if (ierr /= 0) then + ! print *,'File lats.'//ext,'not found' + ! stop + !endif + do k=1,NP + !read(21,*) latitude(k),albedo(k),thIn(k),htopo(k) + ! empirical relation from Mellon & Jakosky + rhoc(k) = 800.*(150.+100.*sqrt(34.2+0.714*thIn(k))) + enddo + !close(21) + + ! set eternal grid + zmax = 25. + !zfac = smartzfac(nz,zmax,6,0.032d0) + !call setgrid(nz,z,zmax,zfac) + l1=2.e-4 + do iloop=0,nz + z(iloop) = l1*(1+iloop**2.9*(1-exp(-real(iloop)/20.))) + enddo + + + !open(unit=30,file='z.'//ext,action='write',status='unknown') + !write(30,'(999(f8.5,1x))') z(1:nz) + !close(30) + + !ecc = ecc_in; eps = obl_in*d2r; omega = Lp_in*d2r ! today + ! total atmospheric pressure + !p0(:) = 600. + ! presently 520 Pa at zero elevation (Smith & Zuber, 1998) + ! do k=1,NP + ! p0(k)=520*exp(-htopo(k)/10800.) + ! enddo + timestep = 1 ! must be integer fraction of 1 ka + icetime = -tmax-timestep ! earth years + + ! initializations + !Tb = -9999. + zdepthF(:) = -9999. + + !zdepthT(1:NP) = -9999. ! reset again below + ! zdepthT(1:NP) = 0. + + ! print *,'RUNNING MARS_FAST' + ! print *,'Global model parameters:' + ! print *,'nz=',nz,' zfac=',zfac,'zmax=',zmax + ! print *,'porosity=',porosity + ! print *,'starting at time',icetime,'years' + ! print *,'time step=',timestep,'years' + ! print *,'eps=',eps/d2r,'ecc=',ecc,'omega=',omega/d2r + ! print *,'number of sites=',NP + ! print *,'Site specific parameters:' + do k=1,NP + if (NP>1) print *,' Site ',k + ! print *,' latitude (deg)',latitude(k),' rho*c (J/m^3/K)',rhoc(k),' thIn=',thIn(k) + ! print *,' total pressure=',p0(k),'partial pressure=',pfrost(k) + delta = thIn(k)/rhoc(k)*sqrt(marsDay/pi) + ! print *,' skin depths (m)',delta,delta*sqrt(solsperyear) + call soilthprop(porosity,1.d0,rhoc(k),thIn(k),1,newrhoc,newti,icefrac) + stretch = (newti/thIn(k))*(rhoc(k)/newrhoc) + do i=1,nz + if (z(i)=tlast) exit + enddo + + ! close(34) + ! close(36); close(37) + +!501 format (f10.0,2x,f7.3,2x,f10.4,2(2x,f6.2),2x,f9.3,2x,g11.4) + +end subroutine dyn_ss_ice_m + + + Index: trunk/LMDZ.COMMON/libf/evolution/NS_fast_modules.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/NS_fast_modules.F90 (revision 3493) +++ trunk/LMDZ.COMMON/libf/evolution/NS_fast_modules.F90 (revision 3493) @@ -0,0 +1,265 @@ + +module allinterfaces + ! interfaces from Fortran 90 subroutines and functions + ! comments have been stripped + + !begin fast_subs_mars.f90 + + interface + subroutine icelayer_mars(bigstep,nz,NP,thIn,rhoc,z,porosity,pfrost, & + & Tb,zdepthF,zdepthE,porefill,Tmean1,Tmean3,zdepthG, & + & albedo,p0,icefrac,zdepthT,avrho1, & + & avrho1prescribed) + use miscparameters, only : d2r, NMAX, icedensity + implicit none + integer, intent(IN) :: nz, NP + real(8), intent(IN) :: bigstep + real(8), intent(IN), dimension(NP) :: thIn, rhoc + real(8), intent(IN) :: z(NMAX) + real(8), intent(IN) :: porosity, pfrost(NP) + real(8), intent(INOUT) :: Tb(NP), porefill(nz,NP), zdepthF(NP), zdepthT(NP) + real(8), intent(OUT), dimension(NP) :: zdepthE, Tmean1, Tmean3, zdepthG + real(8), intent(IN), dimension(NP) :: albedo, p0 + real(8), intent(IN) :: icefrac + real(8), intent(OUT) :: avrho1(NP) + real(8), intent(IN), optional :: avrho1prescribed(NP) + end subroutine icelayer_mars + end interface + + interface + subroutine ajsub_mars(typeT, albedo0, pfrost, nz, z, ti, rhocv, & + & fracIR, fracDust, p0, avdrho, avdrhoP, avrho1, & + & Tb, zdepthE, typeF, zdepthF, ypp, porefill, Tmean1, Tmean3, & + & B, typeG, zdepthG, avrho1prescribed) + use miscparameters + implicit none + integer, intent(IN) :: nz, typeT + real(8), intent(IN) :: albedo0, pfrost, z(NMAX) + real(8), intent(IN) :: ti(NMAX), rhocv(NMAX), fracIR, fracDust, p0 + real(8), intent(IN) :: porefill(nz) + real(8), intent(OUT) :: avdrho, avdrhoP + real(8), intent(OUT) :: avrho1 + real(8), intent(INOUT) :: Tb, Tmean1 + integer, intent(OUT) :: typeF + real(8), intent(OUT) :: zdepthE, zdepthF + real(8), intent(OUT) :: ypp(nz) + real(8), intent(OUT) :: Tmean3, zdepthG + real(8), intent(IN) :: B + integer, intent(OUT) :: typeG + real(8), intent(IN), optional :: avrho1prescribed + end subroutine ajsub_mars + end interface + + !end of fast_subs_mars.f90 + !begin fast_subs_univ.f90 + + interface + pure function zint(y1,y2,z1,z2) + implicit none + real(8), intent(IN) :: y1,y2,z1,z2 + real(8) zint + end function zint + end interface + + interface + pure function equildepth(nz, z, rhosatav, rhosatav0, avrho1) + implicit none + integer, intent(IN) :: nz + real(8), intent(IN) :: z(nz), rhosatav(nz) + real(8), intent(IN) :: rhosatav0, avrho1 + real(8) zint, equildepth + external zint + end function equildepth + end interface + + interface + subroutine depths_avmeth(typeT, nz, z, rhosatav, rhosatav0, rlow, avrho1, & + & porefill, typeF, zdepthF, B, ypp, typeG, zdepthG) + use miscparameters, only : icedensity + implicit none + integer, intent(IN) :: nz, typeT + real(8), intent(IN), dimension(nz) :: z, rhosatav, porefill + real(8), intent(IN) :: rhosatav0, rlow, avrho1 + integer, intent(OUT) :: typeF + real(8), intent(INOUT) :: zdepthF + real(8), intent(IN) :: B + real(8), intent(OUT) :: ypp(nz), zdepthG + integer, intent(INOUT) :: typeG + real(8), external :: zint, deriv1_onesided, colint + end subroutine depths_avmeth + end interface + + interface + pure function constriction(porefill) + implicit none + real(8), intent(IN) :: porefill + real(8) constriction + end function constriction + end interface + + interface + subroutine assignthermalproperties(nz,thIn,rhoc,ti,rhocv, & + & typeT,icefrac,porosity,porefill) + implicit none + integer, intent(IN) :: nz + integer, intent(IN), optional :: typeT + real(8), intent(IN), optional :: icefrac + real(8), intent(IN) :: thIn, rhoc + real(8), intent(IN), optional :: porosity, porefill(nz) + real(8), intent(OUT) :: ti(nz), rhocv(nz) + end subroutine assignthermalproperties + end interface + + interface + pure subroutine icechanges_poreonly(nz,z,typeF,typeG,avdrhoP,ypp,B,porefill) + implicit none + integer, intent(IN) :: nz, typeF, typeG + real(8), intent(IN) :: z(nz), ypp(nz), avdrhoP, B + real(8), intent(INOUT) :: porefill(nz) + end subroutine icechanges_poreonly + end interface + + interface + pure subroutine icechanges(nz,z,typeF,avdrho,avdrhoP,ypp, & + & Diff,porosity,icefrac,bigstep,zdepthT,porefill,typeG) + implicit none + integer, intent(IN) :: nz, typeF, typeG + real(8), intent(IN) :: z(nz), ypp(nz), avdrho, avdrhoP + real(8), intent(IN) :: Diff, porosity, icefrac, bigstep + real(8), intent(INOUT) :: zdepthT, porefill(nz) + end subroutine icechanges + end interface + + interface + subroutine compactoutput(unit,porefill,nz) + implicit none + integer, intent(IN) :: unit,nz + real(8), intent(IN) :: porefill(nz) + end subroutine compactoutput + end interface + + !end of fast_subs_univ + !begin derivs.f90 + + interface + subroutine deriv1(z,nz,y,y0,yNp1,yp) + implicit none + integer, intent(IN) :: nz + real(8), intent(IN) :: z(nz),y(nz),y0,yNp1 + real(8), intent(OUT) :: yp(nz) + end subroutine deriv1 + end interface + + interface + subroutine deriv2_full(z,nz,a,b,a0,b0,bNp1,yp2) + implicit none + integer, intent(IN) :: nz + real(8), intent(IN) :: z(nz),a(nz),b(nz),a0,b0,bNp1 + real(8), intent(OUT) :: yp2(nz) + end subroutine deriv2_full + end interface + + interface + subroutine deriv2_simple(z,nz,y,y0,yNp1,yp2) + implicit none + integer, intent(IN) :: nz + real(8), intent(IN) :: z(nz),y(nz),y0,yNp1 + real(8), intent(OUT) :: yp2(nz) + end subroutine deriv2_simple + end interface + interface + real(8) pure function deriv1_onesided(j,z,nz,y) + implicit none + integer, intent(IN) :: nz,j + real(8), intent(IN) :: z(nz),y(nz) + end function deriv1_onesided + end interface + + !end of derivs.f90 + !begin fast_subs_exper.f90 + + interface + subroutine icelayer_exper(bigstep, nz, thIn, rhoc, z, porosity, pfrost, & + & zdepthT, icefrac, zdepthF, zdepthE, porefill, Tmean, Tampl, Diff, zdepthG) + use miscparameters, only : d2r, NMAX, icedensity + implicit none + integer, intent(IN) :: nz + real(8), intent(IN) :: bigstep + real(8), intent(IN) :: thIn, rhoc, z(NMAX), porosity, pfrost + real(8), intent(INOUT) :: zdepthT, zdepthF, porefill(nz) + real(8), intent(OUT) :: zdepthE + real(8), intent(IN) :: icefrac, Diff, Tmean, Tampl + real(8), intent(OUT) :: zdepthG + end subroutine icelayer_exper + end interface + + interface + subroutine ajsub_exper(typeT, nz, z, ti, rhocv, pfrost, Tmean, Tampl, & + & avdrho, avdrhoP, zdepthE, typeF, zdepthF, ypp, porefill, & + & B, typeG, zdepthG) + use miscparameters, only : NMAX, solsperyear, marsDay + implicit none + integer, intent(IN) :: nz, typeT + real(8), intent(IN) :: z(NMAX), ti(NMAX), rhocv(NMAX), pfrost + real(8), intent(IN) :: Tmean, Tampl + real(8), intent(OUT) :: avdrho, avdrhoP + real(8), intent(OUT) :: zdepthE + integer, intent(OUT) :: typeF + real(8), intent(INOUT) :: zdepthF + real(8), intent(OUT) :: ypp(nz) + real(8), intent(IN) :: porefill(nz) + real(8), intent(IN) :: B + integer, intent(OUT) :: typeG + real(8), intent(OUT) :: zdepthG + real(8), external :: Tsurface, psv + real(8), external :: equildepth + end subroutine ajsub_exper + end interface + + !end fast_subs_exper + + ! Other + interface + pure function flux_mars77(R,decl,HA,albedo,fracir,fracscat) + implicit none + real*8 flux_mars77 + real*8, intent(IN) :: R,decl,HA,albedo,fracIR,fracScat + end function flux_mars77 + end interface + + interface + pure function colint(y,z,nz,i1,i2) + implicit none + integer, intent(IN) :: nz, i1, i2 + real(8), intent(IN) :: y(nz),z(nz) + real(8) colint + end function colint + end interface + + interface + subroutine dyn_ss_ice_m(ssi_depth_in,T1,T_in,nsoil,thIn,p0,pfrost,porefill_in,porefill,ssi_depth) + implicit none + integer, intent(IN) :: nsoil + real(8), intent(IN) :: thIn,ssi_depth_in,T1 + real(8), intent(IN) :: p0(1), pfrost(1) + real(8), intent(IN) :: T_in(nsoil) + real(8), intent(INOUT) :: porefill(nsoil,1) + real(8), intent(INOUT) :: porefill_in(nsoil,1) + real(8), intent(INOUT) :: ssi_depth + end subroutine dyn_ss_ice_m + end interface + + interface + subroutine dyn_ss_ice_m_wrapper(ngrid,nsoil,tHIn,p0,pfrost,T_in,ssi_depth_in,porefill_in,porefill,ssi_depth) + implicit none + integer, intent(IN) :: nsoil,ngrid + real(8), intent(IN) :: thIn(ngrid),ssi_depth_in(ngrid) + real(8), intent(IN) :: p0(ngrid), pfrost(ngrid) + real(8), intent(IN) :: T_in(nsoil,ngrid) + real(8), intent(OUT) :: porefill(nsoil,ngrid) + real(8), intent(IN) :: porefill_in(nsoil,ngrid) + real(8), intent(OUT) :: ssi_depth(ngrid) + end subroutine dyn_ss_ice_m_wrapper + end interface + +end module allinterfaces Index: trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_mars.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_mars.F90 (revision 3493) +++ trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_mars.F90 (revision 3493) @@ -0,0 +1,309 @@ +module thermalmodelparam_mars + ! parameters for thermal model + ! they are only used in the subroutines below + real(8), parameter :: dt = 0.02 ! in units of Mars solar days + !real(8), parameter :: Fgeotherm = 0. + real(8), parameter :: Fgeotherm = 0 !0.028 ! [W/m^2] + real(8), parameter :: Lco2frost=6.0e5, co2albedo=0.60, co2emiss=1. + real(8), parameter :: emiss0 = 1. ! emissivity of dry surface + integer, parameter :: EQUILTIME =15 ! [Mars years] +end module thermalmodelparam_mars + + +!***************************************************** +! Subroutines for fast method +! written by Norbert Schorghofer 2007-2011 +!***************************************************** + + +subroutine icelayer_mars(bigstep,nz,NP,thIn,rhoc,z,porosity,pfrost, & + & Tb,zdepthF,zdepthE,porefill,Tmean1,Tmean3,zdepthG, & + & albedo,p0,icefrac,zdepthT,avrho1, & + & avrho1prescribed) +!************************************************************************* +! bigstep = time step [Earth years] +! latitude [degree] +!************************************************************************* + use miscparameters, only : d2r, NMAX, icedensity + use allinterfaces, except_this_one => icelayer_mars + !use omp_lib + implicit none + integer, intent(IN) :: nz, NP + real(8), intent(IN) :: bigstep + real(8), intent(IN) :: thIn(NP), rhoc(NP), z(NMAX), porosity, pfrost(NP) + real(8), intent(INOUT) :: porefill(nz,NP), zdepthF(NP), zdepthT(NP) + real(8), intent(INOUT) :: Tb(nz) + real(8), intent(OUT), dimension(NP) :: zdepthE, Tmean1, Tmean3, zdepthG + real(8), intent(IN), dimension(NP) :: albedo, p0 + real(8), intent(IN) :: icefrac + real(8), intent(OUT) :: avrho1(NP) + real(8), intent(IN), optional :: avrho1prescribed(NP) ! <0 means absent + + integer k, typeF, typeG, typeT, j, jump, typeP + real(8) fracIR, fracDust, ti(NMAX), rhocv(NMAX) + real(8) Diff, ypp(nz), avdrho(NP), avdrhoP(NP), B, deltaz + real(8), SAVE :: avdrho_old(100), zdepth_old(100) ! NP<=100 + logical mode2 + + !$omp parallel & + !$omp private(Diff,fracIR,fracDust,B,typeT,j,ti,rhocv,typeF,jump,typeG) + !$omp do + do k=1,NP ! big loop + + Diff = 4e-4*600./p0(k) + fracIR = 0.04*p0(k)/600.; fracDust = 0.02*p0(k)/600. + B = Diff*bigstep*86400.*365.24/(porosity*icedensity) + + typeT = -9 + if (zdepthT(k)>=0. .and. zdepthT(k)zdepthT(k)) then ! ice + typeT = j ! first point in ice + exit + endif + enddo + endif + + ! call assignthermalproperties(nz,thIn(k),rhoc(k),ti,rhocv, & + ! & typeT,icefrac,porosity,porefill(:,k)) + + !----run thermal model + call ajsub_mars(typeT, albedo(k), pfrost(k), nz, z, & + & ti, rhocv, fracIR, fracDust, p0(k), & + & avdrho(k), avdrhoP(k), avrho1(k), Tb(k), zdepthE(k), typeF, & + & zdepthF(k), ypp, porefill(:,k), Tmean1(k), Tmean3(k), B, & + & typeG, zdepthG(k), avrho1prescribed(k)) + + if (typeF*zdepthF(k)<0.) stop 'error in icelayer_mars' + ! diagnose + if (zdepthT(k)>=0.) then + jump = 0 + do j=1,nz + if (zdepth_old(k)z(j)) jump = jump+1 + enddo + else + jump = -9 + endif + if (zdepthT(k)>=0. .and. avdrho(k)*avdrho_old(k)<0.) then + write(34,*) '# zdepth arrested' + if (jump>1) write(34,*) '# previous step was too large',jump + endif +! write(34,'(f8.3,1x,f6.2,1x,f11.5,2(1x,g11.4),1x,i3)') & +! & bigstep,latitude(k),zdepthT(k),avdrho(k),avdrhoP(k),jump + zdepth_old(k) = zdepthT(k) + avdrho_old(k) = avdrho(k) + +!----mode 2 growth + typeP = -9; mode2 = .FALSE. + do j=1,nz + if (porefill(j,k)>0.) then + typeP = j ! first point with ice + exit + endif + enddo + if (typeT>0 .and. typeP>2 .and. zdepthE(k)>0.) then + if (porefill(typeP,k)>=1. .and. porefill(typeP-1,k)==0. .and. & + & zdepthE(k)2*(z(typeP)-z(typeP-1))) then ! trick that avoids oscillations + deltaz = -avdrhoP(k)/z(typeP)*18./8314.*B ! conservation of mass + if (deltaz>z(typeP)-z(typeP-1)) then ! also implies avdrhoP<0. + mode2 = .TRUE. + endif + endif + endif + + !call icechanges_poreonly(nz,z,typeF,typeG,avdrhoP(k),ypp,B,porefill(:,k)) + call icechanges(nz,z(:),typeF,avdrho(k),avdrhoP(k),ypp(:), & + & Diff,porosity,icefrac,bigstep,zdepthT(k),porefill(:,k),typeG) + + if (mode2 .and. porefill(typeP,k)>=1. .and. porefill(typeP-1,k)==0.) then ! nothing changed + porefill(typeP-1,k)=1. ! paste a layer + ! write(34,*) '# mode 2 growth occurred',typeP,typeF,typeT + endif + + enddo ! end of big loop + !$omp end do + !$omp end parallel +end subroutine icelayer_mars + + + +subroutine ajsub_mars(typeT, albedo0, pfrost, nz, z, ti, rhocv, & + & fracIR, fracDust, patm, avdrho, avdrhoP, avrho1, & + & Tb, zdepthE, typeF, zdepthF, ypp, porefill, Tmean1, Tmean3, & + & B, typeG, zdepthG, avrho1prescribed) +!*********************************************************************** +! A 1D thermal model that returns various averaged quantities +! +! Tb<0 initializes temperatures +! Tb>0 initializes temperatures with Tb +!*********************************************************************** + use miscparameters + use thermalmodelparam_mars + use allinterfaces, except_this_one => ajsub_mars + implicit none + integer, intent(IN) :: nz, typeT +! real(8), intent(IN) :: latitude ! in radians + real(8), intent(IN) :: albedo0, pfrost, z(NMAX) + real(8), intent(IN) :: ti(NMAX), rhocv(NMAX), fracIR, fracDust, patm + real(8), intent(IN) :: porefill(nz) + real(8), intent(OUT) :: avdrho, avdrhoP ! difference in annual mean vapor density + real(8), intent(OUT) :: avrho1 ! mean vapor density on surface + real(8), intent(INOUT) :: Tmean1 + real(8), intent(INOUT) :: Tb(nz) + integer, intent(OUT) :: typeF ! index of depth below which filling occurs + real(8), intent(OUT) :: zdepthE, zdepthF + real(8), intent(OUT) :: ypp(nz) ! (d rho/dt)/Diff + real(8), intent(OUT) :: Tmean3, zdepthG + real(8), intent(IN) :: B ! just passing through + integer, intent(OUT) :: typeG + real(8), intent(IN), optional :: avrho1prescribed + real(8), parameter :: a = 1.52366 ! Mars semimajor axis in A.U. + integer nsteps, n, i, nm, typeP + real(8) tmax, time, Qn, Qnp1, tdays + real(8) marsR, marsLs, marsDec, HA + real(8) Tsurf, Tco2frost, albedo, Fsurf, m, dE, emiss, T(NMAX) + real(8) Told(nz), Fsurfold, Tsurfold, Tmean0, avrho2 + real(8) rhosatav0, rhosatav(nz), rlow + real(8), external :: psv, tfrostco2 + + tmax = EQUILTIME*solsperyear + nsteps = int(tmax/dt) ! calculate total number of timesteps + + Tco2frost = tfrostco2(patm) + + !if (Tb<=0.) then ! initialize + !Tmean0 = 210.15 ! black-body temperature of planet + !Tmean0 = (589.*(1.-albedo0)*cos(latitude)/(pi*emiss0*sigSB))**0.25 ! better estimate + !Tmean0 = Tmean0-5. + !write(34,*) '# initialized with temperature estimate at',latitude/d2r,'of',Tmean0,'K' + !T(1:nz) = Tmean0 + !else + T(1:nz) = Tb + ! not so good when Fgeotherm is on + !endif + + albedo = albedo0 + emiss = emiss0 + do i=1,nz + if (T(i)Tco2frost) then + ! call conductionQ(nz,z,dt*marsDay,Qn,Qnp1,T,ti,rhocv,emiss, & + ! & Tsurf,Fgeotherm,Fsurf) + endif + if (Tsurf0.) then ! CO2 condensation + T(1:nz) = Told(1:nz) + !call conductionT(nz,z,dt*marsDay,T,Tsurfold,Tco2frost,ti, & + !& rhocv,Fgeotherm,Fsurf) + Tsurf = Tco2frost + ! dE = (- Qn - Qnp1 + Fsurfold + Fsurf + & + ! & emiss*sigSB*(Tsurfold**4+Tsurf**4)/2. + m = m + dt*marsDay*dE/Lco2frost + endif + if (Tsurf>Tco2frost .or. m<=0.) then + albedo = albedo0 + emiss = emiss0 + else + albedo = co2albedo + emiss = co2emiss + endif + !Qn=Qnp1 + + if (time>=tmax-solsperyear) then + Tmean1 = Tmean1 + Tsurf + Tmean3 = Tmean3 + T(nz) + avrho1 = avrho1 + min(psv(Tsurf),pfrost)/Tsurf + rhosatav0 = rhosatav0 + psv(Tsurf)/Tsurf + do i=1,nz + rhosatav(i) = rhosatav(i) + psv(T(i))/T(i) + enddo + nm = nm+1 + endif + + enddo ! end of time loop + + avrho1 = avrho1/nm + if (present(avrho1prescribed)) then + if (avrho1prescribed>=0.) avrho1=avrho1prescribed + endif + Tmean1 = Tmean1/nm; Tmean3 = Tmean3/nm + rhosatav0 = rhosatav0/nm; rhosatav(:)=rhosatav(:)/nm + if (typeT>0) then + avrho2 = rhosatav(typeT) + else + avrho2 = rhosatav(nz) ! no ice + endif + avdrho = avrho2-avrho1 + typeP = -9 + do i=1,nz + if (porefill(i)>0.) then + typeP = i ! first point with ice + exit + endif + enddo + if (typeP>0) then + avdrhoP = rhosatav(typeP) - avrho1 + else + avdrhoP = -9999. + end if + + zdepthE = equildepth(nz, z, rhosatav, rhosatav0, avrho1) + + if (Fgeotherm>0.) then + Tb = Tmean1 + typeG = 1 ! will be overwritten by depths_avmeth + rlow = 2*rhosatav(nz)-rhosatav(nz-1) + else + Tb = T(nz) + typeG = -9 + rlow = rhosatav(nz-1) + endif + call depths_avmeth(typeT, nz, z, rhosatav(:), rhosatav0, rlow, avrho1, & + & porefill(:), typeF, zdepthF, B, ypp(:), typeG, zdepthG) + +end subroutine ajsub_mars + + + +subroutine outputmoduleparameters + use miscparameters + use thermalmodelparam_mars + implicit none +! print *,'Parameters stored in modules' +! print *,' Ice bulk density',icedensity,'kg/m^3' +! print *,' dt=',dt,'Mars solar days' +! print *,' Fgeotherm=',Fgeotherm,'W/m^2' +! write(6,'(2x,a27,1x,f5.3)') 'Emissivity of dry surface=',emiss0 +! write(6,'(2x,a12,1x,f5.3,2x,a16,1x,f5.3)') 'CO2 albedo=',co2albedo,'CO2 emissivity=',co2emiss +! print *,' Thermal model equilibration time',EQUILTIME,'Mars years' +end subroutine outputmoduleparameters + + + Index: trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_univ.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_univ.F90 (revision 3493) +++ trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_univ.F90 (revision 3493) @@ -0,0 +1,383 @@ +!***************************************************** +! Commonly used subroutines for fast method +! written by Norbert Schorghofer 2007-2011 +!***************************************************** + +pure function zint(y1,y2,z1,z2) + ! interpolate between two values, y1*y2<0 + implicit none + real(8), intent(IN) :: y1,y2,z1,z2 + real(8) zint + zint = (y1*z2 - y2*z1)/(y1-y2) +end function zint + + + +pure function equildepth(nz, z, rhosatav, rhosatav0, avrho1) +!*********************************************************************** +! returns equilibrium depth for given ice content +! this is not the true (final) equilibrium depth +!*********************************************************************** + use allinterfaces, except_this_one => equildepth + implicit none + integer, intent(IN) :: nz + real(8), intent(IN) :: z(nz), rhosatav(nz) + real(8), intent(IN) :: rhosatav0, avrho1 + integer i, typeE + real(8) equildepth ! =zdepthE + !real(8), external :: zint ! defined in allinterfaces.mod + + typeE = -9; equildepth = -9999. + do i=1,nz + if (rhosatav(i) <= avrho1) then + typeE=i + exit + endif + enddo + if (typeE>1) then ! interpolate + equildepth=zint(avrho1-rhosatav(typeE-1),avrho1-rhosatav(typeE),z(typeE-1),z(typeE)) + end if + if (typeE==1) equildepth=zint(avrho1-rhosatav0,avrho1-rhosatav(1),0.d0,z(1)) + if (equildepth>z(nz)) equildepth=-9999. ! happens very rarely +end function equildepth + + + +subroutine depths_avmeth(typeT, nz, z, rhosatav, rhosatav0, rlow, avrho1, & + & porefill, typeF, zdepthF, B, ypp, typeG, zdepthG) +!*********************************************************************** +! returns interface depth and ypp +! also returns lower geothermally limited boundary, if applicable +! +! this is a crucial subroutine +! +! B = Diff*bigstep/(porosity*icedensity) [SI units] +!*********************************************************************** + use allinterfaces, except_this_one => depths_avmeth + implicit none + integer, intent(IN) :: nz, typeT + real(8), intent(IN), dimension(nz) :: z, rhosatav, porefill + real(8), intent(IN) :: rhosatav0, rlow, avrho1 + integer, intent(OUT) :: typeF ! index of depth below which filling occurs + real(8), intent(INOUT) :: zdepthF + real(8), intent(IN) :: B + real(8), intent(OUT) :: ypp(nz), zdepthG + integer, intent(INOUT) :: typeG ! positive on input when Fgeotherm>0 + + integer i, typeP, nlb, newtypeG + real(8) eta(nz), Jpump1, help(nz), yp(nz), zdepthFold, ap_one, ap(nz) + real(8) leak, cumfill, cumfillabove + + if (typeT<0) then + nlb = nz + do i=1,nz + eta(i) = constriction(porefill(i)) + enddo + else + !nlb = typeT-1 + nlb = typeT ! changed 2010-09-29 + do i=1,typeT-1 + eta(i) = constriction(porefill(i)) + enddo + do i=typeT,nz + eta(i)=0. + enddo + end if + +!-fill depth + zdepthFold = zdepthF + typeF = -9; zdepthF = -9999. + call deriv1(z,nz,rhosatav,rhosatav0,rlow,yp) ! yp also used below + do i=1,nlb + Jpump1 = (rhosatav(i)-avrho1)/z(i) ! <0 when stable + ! yp is always <0 + help(i) = Jpump1 - eta(i)*yp(i) + leak = porefill(i)/B*(z(i)-zdepthFold)/(18./8314.) + !help(i) = help(i)-leak ! optional + if (help(i) <= 0.) then + typeF=i + !print *,'#',typeF,Jpump1,eta(typeF)*yp(typeF),leak + exit + endif + enddo + if (typeF>1) zdepthF = zint(help(typeF-1),help(typeF),z(typeF-1),z(typeF)) + if (typeF==1) zdepthF=z(1) + + +!-depth to shallowest perennial ice + typeP = -9 + do i=1,nz + if (porefill(i)>0.) then + typeP = i ! first point with ice + exit + endif + enddo + +!-calculate ypp + !call deriv1(z,nz,rhosatav,rhosatav0,rlow,yp) + call deriv1(z,nz,eta(:),1.d0,eta(nz-1),ap) + if (typeP>0 .and. typeP0) typeG=-9 + if (typeG<0) zdepthG=-9999. + if (typeG>0 .and. typeT<0) then + typeG=-9 + do i=2,nz + if (yp(i)>0.) then ! first point with reversed flux + typeG=i + zdepthG=zint(yp(i-1),yp(i),z(i-1),z(i)) + !zdepthG=z(typeG) + exit + endif + enddo + else + typeG = -9 + endif + if (typeG>0 .and. typeT<0) then + cumfillabove = colint(porefill(:)/eta(:),z,nz,typeG-1,nz) + newtypeG = -9 + do i=typeG,nz + if (minval(eta(i:nz))<=0.) cycle ! eta=0 means completely full + cumfill=colint(porefill(:)/eta(:),z,nz,i,nz) + if (cumfilltypeG) then + write(34,*) '# adjustment to geotherm depth by',i-typeG + zdepthG = zint(yp(i-1)*18./8314.*B-cumfillabove, & ! no good + & yp(i)*18./8314.*B-cumfill,z(i-1),z(i)) + if (zdepthG>z(i) .or. zdepthG0) typeG=newtypeG + end if + ! if typeG>0, then all ice at and below typeG should be erased +end subroutine depths_avmeth + + + +pure function constriction(porefill) +! specify constriction function here, 0<=eta<=1 + implicit none + real(8), intent(IN) :: porefill + real(8) eta, constriction + if (porefill<=0.) eta = 1. + if (porefill>0. .and. porefill<1.) then + ! eta = 1. + ! eta = 1-porefill + eta = (1-porefill)**2 ! Hudson et al., JGR, 2009 + endif + if (porefill>=1.) eta = 0. + constriction = eta +end function constriction + + + +pure subroutine icechanges_poreonly(nz,z,typeF,typeG,avdrhoP,ypp,B,porefill) + use allinterfaces, except_this_one => icechanges_poreonly + implicit none + integer, intent(IN) :: nz, typeF, typeG + real(8), intent(IN) :: z(nz), ypp(nz), avdrhoP, B + real(8), intent(INOUT) :: porefill(nz) + integer j, erase, newtypeP, ub + real(8) integ + + !----retreat + ! avdrhoP>0 is outward loss from zdepthP + ! avdrhoP<0 means gain at zdepthP or no ice anywhere + if (avdrhoP>0.) then + erase=0 + do j=1,nz + if (typeF>0 .and. j>=typeF) exit ! don't retreat beyond typeF + integ = colint(porefill(1:nz)*z(1:nz),z(1:nz),nz,1,j) + erase = j + if (integ>B*avdrhoP*18./8314.) exit + end do + if (erase>0) porefill(1:erase)=0. + endif + + ! new depth + newtypeP = -9 + do j=1,nz + if (porefill(j)>0.) then + newtypeP = j ! first point with ice + exit + endif + enddo + + !----diffusive filling + ub = typeF + if (newtypeP>0 .and. typeF>0 .and. newtypeP0) then + do j=ub,nz + ! B=Diff/(porosity*icedensity)*86400*365.24*bigstep + porefill(j) = porefill(j) + B*ypp(j) + if (porefill(j)<0.) porefill(j)=0. + if (porefill(j)>1.) porefill(j)=1. + enddo + end if + + !----enact bottom boundary + if (typeG>0) porefill(typeG:nz)=0. + +end subroutine icechanges_poreonly + + + +pure subroutine icechanges(nz,z,typeF,avdrho,avdrhoP,ypp, & + & Diff,porosity,icefrac,bigstep,zdepthT,porefill,typeG) +!*********************************************************** +! advances ice table, advances interface, grows pore ice +! +! a crucial subroutine +!*********************************************************** + use miscparameters, only : icedensity + use allinterfaces, except_this_one => icechanges + implicit none + integer, intent(IN) :: nz, typeF, typeG + real(8), intent(IN) :: z(nz), ypp(nz), avdrho, avdrhoP + real(8), intent(IN) :: Diff, porosity, icefrac, bigstep + real(8), intent(INOUT) :: zdepthT, porefill(nz) + integer j, erase, newtypeP, ub, typeP, typeT + real(8) B, beta, integ + + B = Diff*bigstep*86400.*365.24/(porosity*icedensity) + + ! advance ice table, avdrho>0 is retreat + if (zdepthT>=0. .and. avdrho>0.) then + typeP=-9999; typeT=-9999 + do j=1,nz + if (z(j)>zdepthT) then + typeT=j + exit + endif + enddo + do j=1,nz + if (porefill(j)>0.) then + typeP=j + exit + endif + enddo + if (typeP==typeT) then ! new 2011-09-01 + beta = (1-icefrac)/(1-porosity)/icefrac + beta = Diff*bigstep*beta*86400*365.24/icedensity + zdepthT = sqrt(2*beta*avdrho*18./8314. + zdepthT**2) + endif + endif + if (zdepthT>z(nz)) zdepthT=-9999. + + ! advance interface, avdrhoP>0 is loss from zdepthP + if (avdrhoP>0.) then + erase=0 + do j=1,nz + if (typeF>0 .and. j>=typeF) exit ! don't retreat beyond typeF + if (zdepthT>=0. .and. z(j)>zdepthT) exit + integ = colint(porefill(1:nz)*z(1:nz),z(1:nz),nz,1,j) + erase = j + if (integ>B*avdrhoP*18./8314.) exit + end do + if (erase>0) porefill(1:erase)=0. + endif + + ! new depth + newtypeP = -9 + do j=1,nz + if (zdepthT>=0. .and. z(j)>zdepthT) exit + if (porefill(j)>0.) then + newtypeP = j ! first point with pore ice + exit + endif + enddo + + ! diffusive filling + ub = typeF + if (newtypeP>0 .and. typeF>0 .and. newtypeP0) then + do j=ub,nz + porefill(j) = porefill(j) + B*ypp(j) + if (porefill(j)<0.) porefill(j)=0. + if (porefill(j)>1.) porefill(j)=1. + if (zdepthT>=0. .and. z(j)>zdepthT) exit + enddo + end if + + ! below icetable + if (zdepthT>=0.) then + do j=1,nz + if (z(j)>zdepthT) porefill(j) = icefrac/porosity + enddo + else + ! geothermal lower boundary + if (typeG>0) porefill(typeG:nz)=0. + end if +end subroutine icechanges + + +subroutine assignthermalproperties(nz,thIn,rhoc, & + & ti,rhocv,typeT,icefrac,porosity,porefill) +!********************************************************* +! assign thermal properties of soil +!********************************************************* + implicit none + integer, intent(IN) :: nz + integer, intent(IN), optional :: typeT + real(8), intent(IN), optional :: icefrac + real(8), intent(IN) :: thIn, rhoc + real(8), intent(IN), optional :: porosity, porefill(nz) + real(8), intent(OUT) :: ti(nz), rhocv(nz) + integer j + real(8) newrhoc, newti, fill + real(8), parameter :: NULL=0. + + ti(1:nz) = thIn + rhocv(1:nz) = rhoc + if (typeT>0) then + call soilthprop(porosity,NULL,rhoc,thIn,2,newrhoc,newti,icefrac) + rhocv(typeT:nz) = newrhoc + ti(typeT:nz) = newti + endif + do j=1,nz + fill = porefill(j) ! off by half point + if (fill>0. .and. (typeT<0 .or. (typeT>0 .and. j0.d0) then + solarAttenuation = (1.- fracIR - fracScat)**(1./max(sinbeta,0.04)) + Qdir = Qo*sinbeta*solarAttenuation + Qscat = 0.5*Qo*fracScat + Q = (1.-albedo)*(Qdir + Qscat) + else + Q = 0. + endif + + ! atmospheric IR contribution based on Kieffer et al. (1977), JGR 82, 4249 + Fout = 1.*sigSB*150**4 ! matters only in polar region + Qlw = fracIR*max(Qo*sinbetaNoon,Fout) + + ! net flux = short-wavelength + long-wavelength irradiance + flux_mars77 = Q + Qlw +end function flux_mars77 + + + +pure subroutine flux_mars2(R,decl,latitude,HA,fracIR,fracScat, & + & SlopeAngle,azFac,emax,Qdir,Qscat,Qlw) +!***************************************************************************** +! flux_mars2: Insolation (incoming solar radiation) at Mars on sloped surface; +! returns several irradiances +! +! INPUTS: +! R: distance from sun [AU] +! decl: planetocentric solar declination [radians] +! latitude: [radians] +! HA: hour angle (radians from noon, clockwise) +! fracIR: fraction of absorption +! fracScat: fraction of scattering +! SlopeAngle: >0, [radians] +! azFac: azimuth of topographic gradient (radians east of north) +! azFac=0 is south-facing +! emax: maximum horizon elevation in direction of azimuth [radians] +! OUTPUTS: +! Qdir: direct incoming short-wavelength irradiance [W/m^2] +! Qscat: diffuse short-wavelength irradiance from atmosphere [W/m^2] +! Qlw: diffuse long-wavelength irradiance from atmosphere [W/m^2] +!***************************************************************************** + implicit none + real(8), parameter :: pi=3.1415926535897932, So=1365. + real(8), parameter :: sigSB=5.6704d-8 + real(8), intent(IN) :: R,decl,latitude,HA,SlopeAngle,azFac,emax + real(8), intent(IN) :: fracIR,fracScat + real(8), intent(OUT) :: Qdir,Qscat,Qlw + real(8) c1,s1,solarAttenuation,Qo + real(8) sinbeta,sinbetaNoon,cosbeta,sintheta,azSun,buf,Fout + + c1 = cos(latitude)*cos(decl) + s1 = sin(latitude)*sin(decl) + ! beta = 90 minus incidence angle for horizontal surface + ! beta = elevation of sun above (horizontal) horizon + sinbeta = c1*cos(HA) + s1 + sinbetaNoon = c1 + s1 + sinbetaNoon = max(sinbetaNoon,0.d0) ! horizon + + cosbeta = sqrt(1-sinbeta**2) + buf = (sin(decl)-sin(latitude)*sinbeta)/(cos(latitude)*cosbeta) + if (buf>+1.) buf=+1.d0 ! roundoff + if (buf<-1.) buf=-1.d0 ! roundoff + azSun = acos(buf) + if (sin(HA)>=0) azSun = 2*pi-azSun + +! theta = 90 minus incidence angle for sloped surface + sintheta = cos(SlopeAngle)*sinbeta - & + & sin(SlopeAngle)*cosbeta*cos(azSun-azFac) + if (cosbeta==0) sintheta=cos(SlopeAngle)*sinbeta ! zenith, where azSun=NaN + + if (sintheta<0.) sintheta=0. ! local horizon (self-shadowing) + if (sinbeta<0.) sintheta=0. ! horizontal horizon at infinity + if (sinbeta0.d0) then + Qscat = 0.5*fracScat*Qo + else + Qscat = 0. + endif + +! For a horizontal and unobstructed surface +! absorbed flux = (1-albedo)*(Qdir+Qscat) + emiss*Qlw +! +! For a tilted surface +! absorbed flux = (1-albedo)*(Qdir+Qscat*viewfactor) + emiss*Qlw*viewfactor +! then add irradiance from land in field of view +! in the case of a planar slope, viewfactor = cos(SlopeAngle/2.)**2 +end subroutine flux_mars2 Index: trunk/LMDZ.COMMON/libf/evolution/NS_generalorbit.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/NS_generalorbit.F90 (revision 3493) +++ trunk/LMDZ.COMMON/libf/evolution/NS_generalorbit.F90 (revision 3493) @@ -0,0 +1,54 @@ +!============================================================= +! Subroutine to return distance, longitude, and declination of +! the sun in planetocentric coordinates from orbital elements +! +! INPUTS: +! edays = time since perihelion (earth days) +! a = semimajor axis (AU) +! ecc = eccentricity +! omega = Ls of perihelion, relative to equinox (radians) +! eps = obliquity (radians) +! +! OUTPUTS: +! Ls = areocentric longitude (radians) +! dec = planetocentric solar declination (radians) +! r = heliocentric distance (AU) +!============================================================= + + SUBROUTINE generalorbit(edays,a,ecc,omega,eps,Ls,dec,r) + implicit none + real*8 edays,a,ecc,omega,eps ! input + real*8 Ls,dec,r ! output + real*8 pi,d2r + parameter (pi=3.1415926535897932,d2r=pi/180.d0) + integer j + real*8 M,E,nu,T,Eold + +! T = orbital period (days) + T = sqrt(4*pi**2/(6.674e-11*1.989e30)*(a*149.598e9)**3)/86400. + +! M = mean anomaly (radians) + M = 2.*pi*edays/T ! M=0 at perihelion + +! E = eccentric anomaly +! solve M = E - ecc*sin(E) by newton method + E = M + do j=1,10 + Eold = E + E = E - (E - ecc*sin(E) - M)/(1.-ecc*cos(E)) + if (abs(E-Eold)<1.e-8) exit + enddo + +! nu = true anomaly + !nu = acos(cos(E)-ecc/(1.-ecc*cos(E))) + !nu = sqrt(1-ecc^2)*sin(E)/(1.-ecc*cos(E)) + !nu = atan(sqrt(1-ecc^2)*sin(E)/(1-cos(E))) + nu = 2.*atan(sqrt((1.+ecc)/(1.-ecc))*tan(E/2.)) + + !r = a*(1.-ecc**2)/(1.+ecc*cos(nu)) + r = a*(1-ecc*cos(E)) + Ls = mod(nu + omega,2.*pi) + dec = asin(sin(eps)*sin(Ls)) + + END + Index: trunk/LMDZ.COMMON/libf/evolution/NS_grids.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/NS_grids.F90 (revision 3493) +++ trunk/LMDZ.COMMON/libf/evolution/NS_grids.F90 (revision 3493) @@ -0,0 +1,96 @@ +subroutine setgrid(nz,z,zmax,zfac) + ! construct regularly or geometrically spaced 1D grid + ! z(n)-z(1) = 2*z(1)*(zfac**(n-1)-1)/(zfac-1) + ! choice of z(1) and z(2) is compatible with conductionQ + implicit none + integer, intent(IN) :: nz + real(8), intent(IN) :: zfac, zmax + real(8), intent(OUT) :: z(nz) + integer i + real(8) dz + + dz = zmax/nz + do i=1,nz + z(i) = (i-0.5)*dz + enddo + + if (zfac>1.) then + dz = zmax/(3.+2.*zfac*(zfac**(nz-2)-1.)/(zfac-1.)) + z(1) = dz + z(2) = 3*z(1) + do i=3,nz + z(i) = (1.+zfac)*z(i-1) - zfac*z(i-2) + ! z(i) = z(i-1) + zfac*(z(i-1)-z(i-2)) ! equivalent + enddo + endif + +end subroutine setgrid + + + +real(8) function smartzfac(nz_max,zmax,nz_delta,delta) + ! output can be used as input to setgrid + ! produces zfac with desired number of points within depth delta + implicit none + integer, intent(IN) :: nz_max, nz_delta + real(8), intent(IN) :: zmax, delta + integer j + real(8) f, g, gprime + + if (nz_max2.) then + print *,'zfac=',smartzfac + stop 'unwanted result in smartzfac' + endif +end function smartzfac + + + +!-grid-dependent utility functions + +function colint(y,z,nz,i1,i2) + ! column integrates y + implicit none + integer, intent(IN) :: nz, i1, i2 + real(8), intent(IN) :: y(nz), z(nz) + real(8) colint + integer i + real(8) dz(nz) + + dz(1) = (z(2)-0.)/2 + do i=2,nz-1 + dz(i) = (z(i+1)-z(i-1))/2. + enddo + dz(nz) = z(nz)-z(nz-1) + colint = sum(y(i1:i2)*dz(i1:i2)) +end function colint + + + +subroutine heatflux_from_temperature(nz,z,T,k,H) + ! calculates heat flux from temperature profile + ! like k, the heat flux H is defined mid-point + implicit none + integer, intent(IN) :: nz + real(8), intent(IN) :: z(nz), T(nz), k(nz) + real(8), intent(OUT) :: H(nz) + integer j + + ! H(1) = -k(1)*(T(1)-Tsurf)/z(1) + H(1) = 0. ! to avoid ill-defined value + do j=2,nz + H(j) = -k(j)*(T(j)-T(j-1))/(z(j)-z(j-1)) + enddo +end subroutine heatflux_from_temperature Index: trunk/LMDZ.COMMON/libf/evolution/NS_misc_params.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/NS_misc_params.F90 (revision 3493) +++ trunk/LMDZ.COMMON/libf/evolution/NS_misc_params.F90 (revision 3493) @@ -0,0 +1,14 @@ +module miscparameters + ! miscellaneous parameters that are very constant + real(8), parameter :: pi=3.1415926535897932, d2r=pi/180. + integer, parameter :: NMAX=1000 + real(8), parameter :: marsDay=88775.244, solsperyear=668.60 + real(8), parameter :: icedensity=927. + real(8), parameter :: earthDay=86400. + real(8), parameter :: sigSB=5.6704e-8 + ! for reference here are some parameters that are hard coded + ! mass of H2O molecule = 18 + ! universal gas constant = 8314 + ! length of Earth year in days = 365.24 +end module miscparameters + Index: trunk/LMDZ.COMMON/libf/evolution/NS_psv.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/NS_psv.F90 (revision 3493) +++ trunk/LMDZ.COMMON/libf/evolution/NS_psv.F90 (revision 3493) @@ -0,0 +1,54 @@ + real*8 function psv(T) +! saturation vapor pressure of H2O ice [Pascal] +! input is temperature [Kelvin] + implicit none + real*8 T + +!-----parametrization 1 +! real*8 DHmelt,DHvap,DHsub,R,pt,Tt,C +! parameter (DHmelt=6008.,DHvap=45050.) +! parameter (DHsub=DHmelt+DHvap) ! sublimation enthalpy [J/mol] +! parameter (R=8.314,pt=6.11e2,Tt=273.16) +! C = (DHsub/R)*(1./T - 1./Tt) +! psv = pt*exp(-C) + +!-----parametrization 2 +! eq. (2) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005) +! differs from parametrization 1 by only 0.1% + real*8 A,B + parameter (A=-6143.7, B=28.9074) + psv = exp(A/T+B) ! Clapeyron + +!-----parametrization 3 +! eq. (7) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005) +! psv = exp(9.550426 - 5723.265/T + 3.53068*log(T) - 0.00728332*T) + + end + + + + real*8 function frostpoint(p) +! inverse of psv +! input is partial pressure [Pascal] +! output is temperature [Kelvin] + implicit none + real*8 p + +!-----inverse of parametrization 1 +! real*8 DHmelt,DHvap,DHsub,R,pt,Tt +! parameter (DHmelt=6008.,DHvap=45050.) +! parameter (DHsub=DHmelt+DHvap) +! parameter (R=8.314,pt=6.11e2,Tt=273.16) +! frostpoint = 1./(1./Tt-R/DHsub*log(p/pt)) + +!-----inverse of parametrization 2 +! inverse of eq. (2) in Murphy & Koop (2005) + real*8 A,B + parameter (A=-6143.7, B=28.9074) + frostpoint = A / (log(p) - B) + +!-----approximate inverse of parametrization 3 +! eq. (8) in Murphy & Koop (2005) +! frostpoint = (1.814625*log(p) + 6190.134)/(29.120 - log(p)) + + end Index: trunk/LMDZ.COMMON/libf/evolution/NS_psvco2.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/NS_psvco2.F90 (revision 3493) +++ trunk/LMDZ.COMMON/libf/evolution/NS_psvco2.F90 (revision 3493) @@ -0,0 +1,50 @@ + real*8 function psvco2(T) +! solid-vapor transition for CO2 +! returns saturation pressure in Pascal, input is temperature in Kelvin + implicit none + real*8 T + + psvco2 = exp(23.3494 - 3182.48/T)*100. + end + + + real*8 function tfrostco2(p) +! the inverse of function psvco2 +! input is pressure in Pascal, output is temperature in Kelvin + implicit none + real*8 p + tfrostco2 = 3182.48/(23.3494+log(100./p)) + end + + + + +!------------------------------------------------------------------------------ +! Antoine equation parameters from NIST, 154K-196K +! based on Giauque and Egan (1937) +! A=6.81228, B=1301.679, C=-3.494 +! p = 10**(A-(B/(T+C)))*1.e5 + +! Expressions from Int. Crit. Tabl. 3.207, based on many references +! mm2Pa = 133.32 +! -135C to -56.7C (138K-216K) +! p = 10**(-0.05223*26179.3/T + 9.9082)*mm2Pa +! -183C to -135C (90K-138K) +! p = 10**(-1275.62/T + 0.006833*T + 8.3071)*mm2Pa +! Expressions from Int. Crit. Tabl. 3.208, based on Henning +! -110 to -80C (163K-193K) +! p = 10**(- 1279.11/T + 1.75*log10(T) - 0.0020757*T + 5.85242)*mm2Pa +! p = 10**(- 1352.6/T + 9.8318)*mm2Pa + +! Mars book (1992), p959, fit by chapter authors +! p = exp(23.3494 - 3182.48/T)*100. ! 120-160K +! p = exp(25.2194 - 3311.57/T - 6.71e-3*T)*100 ! 100-170K +! Mars book (1992), p960, based on Miller & Smythe (1970) +! p = exp(26.1228 - 3385.26/T - 9.4461e-3*T)*100 ! 123-173K + +! Fray & Schmitt, PSS 57, 2053 (2009) +! A0=1.476e1, A1=-2.571e3, A2=-7.781e4, A3=4.325e6, A4=-1.207e8, A5=1.350e9 +! p = exp(A0+A1/T+A2/T**2+A3/T**3+A4/T**4+A5/T**5)*1e5 ! 40K-194.7K +! A0=1.861e1, A1=-4.154e3, A2=1.041e5 +! p = exp(A0 + A1/T + A2/T**2)*1e5 ! 194.7K-216.58K +!------------------------------------------------------------------------------ Index: trunk/LMDZ.COMMON/libf/evolution/NS_soilthprop_mars.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/NS_soilthprop_mars.F90 (revision 3493) +++ trunk/LMDZ.COMMON/libf/evolution/NS_soilthprop_mars.F90 (revision 3493) @@ -0,0 +1,152 @@ +subroutine soilthprop(porosity,fill,rhocobs,tiobs,layertype, & + & newrhoc,newti,icefrac) +!*********************************************************************** +! soilthprop: assign thermal properties of icy soil or dirty ice +! +! porositiy = void space / total volume +! rhof = density of free ice in space not occupied by regolith [kg/m^3] +! fill = rhof/icedensity <=1 (only relevant for layertype 1) +! rhocobs = heat capacity per volume of dry regolith [J/m^3] +! tiobs = thermal inertia of dry regolith [SI-units] +! layertype: 1=interstitial ice, 2=pure ice or ice with dirt +! 3=pure ice, 4=ice-cemented soil, 5=custom values +! icefrac: fraction of ice in icelayer (only relevant for layertype 2) +! output are newti and newrhoc +!*********************************************************************** + implicit none + integer, intent(IN) :: layertype + real(8), intent(IN) :: porosity, fill, rhocobs, tiobs + real(8), intent(OUT) :: newti, newrhoc + real(8), intent(IN) :: icefrac + real(8) kobs, cice, icedensity, kice + !parameter (cice=2000.d0, icedensity=926.d0, kice=2.4d0) ! unaffected by scaling + parameter (cice=1540.d0, icedensity=927.d0, kice=3.2d0) ! at 198 Kelvin + real(8) fA, ki0, ki, k + real(8), parameter :: kw=3. ! Mellon et al., JGR 102, 19357 (1997) + + kobs = tiobs**2/rhocobs + ! k, rhoc, and ti are defined in between grid points + ! rhof and T are defined on grid points + + newrhoc = -9999. + newti = -9999. + + select case (layertype) + case (1) ! interstitial ice + newrhoc = rhocobs + porosity*fill*icedensity*cice + if (fill>0.) then + !--linear addition (option A) + k = porosity*fill*kice + kobs + !--Mellon et al. 1997 (option B) + ki0 = porosity/(1/kobs-(1-porosity)/kw) + fA = sqrt(fill) + ki = (1-fA)*ki0 + fA*kice + !k = kw*ki/((1-porosity)*ki+porosity*kw) + else + k = kobs + endif + newti = sqrt(newrhoc*k) + + case (2) ! massive ice (pure or dirty ice) + newrhoc = rhocobs*(1.-icefrac)/(1.-porosity) + icefrac*icedensity*cice + k = icefrac*kice + (1.-icefrac)*kw + newti = sqrt(newrhoc*k) + + case (3) ! all ice, special case of layertype 2, which doesn't use porosity + newrhoc = icedensity*cice + k = kice + newti = sqrt(newrhoc*k) + + case (4) ! pores completely filled with ice, special case of layertype 1 + newrhoc = rhocobs + porosity*icedensity*cice + k = porosity*kice + kobs ! option A, end-member case of type 1, option A + !k = kw*kice/((1-porosity)*kice+porosity*kw) ! option B, harmonic average + newti = sqrt(newrhoc*k) + + case (5) ! custom values + ! values from Mellon et al. (2004) for ice-cemented soil + newrhoc = 2018.*1040. + k = 2.5 + newti = sqrt(newrhoc*k) + + case default + error stop 'invalid layer type' + + end select + +end subroutine soilthprop + + + +subroutine smartgrid(nz,z,zdepth,thIn,rhoc,porosity,ti,rhocv,layertype,icefrac) +!*********************************************************************** +! smartgrid: returns intelligently spaced grid and appropriate +! values of thermal inertia ti and rho*c in icy layer +! +! INPUTS: +! nz = number of grid points +! z = grid spacing for dry regolith +! (will be partially overwritten) +! zdepth = depth where ice table starts +! negative values indicate no ice +! rhoc = heat capacity per volume of dry regolith [J/m^3] +! thIn = thermal inertia of dry regolith [SI-units] +! porosity = void space / total volume +! layertypes are explained below +! icefrac = fraction of ice in icelayer +! +! OUTPUTS: z = grid coordinates +! vectors ti and rhocv +!*********************************************************************** + implicit none + integer, intent(IN) :: nz, layertype + real(8), intent(IN) :: zdepth, thIn, rhoc, porosity, icefrac + real(8), intent(INOUT) :: z(nz) + real(8), intent(OUT) :: ti(nz), rhocv(nz) + integer j, b + real(8) stretch, newrhoc, newti + real(8), parameter :: NULL=0. + + if (zdepth>0 .and. zdepth',n +! stop +! endif + bet=b(1) + u(1)=r(1)/bet + do j=2,n + gam(j)=c(j-1)/bet + bet=b(j)-a(j)*gam(j) +! if(bet.eq.0.)pause 'tridag failed' + u(j)=(r(j)-a(j)*u(j-1))/bet + enddo + do j=n-1,1,-1 + u(j)=u(j)-gam(j+1)*u(j+1) + enddo +! return + END +! (C) Copr. 1986-92 Numerical Recipes Software 0(9p#31(+. Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3492) +++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3493) @@ -460,2 +460,8 @@ == 05/11/2024 == JBC Correction in the formula to update the potential temperature introduced in r3442. + +== 05/11/2024 == JBC +- 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. Index: trunk/LMDZ.COMMON/libf/evolution/conductionQ.f90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/conductionQ.f90 (revision 3492) +++ (revision ) @@ -1,113 +1,0 @@ -subroutine conductionQ(nz,z,dt,Qn,Qnp1,T,ti,rhoc,emiss,Tsurf, & - & Fgeotherm,Fsurf) -!*********************************************************************** -! conductionQ: program to calculate the diffusion of temperature -! into the ground and thermal emission at the surface -! with variable thermal properties on irregular grid -! Crank-Nicolson scheme, flux conservative -! uses Samar's radiation formula -! Eqn: rhoc*T_t = (k*T_z)_z -! BC (z=0): Q(t) + kT_z = em*sig*T^4 -! BC (z=L): heat flux = Fgeotherm -! -! nz = number of grid points -! dt = time step -! Qn,Qnp1 = net solar insolation at time steps n and n+1 [W/m^2] -! T = vertical temperature profile [K] (in- and output) -! ti = thermal inertia [J m^-2 K^-1 s^-1/2] VECTOR -! rhoc = rho*c VECTOR where rho=density [kg/m^3] and -! c=specific heat [J K^-1 kg^-1] -! ti and rhoc are not allowed to vary in the layers immediately -! adjacent to the surface or the bottom -! emiss = emissivity -! Tsurf = surface temperature [K] (in- and output) -! Fgeotherm = geothermal heat flux at bottom boundary [W/m^2] -! Fsurf = heat flux at surface [W/m^2] (output) -! -! Grid: surface is at z=0 -! z(0)=0, z(2)=3*z(1), i.e., the top layer has half the width -! T(1) is at z(1); ...; T(i) is at z(i) -! k(i), rhoc(i), ti(i) are midway between z(i-1) and z(i) -! -! originally written by Samar Khatiwala, 2001 -! extended to variable thermal properties -! and irregular grid by Norbert Schorghofer -! added predictor-corrector 9/2019 -!*********************************************************************** - - implicit none - real*8, parameter :: sigSB=5.6704d-8 - - integer, intent(IN) :: nz - real*8, intent(IN) :: z(nz), dt, Qn, Qnp1, ti(nz),rhoc(nz) - real*8, intent(IN) :: emiss, Fgeotherm - real*8, intent(INOUT) :: T(nz), Tsurf - real*8, intent(OUT) :: Fsurf - integer i, iter - real*8 k(nz), k1, alpha(nz), gamma(nz), Tr - real*8 a(nz), b(nz), c(nz), r(nz), Told(nz) - real*8 arad, brad, ann, annp1, bn, buf, dz, beta - - ! set some constants - k(:) = ti(:)**2/rhoc(:) ! thermal conductivity - dz = 2.*z(1) - beta = dt/rhoc(1)/(2.*dz**2) ! assumes rhoc(0)=rhoc(1) - alpha(1) = beta*k(2) - gamma(1) = beta*k(1) - do i=2,nz-1 - buf = dt/(z(i+1)-z(i-1)) - alpha(i) = 2.*k(i+1)*buf/(rhoc(i)+rhoc(i+1))/(z(i+1)-z(i)) - gamma(i) = 2.*k(i)*buf/(rhoc(i)+rhoc(i+1))/(z(i)-z(i-1)) - enddo - buf = dt/(z(nz)-z(nz-1))**2 - gamma(nz) = k(nz)*buf/(2*rhoc(nz)) ! assumes rhoc(nz+1)=rhoc(nz) - - k1 = k(1)/dz - - ! elements of tridiagonal matrix - a(:) = -gamma(:) ! a(1) is not used - b(:) = 1. + alpha(:) + gamma(:) ! b(1) has to be reset at every timestep - c(:) = -alpha(:) ! c(nz) is not used - b(nz) = 1. + gamma(nz) - - Tr = Tsurf ! 'reference' temperature - iter = 0 - Told(:) = T(:) -30 continue ! in rare case of iterative correction - - ! Emission - arad = -3.*emiss*sigSB*Tr**4 - brad = 2.*emiss*sigSB*Tr**3 - ann = (Qn-arad)/(k1+brad) - annp1 = (Qnp1-arad)/(k1+brad) - bn = (k1-brad)/(k1+brad) - b(1) = 1. + alpha(1) + gamma(1) - gamma(1)*bn - - ! Set RHS - r(1) = gamma(1)*(annp1+ann) + & - & (1.-alpha(1)-gamma(1)+gamma(1)*bn)*T(1) + alpha(1)*T(2) - do concurrent (i=2:nz-1) - r(i) = gamma(i)*T(i-1) + (1.-alpha(i)-gamma(i))*T(i)+ alpha(i)*T(i+1) - enddo - r(nz) = gamma(nz)*T(nz-1) + (1.-gamma(nz))*T(nz) & - & + dt/rhoc(nz)*Fgeotherm/(z(nz)-z(nz-1)) ! assumes rhoc(nz+1)=rhoc(nz) - - ! Solve for T at n+1 - call tridag(a,b,c,r,T,nz) ! update by tridiagonal inversion - - Tsurf = 0.5*(annp1 + bn*T(1) + T(1)) ! (T0+T1)/2 - - ! iterative predictor-corrector - if ((Tsurf > 1.2*Tr .or. Tsurf < 0.8*Tr) .and. iter<10) then ! linearization error expected - ! redo until Tr is within 20% of new surface temperature - ! (under most circumstances, the 20% threshold is never exceeded) - iter = iter+1 - Tr = sqrt(Tr*Tsurf) ! linearize around an intermediate temperature - T(:) = Told(:) - goto 30 - endif - !if (iter>=5) print *,'consider taking shorter time steps',iter - !if (iter>=10) print *,'warning: too many iterations' - - Fsurf = - k(1)*(T(1)-Tsurf)/z(1) ! heat flux into surface -end subroutine conductionQ Index: trunk/LMDZ.COMMON/libf/evolution/conductionT.f90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/conductionT.f90 (revision 3492) +++ (revision ) @@ -1,70 +1,0 @@ -subroutine conductionT(nz,z,dt,T,Tsurf,Tsurfp1,ti,rhoc,Fgeotherm,Fsurf) -!*********************************************************************** -! conductionT: program to calculate the diffusion of temperature -! into the ground with prescribed surface temperature -! and variable thermal properties on irregular grid -! Crank-Nicholson scheme, flux conservative -! -! Eqn: rhoc*T_t = (k*T_z)_z -! BC (z=0): T=T(t) -! BC (z=L): heat flux = Fgeotherm -! -! nz = number of grid points -! dt = time step -! T = vertical temperature profile [K] (in- and output) -! Tsurf, Tsurfp1 = surface temperatures at times n and n+1 -! ti = thermal inertia [J m^-2 K^-1 s^-1/2] VECTOR -! rhoc = rho*c heat capacity per volume [J m^-3 K^-1] VECTOR -! ti and rhoc are not allowed to vary in the layers immediately -! adjacent to the surface or the bottom -! Fgeotherm = geothermal heat flux at bottom boundary [W/m^2] -! Fsurf = heat flux at surface [W/m^2] (output) -! -! Grid: surface is at z=0 -! T(1) is at z(1); ...; T(i) is at z(i) -! k(i) is midway between z(i-1) and z(i) -! rhoc(i) is midway between z(i-1) and z(i) -!*********************************************************************** - implicit none - - integer, intent(IN) :: nz - real*8, intent(IN) :: z(nz), dt, Tsurf, Tsurfp1, ti(nz), rhoc(nz) - real*8, intent(IN) :: Fgeotherm - real*8, intent(INOUT) :: T(nz) - real*8, intent(OUT) :: Fsurf - integer i - real*8 alpha(nz), k(nz), gamma(nz), buf - real*8 a(nz), b(nz), c(nz), r(nz) - - ! set some constants - k(:) = ti(:)**2/rhoc(:) ! thermal conductivity - alpha(1) = k(2)*dt/rhoc(1)/(z(2)-z(1))/z(2) - gamma(1) = k(1)*dt/rhoc(1)/z(1)/z(2) - do i=2,nz-1 - buf = dt/(z(i+1)-z(i-1)) - alpha(i) = k(i+1)*buf*2./(rhoc(i)+rhoc(i+1))/(z(i+1)-z(i)) - gamma(i) = k(i)*buf*2./(rhoc(i)+rhoc(i+1))/(z(i)-z(i-1)) - enddo - buf = dt/(z(nz)-z(nz-1))**2 - gamma(nz) = k(nz)*buf/(rhoc(nz)+rhoc(nz)) ! assumes rhoc(nz+1)=rhoc(nz) - - ! elements of tridiagonal matrix - a(:) = -gamma(:) ! a(1) is not used - b(:) = 1. + alpha(:) + gamma(:) - c(:) = -alpha(:) ! c(nz) is not used - b(nz) = 1. + gamma(nz) - - ! Set RHS - r(1)= alpha(1)*T(2) + (1.-alpha(1)-gamma(1))*T(1) + gamma(1)*(Tsurf+Tsurfp1) - do concurrent (i=2:nz-1) - r(i) = gamma(i)*T(i-1) + (1.-alpha(i)-gamma(i))*T(i) + alpha(i)*T(i+1) - enddo - r(nz) = gamma(nz)*T(nz-1) + (1.-gamma(nz))*T(nz) + & - & dt/rhoc(nz)*Fgeotherm/(z(nz)-z(nz-1)) ! assumes rhoc(nz+1)=rhoc(nz) - - ! Solve for T at n+1 - call tridag(a,b,c,r,T,nz) ! update by tridiagonal inversion - - Fsurf = -k(1)*(T(1)-Tsurfp1)/z(1) ! heat flux into surface - -end subroutine conductionT Index: trunk/LMDZ.COMMON/libf/evolution/derivs.f90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/derivs.f90 (revision 3492) +++ (revision ) @@ -1,115 +1,0 @@ - -subroutine deriv1(z,nz,y,y0,yNp1,yp) - ! first derivative of a function y(z) on irregular grid - ! upper b.c.: y(0)=y0 - ! lower b.c.: yp =0 - implicit none - integer, intent(IN) :: nz - real(8), intent(IN) :: z(nz),y(nz),y0,yNp1 - real(8), intent(OUT) :: yp(nz) - integer j - real(8) hm,hp,c1,c2,c3 - - hp = z(2)-z(1) - !hm = z(1) - c1 = z(1)/(hp*z(2)) - c2 = 1/z(1) - 1/(z(2)-z(1)) - c3 = -hp/(z(1)*z(2)) - yp(1) = c1*y(2) + c2*y(1) + c3*y0 - do j=2,nz-1 - hp = z(j+1)-z(j) - hm = z(j)-z(j-1) - c1 = +hm/(hp*(z(j+1)-z(j-1))) - c2 = 1/hm - 1/hp - c3 = -hp/(hm*(z(j+1)-z(j-1))) - yp(j) = c1*y(j+1) + c2*y(j) + c3*y(j-1) - enddo - yp(nz) = (yNp1 - y(nz-1))/(2.*(z(nz)-z(nz-1))) -end subroutine deriv1 - - - -subroutine deriv2_full(z,nz,a,b,a0,b0,bNp1,yp2) - ! second derivative (a*b_z)_z on irregular grid - ! upper b.c.: a(0)=a0, b(0)=b0 - ! 2nd order, without cross-terms - implicit none - integer, intent(IN) :: nz - real(8), intent(IN) :: z(nz),a(nz),b(nz),a0,b0,bNp1 - real(8), intent(OUT) :: yp2(nz) - integer j - real(8) hm,hp,c1,c2,c3 - - c1 = 1./(z(1)*z(2)) - c2 = 1./((z(2)-z(1))*z(1)) - c3 = 1./((z(2)-z(1))*z(2)) - yp2(1) = -c2*a(1)*b(1) & - & -c1*a0*b(1) + c1*(a(1)+a0)*b0 & - & -c3*a(2)*b(1) + c3*(a(1)+a(2))*b(2) - do j=2,nz-1 - hp = z(j+1)-z(j) - hm = z(j)-z(j-1) - c1 = 1./(hm*(z(j+1)-z(j-1))) - c2 = 1./(hp*hm) - c3 = 1./(hp*(z(j+1)-z(j-1))) - yp2(j) = -c2*a(j)*b(j) & - & -c1*a(j-1)*b(j) + c1*(a(j)+a(j-1))*b(j-1) & - & -c3*a(j+1)*b(j) + c3*(a(j)+a(j+1))*b(j+1) - enddo - - ! lower bc: b_z = 0 - !yp(nz)= 2.*a(nz)*(b(nz-1)-b(nz))/(z(nz)-z(nz-1))**2 - !more general lower bc - yp2(nz) = a(nz)*(bNp1 - 2*b(nz) + b(nz-1))/(z(nz)-z(nz-1))**2 -end subroutine deriv2_full - - - -subroutine deriv2_simple(z,nz,y,y0,yNp1,yp2) - ! second derivative y_zz on irregular grid - ! boundary conditions: y(0)=y0, y(nz+1)=yNp1 - implicit none - integer, intent(IN) :: nz - real(8), intent(IN) :: z(nz),y(nz),y0,yNp1 - real(8), intent(OUT) :: yp2(nz) - integer j - real(8) hm,hp,c1,c2,c3 - - c1 = +2./((z(2)-z(1))*z(2)) - c2 = -2./((z(2)-z(1))*z(1)) - c3 = +2./(z(1)*z(2)) - yp2(1) = c1*y(2) + c2*y(1) + c3*y0 - - do j=2,nz-1 - hp = z(j+1)-z(j) - hm = z(j)-z(j-1) - c1 = +2./(hp*(z(j+1)-z(j-1))) - c2 = -2./(hp*hm) - c3 = +2./(hm*(z(j+1)-z(j-1))) - yp2(j) = c1*y(j+1) + c2*y(j) + c3*y(j-1) - enddo - - yp2(nz) = (yNp1 - 2*y(nz) + y(nz-1))/(z(nz)-z(nz-1))**2 -end subroutine deriv2_simple - - - -real(8) function deriv1_onesided(j,z,nz,y) - ! first derivative of function y(z) at z(j) - ! one-sided derivative on irregular grid - implicit none - integer, intent(IN) :: nz,j - real(8), intent(IN) :: z(nz),y(nz) - real(8) h1,h2,c1,c2,c3 - if (j<1 .or. j>nz-2) then - deriv1_onesided = -9999. - else - h1 = z(j+1)-z(j) - h2 = z(j+2)-z(j+1) - c1 = -(2*h1+h2)/(h1*(h1+h2)) - c2 = (h1+h2)/(h1*h2) - c3 = -h1/(h2*(h1+h2)) - deriv1_onesided = c1*y(j) + c2*y(j+1) + c3*y(j+2) - endif -end function deriv1_onesided - Index: trunk/LMDZ.COMMON/libf/evolution/dyn_ss_ice_m.f90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/dyn_ss_ice_m.f90 (revision 3492) +++ (revision ) @@ -1,177 +1,0 @@ - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!! -!!! Purpose: Retreat and growth of subsurface ice on Mars -!!! orbital elements remain constant -!!! -!!! -!!! Author: EV, updated NS MSIM dynamical program for the PEM -!!! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - - -SUBROUTINE dyn_ss_ice_m(ssi_depth_in,T1,T_in,nsoil,thIn,p0,pfrost,porefill_in,porefill,ssi_depth) - -!*********************************************************************** -! Retreat and growth of subsurface ice on Mars -! orbital elements remain constant -!*********************************************************************** - use miscparameters, only : pi, d2r, NMAX, marsDay, solsperyear - !use allinterfaces - implicit none - integer, parameter :: NP=1 ! # of sites - integer nz, i, k, nsoil, iloop - real(8) zmax, delta, z(NMAX), icetime, porosity, icefrac - real(8), dimension(NP) :: albedo, thIn, rhoc - real(8), dimension(NP) :: pfrost, p0 - real(8) newti, stretch, newrhoc, ecc, omega, eps, timestep - real(8) ssi_depth_in, ssi_depth, T1 - real(8), dimension(NP) :: zdepthF, zdepthE, zdepthT, zdepthG - real(8), dimension(NMAX,NP) :: porefill, porefill_in - real(8), dimension(NP) :: Tb, T_in, Tmean1, Tmean3, avrho1 - real(8) tmax, tlast, avrho1prescribed(NP), l1 - real(8), external :: smartzfac - - !if (iargc() /= 1) then - ! stop 'USAGE: icages ext' - !endif - !call getarg( 1, ext ) - - if (NP>100) stop 'subroutine icelayer_mars cannot handle this many sites' - - ! parameters that never ever change - nz=nsoil - porosity = 0.4d0 ! porosity of till - !rhoc(:) = 1500.*800. ! will be overwritten - icefrac = 0.98 - tmax = 1 - tlast = 0. - avrho1prescribed(:) = pfrost/T1 ! <0 means absent - albedo=0.23 - !avrho1prescribed(:) = 0.16/200. ! units are Pa/K - - !open(unit=21,file='lats.'//ext,action='read',status='old',iostat=ierr) - !if (ierr /= 0) then - ! print *,'File lats.'//ext,'not found' - ! stop - !endif - do k=1,NP - !read(21,*) latitude(k),albedo(k),thIn(k),htopo(k) - ! empirical relation from Mellon & Jakosky - rhoc(k) = 800.*(150.+100.*sqrt(34.2+0.714*thIn(k))) - enddo - !close(21) - - ! set eternal grid - zmax = 25. - !zfac = smartzfac(nz,zmax,6,0.032d0) - !call setgrid(nz,z,zmax,zfac) - l1=2.e-4 - do iloop=0,nsoil - z(iloop) = l1*(1+iloop**2.9*(1-exp(-real(iloop)/20.))) - enddo - - - !open(unit=30,file='z.'//ext,action='write',status='unknown') - !write(30,'(999(f8.5,1x))') z(1:nz) - !close(30) - - !ecc = ecc_in; eps = obl_in*d2r; omega = Lp_in*d2r ! today - ! total atmospheric pressure - !p0(:) = 600. - ! presently 520 Pa at zero elevation (Smith & Zuber, 1998) - ! do k=1,NP - ! p0(k)=520*exp(-htopo(k)/10800.) - ! enddo - timestep = 1 ! must be integer fraction of 1 ka - icetime = -tmax-timestep ! earth years - - ! initializations - Tb=T_in - zdepthF(:) = -9999. - - !zdepthT(1:NP) = -9999. ! reset again below - ! zdepthT(1:NP) = 0. - - ! print *,'RUNNING MARS_FAST' - ! print *,'Global model parameters:' - ! print *,'nz=',nz,' zfac=',zfac,'zmax=',zmax - ! print *,'porosity=',porosity - ! print *,'starting at time',icetime,'years' - ! print *,'time step=',timestep,'years' - ! print *,'eps=',eps/d2r,'ecc=',ecc,'omega=',omega/d2r - ! print *,'number of sites=',NP - ! print *,'Site specific parameters:' - do k=1,NP - if (NP>1) print *,' Site ',k - ! print *,' latitude (deg)',latitude(k),' rho*c (J/m^3/K)',rhoc(k),' thIn=',thIn(k) - ! print *,' total pressure=',p0(k),'partial pressure=',pfrost(k) - delta = thIn(k)/rhoc(k)*sqrt(marsDay/pi) - ! print *,' skin depths (m)',delta,delta*sqrt(solsperyear) - call soilthprop(porosity,1.d0,rhoc(k),thIn(k),1,newrhoc,newti,icefrac) - stretch = (newti/thIn(k))*(rhoc(k)/newrhoc) - do i=1,nz - if (z(i)=tlast) exit - enddo - - ! close(34) - ! close(36); close(37) - -!501 format (f10.0,2x,f7.3,2x,f10.4,2(2x,f6.2),2x,f9.3,2x,g11.4) - -end subroutine dyn_ss_ice_m - - - Index: trunk/LMDZ.COMMON/libf/evolution/dyn_ss_ice_m_wrapper.f90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/dyn_ss_ice_m_wrapper.f90 (revision 3492) +++ (revision ) @@ -1,27 +1,0 @@ -subroutine dyn_ss_ice_m_wrapper(ngrid,nsoil,tHIn,p0,pfrost,T_in,ssi_depth_in,porefill_in,porefill,ssi_depth) - - use allinterfaces, only: dyn_ss_ice_m - implicit none - integer ngrid,nsoil - integer i - real(8), dimension(ngrid) :: p0, pfrost - real(8), dimension(ngrid) :: thIn, ssi_depth_in,T1,ssi_depth - real(8), dimension(nsoil,ngrid) :: T_in - real(8), dimension(nsoil,ngrid) :: porefill, porefill_in - - do i=1,ngrid - !tHIn(i)=250 - !p0(i)=611 - !pfrost(i)=0.31 - !T_in(1:nsoil,i)=180 - !ssi_depth_in(i)=0.3*i - porefill(1:nsoil,i) = -9999 - ssi_depth(i)=-9999 - !porefill_in(1:nsoil,i) =0.01*i - T1(i)=T_in(1,i) - call dyn_ss_ice_m(ssi_depth_in(i),T1(i),T_in(:,i),nsoil,thIn(i),p0(i),pfrost(i),porefill_in(:,i),porefill(:,i),ssi_depth(i)) - - !print *,'new porefill =',porefill(:,i) - enddo - -end subroutine dyn_ss_ice_m_wrapper Index: trunk/LMDZ.COMMON/libf/evolution/fast_modules.f90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/fast_modules.f90 (revision 3492) +++ (revision ) @@ -1,265 +1,0 @@ - -module allinterfaces - ! interfaces from Fortran 90 subroutines and functions - ! comments have been stripped - - !begin fast_subs_mars.f90 - - interface - subroutine icelayer_mars(bigstep,nz,NP,thIn,rhoc,z,porosity,pfrost, & - & Tb,zdepthF,zdepthE,porefill,Tmean1,Tmean3,zdepthG, & - & albedo,p0,icefrac,zdepthT,avrho1, & - & avrho1prescribed) - use miscparameters, only : d2r, NMAX, icedensity - implicit none - integer, intent(IN) :: nz, NP - real(8), intent(IN) :: bigstep - real(8), intent(IN), dimension(NP) :: thIn, rhoc - real(8), intent(IN) :: z(NMAX) - real(8), intent(IN) :: porosity, pfrost(NP) - real(8), intent(INOUT) :: Tb(NP), porefill(nz,NP), zdepthF(NP), zdepthT(NP) - real(8), intent(OUT), dimension(NP) :: zdepthE, Tmean1, Tmean3, zdepthG - real(8), intent(IN), dimension(NP) :: albedo, p0 - real(8), intent(IN) :: icefrac - real(8), intent(OUT) :: avrho1(NP) - real(8), intent(IN), optional :: avrho1prescribed(NP) - end subroutine icelayer_mars - end interface - - interface - subroutine ajsub_mars(typeT, albedo0, pfrost, nz, z, ti, rhocv, & - & fracIR, fracDust, p0, avdrho, avdrhoP, avrho1, & - & Tb, zdepthE, typeF, zdepthF, ypp, porefill, Tmean1, Tmean3, & - & B, typeG, zdepthG, avrho1prescribed) - use miscparameters - implicit none - integer, intent(IN) :: nz, typeT - real(8), intent(IN) :: albedo0, pfrost, z(NMAX) - real(8), intent(IN) :: ti(NMAX), rhocv(NMAX), fracIR, fracDust, p0 - real(8), intent(IN) :: porefill(nz) - real(8), intent(OUT) :: avdrho, avdrhoP - real(8), intent(OUT) :: avrho1 - real(8), intent(INOUT) :: Tb, Tmean1 - integer, intent(OUT) :: typeF - real(8), intent(OUT) :: zdepthE, zdepthF - real(8), intent(OUT) :: ypp(nz) - real(8), intent(OUT) :: Tmean3, zdepthG - real(8), intent(IN) :: B - integer, intent(OUT) :: typeG - real(8), intent(IN), optional :: avrho1prescribed - end subroutine ajsub_mars - end interface - - !end of fast_subs_mars.f90 - !begin fast_subs_univ.f90 - - interface - pure function zint(y1,y2,z1,z2) - implicit none - real(8), intent(IN) :: y1,y2,z1,z2 - real(8) zint - end function zint - end interface - - interface - pure function equildepth(nz, z, rhosatav, rhosatav0, avrho1) - implicit none - integer, intent(IN) :: nz - real(8), intent(IN) :: z(nz), rhosatav(nz) - real(8), intent(IN) :: rhosatav0, avrho1 - real(8) zint, equildepth - external zint - end function equildepth - end interface - - interface - subroutine depths_avmeth(typeT, nz, z, rhosatav, rhosatav0, rlow, avrho1, & - & porefill, typeF, zdepthF, B, ypp, typeG, zdepthG) - use miscparameters, only : icedensity - implicit none - integer, intent(IN) :: nz, typeT - real(8), intent(IN), dimension(nz) :: z, rhosatav, porefill - real(8), intent(IN) :: rhosatav0, rlow, avrho1 - integer, intent(OUT) :: typeF - real(8), intent(INOUT) :: zdepthF - real(8), intent(IN) :: B - real(8), intent(OUT) :: ypp(nz), zdepthG - integer, intent(INOUT) :: typeG - real(8), external :: zint, deriv1_onesided, colint - end subroutine depths_avmeth - end interface - - interface - pure function constriction(porefill) - implicit none - real(8), intent(IN) :: porefill - real(8) constriction - end function constriction - end interface - - interface - subroutine assignthermalproperties(nz,thIn,rhoc,ti,rhocv, & - & typeT,icefrac,porosity,porefill) - implicit none - integer, intent(IN) :: nz - integer, intent(IN), optional :: typeT - real(8), intent(IN), optional :: icefrac - real(8), intent(IN) :: thIn, rhoc - real(8), intent(IN), optional :: porosity, porefill(nz) - real(8), intent(OUT) :: ti(nz), rhocv(nz) - end subroutine assignthermalproperties - end interface - - interface - pure subroutine icechanges_poreonly(nz,z,typeF,typeG,avdrhoP,ypp,B,porefill) - implicit none - integer, intent(IN) :: nz, typeF, typeG - real(8), intent(IN) :: z(nz), ypp(nz), avdrhoP, B - real(8), intent(INOUT) :: porefill(nz) - end subroutine icechanges_poreonly - end interface - - interface - pure subroutine icechanges(nz,z,typeF,avdrho,avdrhoP,ypp, & - & Diff,porosity,icefrac,bigstep,zdepthT,porefill,typeG) - implicit none - integer, intent(IN) :: nz, typeF, typeG - real(8), intent(IN) :: z(nz), ypp(nz), avdrho, avdrhoP - real(8), intent(IN) :: Diff, porosity, icefrac, bigstep - real(8), intent(INOUT) :: zdepthT, porefill(nz) - end subroutine icechanges - end interface - - interface - subroutine compactoutput(unit,porefill,nz) - implicit none - integer, intent(IN) :: unit,nz - real(8), intent(IN) :: porefill(nz) - end subroutine compactoutput - end interface - - !end of fast_subs_univ - !begin derivs.f90 - - interface - subroutine deriv1(z,nz,y,y0,yNp1,yp) - implicit none - integer, intent(IN) :: nz - real(8), intent(IN) :: z(nz),y(nz),y0,yNp1 - real(8), intent(OUT) :: yp(nz) - end subroutine deriv1 - end interface - - interface - subroutine deriv2_full(z,nz,a,b,a0,b0,bNp1,yp2) - implicit none - integer, intent(IN) :: nz - real(8), intent(IN) :: z(nz),a(nz),b(nz),a0,b0,bNp1 - real(8), intent(OUT) :: yp2(nz) - end subroutine deriv2_full - end interface - - interface - subroutine deriv2_simple(z,nz,y,y0,yNp1,yp2) - implicit none - integer, intent(IN) :: nz - real(8), intent(IN) :: z(nz),y(nz),y0,yNp1 - real(8), intent(OUT) :: yp2(nz) - end subroutine deriv2_simple - end interface - interface - real(8) pure function deriv1_onesided(j,z,nz,y) - implicit none - integer, intent(IN) :: nz,j - real(8), intent(IN) :: z(nz),y(nz) - end function deriv1_onesided - end interface - - !end of derivs.f90 - !begin fast_subs_exper.f90 - - interface - subroutine icelayer_exper(bigstep, nz, thIn, rhoc, z, porosity, pfrost, & - & zdepthT, icefrac, zdepthF, zdepthE, porefill, Tmean, Tampl, Diff, zdepthG) - use miscparameters, only : d2r, NMAX, icedensity - implicit none - integer, intent(IN) :: nz - real(8), intent(IN) :: bigstep - real(8), intent(IN) :: thIn, rhoc, z(NMAX), porosity, pfrost - real(8), intent(INOUT) :: zdepthT, zdepthF, porefill(nz) - real(8), intent(OUT) :: zdepthE - real(8), intent(IN) :: icefrac, Diff, Tmean, Tampl - real(8), intent(OUT) :: zdepthG - end subroutine icelayer_exper - end interface - - interface - subroutine ajsub_exper(typeT, nz, z, ti, rhocv, pfrost, Tmean, Tampl, & - & avdrho, avdrhoP, zdepthE, typeF, zdepthF, ypp, porefill, & - & B, typeG, zdepthG) - use miscparameters, only : NMAX, solsperyear, marsDay - implicit none - integer, intent(IN) :: nz, typeT - real(8), intent(IN) :: z(NMAX), ti(NMAX), rhocv(NMAX), pfrost - real(8), intent(IN) :: Tmean, Tampl - real(8), intent(OUT) :: avdrho, avdrhoP - real(8), intent(OUT) :: zdepthE - integer, intent(OUT) :: typeF - real(8), intent(INOUT) :: zdepthF - real(8), intent(OUT) :: ypp(nz) - real(8), intent(IN) :: porefill(nz) - real(8), intent(IN) :: B - integer, intent(OUT) :: typeG - real(8), intent(OUT) :: zdepthG - real(8), external :: Tsurface, psv - real(8), external :: equildepth - end subroutine ajsub_exper - end interface - - !end fast_subs_exper - - ! Other - interface - pure function flux_mars77(R,decl,HA,albedo,fracir,fracscat) - implicit none - real*8 flux_mars77 - real*8, intent(IN) :: R,decl,HA,albedo,fracIR,fracScat - end function flux_mars77 - end interface - - interface - pure function colint(y,z,nz,i1,i2) - implicit none - integer, intent(IN) :: nz, i1, i2 - real(8), intent(IN) :: y(nz),z(nz) - real(8) colint - end function colint - end interface - - interface - subroutine dyn_ss_ice_m(ssi_depth_in,T1,T_in,nsoil,thIn,p0,pfrost,porefill_in,porefill,ssi_depth) - implicit none - integer, intent(IN) :: nsoil - real(8), intent(IN) :: thIn,ssi_depth_in,T1 - real(8), intent(IN) :: p0(1), pfrost(1) - real(8), intent(IN) :: T_in(nsoil) - real(8), intent(INOUT) :: porefill(nsoil,1) - real(8), intent(INOUT) :: porefill_in(nsoil,1) - real(8), intent(INOUT) :: ssi_depth - end subroutine dyn_ss_ice_m - end interface - - interface - subroutine dyn_ss_ice_m_wrapper(ngrid,nsoil,tHIn,p0,pfrost,T_in,ssi_depth_in,porefill_in,porefill,ssi_depth) - implicit none - integer, intent(IN) :: nsoil,ngrid - real(8), intent(IN) :: thIn(ngrid),ssi_depth_in(ngrid) - real(8), intent(IN) :: p0(ngrid), pfrost(ngrid) - real(8), intent(IN) :: T_in(nsoil,ngrid) - real(8), intent(OUT) :: porefill(nsoil,ngrid) - real(8), intent(IN) :: porefill_in(nsoil,ngrid) - real(8), intent(OUT) :: ssi_depth(ngrid) - end subroutine dyn_ss_ice_m_wrapper - end interface - -end module allinterfaces Index: trunk/LMDZ.COMMON/libf/evolution/fast_subs_mars.f90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/fast_subs_mars.f90 (revision 3492) +++ (revision ) @@ -1,307 +1,0 @@ -module thermalmodelparam_mars - ! parameters for thermal model - ! they are only used in the subroutines below - real(8), parameter :: dt = 0.02 ! in units of Mars solar days - !real(8), parameter :: Fgeotherm = 0. - real(8), parameter :: Fgeotherm = 0 !0.028 ! [W/m^2] - real(8), parameter :: Lco2frost=6.0e5, co2albedo=0.60, co2emiss=1. - real(8), parameter :: emiss0 = 1. ! emissivity of dry surface - integer, parameter :: EQUILTIME =15 ! [Mars years] -end module thermalmodelparam_mars - - -!***************************************************** -! Subroutines for fast method -! written by Norbert Schorghofer 2007-2011 -!***************************************************** - - -subroutine icelayer_mars(bigstep,nz,NP,thIn,rhoc,z,porosity,pfrost, & - & Tb,zdepthF,zdepthE,porefill,Tmean1,Tmean3,zdepthG, & - & albedo,p0,icefrac,zdepthT,avrho1, & - & avrho1prescribed) -!************************************************************************* -! bigstep = time step [Earth years] -! latitude [degree] -!************************************************************************* - use miscparameters, only : d2r, NMAX, icedensity - use allinterfaces, except_this_one => icelayer_mars - !use omp_lib - implicit none - integer, intent(IN) :: nz, NP - real(8), intent(IN) :: bigstep - real(8), intent(IN) :: thIn(NP), rhoc(NP), z(NMAX), porosity, pfrost(NP) - real(8), intent(INOUT) :: Tb(NP), porefill(nz,NP), zdepthF(NP), zdepthT(NP) - real(8), intent(OUT), dimension(NP) :: zdepthE, Tmean1, Tmean3, zdepthG - real(8), intent(IN), dimension(NP) :: albedo, p0 - real(8), intent(IN) :: icefrac - real(8), intent(OUT) :: avrho1(NP) - real(8), intent(IN), optional :: avrho1prescribed(NP) ! <0 means absent - - integer k, typeF, typeG, typeT, j, jump, typeP - real(8) fracIR, fracDust, ti(NMAX), rhocv(NMAX) - real(8) Diff, ypp(nz), avdrho(NP), avdrhoP(NP), B, deltaz - real(8), SAVE :: avdrho_old(100), zdepth_old(100) ! NP<=100 - logical mode2 - - !$omp parallel & - !$omp private(Diff,fracIR,fracDust,B,typeT,j,ti,rhocv,typeF,jump,typeG) - !$omp do - do k=1,NP ! big loop - - Diff = 4e-4*600./p0(k) - fracIR = 0.04*p0(k)/600.; fracDust = 0.02*p0(k)/600. - B = Diff*bigstep*86400.*365.24/(porosity*icedensity) - - typeT = -9 - if (zdepthT(k)>=0. .and. zdepthT(k)zdepthT(k)) then ! ice - typeT = j ! first point in ice - exit - endif - enddo - endif - - ! call assignthermalproperties(nz,thIn(k),rhoc(k),ti,rhocv, & - ! & typeT,icefrac,porosity,porefill(:,k)) - - !----run thermal model - call ajsub_mars(typeT, albedo(k), pfrost(k), nz, z, & - & ti, rhocv, fracIR, fracDust, p0(k), & - & avdrho(k), avdrhoP(k), avrho1(k), Tb(k), zdepthE(k), typeF, & - & zdepthF(k), ypp, porefill(:,k), Tmean1(k), Tmean3(k), B, & - & typeG, zdepthG(k), avrho1prescribed(k)) - - if (typeF*zdepthF(k)<0.) stop 'error in icelayer_mars' - ! diagnose - if (zdepthT(k)>=0.) then - jump = 0 - do j=1,nz - if (zdepth_old(k)z(j)) jump = jump+1 - enddo - else - jump = -9 - endif - if (zdepthT(k)>=0. .and. avdrho(k)*avdrho_old(k)<0.) then - write(34,*) '# zdepth arrested' - if (jump>1) write(34,*) '# previous step was too large',jump - endif -! write(34,'(f8.3,1x,f6.2,1x,f11.5,2(1x,g11.4),1x,i3)') & -! & bigstep,latitude(k),zdepthT(k),avdrho(k),avdrhoP(k),jump - zdepth_old(k) = zdepthT(k) - avdrho_old(k) = avdrho(k) - -!----mode 2 growth - typeP = -9; mode2 = .FALSE. - do j=1,nz - if (porefill(j,k)>0.) then - typeP = j ! first point with ice - exit - endif - enddo - if (typeT>0 .and. typeP>2 .and. zdepthE(k)>0.) then - if (porefill(typeP,k)>=1. .and. porefill(typeP-1,k)==0. .and. & - & zdepthE(k)2*(z(typeP)-z(typeP-1))) then ! trick that avoids oscillations - deltaz = -avdrhoP(k)/z(typeP)*18./8314.*B ! conservation of mass - if (deltaz>z(typeP)-z(typeP-1)) then ! also implies avdrhoP<0. - mode2 = .TRUE. - endif - endif - endif - - !call icechanges_poreonly(nz,z,typeF,typeG,avdrhoP(k),ypp,B,porefill(:,k)) - call icechanges(nz,z(:),typeF,avdrho(k),avdrhoP(k),ypp(:), & - & Diff,porosity,icefrac,bigstep,zdepthT(k),porefill(:,k),typeG) - - if (mode2 .and. porefill(typeP,k)>=1. .and. porefill(typeP-1,k)==0.) then ! nothing changed - porefill(typeP-1,k)=1. ! paste a layer - ! write(34,*) '# mode 2 growth occurred',typeP,typeF,typeT - endif - - enddo ! end of big loop - !$omp end do - !$omp end parallel -end subroutine icelayer_mars - - - -subroutine ajsub_mars(typeT, albedo0, pfrost, nz, z, ti, rhocv, & - & fracIR, fracDust, patm, avdrho, avdrhoP, avrho1, & - & Tb, zdepthE, typeF, zdepthF, ypp, porefill, Tmean1, Tmean3, & - & B, typeG, zdepthG, avrho1prescribed) -!*********************************************************************** -! A 1D thermal model that returns various averaged quantities -! -! Tb<0 initializes temperatures -! Tb>0 initializes temperatures with Tb -!*********************************************************************** - use miscparameters - use thermalmodelparam_mars - use allinterfaces, except_this_one => ajsub_mars - implicit none - integer, intent(IN) :: nz, typeT -! real(8), intent(IN) :: latitude ! in radians - real(8), intent(IN) :: albedo0, pfrost, z(NMAX) - real(8), intent(IN) :: ti(NMAX), rhocv(NMAX), fracIR, fracDust, patm - real(8), intent(IN) :: porefill(nz) - real(8), intent(OUT) :: avdrho, avdrhoP ! difference in annual mean vapor density - real(8), intent(OUT) :: avrho1 ! mean vapor density on surface - real(8), intent(INOUT) :: Tb, Tmean1 - integer, intent(OUT) :: typeF ! index of depth below which filling occurs - real(8), intent(OUT) :: zdepthE, zdepthF - real(8), intent(OUT) :: ypp(nz) ! (d rho/dt)/Diff - real(8), intent(OUT) :: Tmean3, zdepthG - real(8), intent(IN) :: B ! just passing through - integer, intent(OUT) :: typeG - real(8), intent(IN), optional :: avrho1prescribed - real(8), parameter :: a = 1.52366 ! Mars semimajor axis in A.U. - integer nsteps, n, i, nm, typeP - real(8) tmax, time, Qn, Qnp1, tdays - real(8) marsR, marsLs, marsDec, HA - real(8) Tsurf, Tco2frost, albedo, Fsurf, m, dE, emiss, T(NMAX) - real(8) Told(nz), Fsurfold, Tsurfold, Tmean0, avrho2 - real(8) rhosatav0, rhosatav(nz), rlow - real(8), external :: psv, tfrostco2 - - tmax = EQUILTIME*solsperyear - nsteps = int(tmax/dt) ! calculate total number of timesteps - - Tco2frost = tfrostco2(patm) - - !if (Tb<=0.) then ! initialize - !Tmean0 = 210.15 ! black-body temperature of planet - !Tmean0 = (589.*(1.-albedo0)*cos(latitude)/(pi*emiss0*sigSB))**0.25 ! better estimate - !Tmean0 = Tmean0-5. - !write(34,*) '# initialized with temperature estimate at',latitude/d2r,'of',Tmean0,'K' - !T(1:nz) = Tmean0 - !else - T(1:nz) = Tb - ! not so good when Fgeotherm is on - !endif - - albedo = albedo0 - emiss = emiss0 - do i=1,nz - if (T(i)Tco2frost) then - ! call conductionQ(nz,z,dt*marsDay,Qn,Qnp1,T,ti,rhocv,emiss, & - ! & Tsurf,Fgeotherm,Fsurf) - endif - if (Tsurf0.) then ! CO2 condensation - T(1:nz) = Told(1:nz) - !call conductionT(nz,z,dt*marsDay,T,Tsurfold,Tco2frost,ti, & - !& rhocv,Fgeotherm,Fsurf) - Tsurf = Tco2frost - ! dE = (- Qn - Qnp1 + Fsurfold + Fsurf + & - ! & emiss*sigSB*(Tsurfold**4+Tsurf**4)/2. - m = m + dt*marsDay*dE/Lco2frost - endif - if (Tsurf>Tco2frost .or. m<=0.) then - albedo = albedo0 - emiss = emiss0 - else - albedo = co2albedo - emiss = co2emiss - endif - !Qn=Qnp1 - - if (time>=tmax-solsperyear) then - Tmean1 = Tmean1 + Tsurf - Tmean3 = Tmean3 + T(nz) - avrho1 = avrho1 + min(psv(Tsurf),pfrost)/Tsurf - rhosatav0 = rhosatav0 + psv(Tsurf)/Tsurf - do i=1,nz - rhosatav(i) = rhosatav(i) + psv(T(i))/T(i) - enddo - nm = nm+1 - endif - - enddo ! end of time loop - - avrho1 = avrho1/nm - if (present(avrho1prescribed)) then - if (avrho1prescribed>=0.) avrho1=avrho1prescribed - endif - Tmean1 = Tmean1/nm; Tmean3 = Tmean3/nm - rhosatav0 = rhosatav0/nm; rhosatav(:)=rhosatav(:)/nm - if (typeT>0) then - avrho2 = rhosatav(typeT) - else - avrho2 = rhosatav(nz) ! no ice - endif - avdrho = avrho2-avrho1 - typeP = -9 - do i=1,nz - if (porefill(i)>0.) then - typeP = i ! first point with ice - exit - endif - enddo - if (typeP>0) then - avdrhoP = rhosatav(typeP) - avrho1 - else - avdrhoP = -9999. - end if - - zdepthE = equildepth(nz, z, rhosatav, rhosatav0, avrho1) - - if (Fgeotherm>0.) then - Tb = Tmean1 - typeG = 1 ! will be overwritten by depths_avmeth - rlow = 2*rhosatav(nz)-rhosatav(nz-1) - else - Tb = T(nz) - typeG = -9 - rlow = rhosatav(nz-1) - endif - call depths_avmeth(typeT, nz, z, rhosatav(:), rhosatav0, rlow, avrho1, & - & porefill(:), typeF, zdepthF, B, ypp(:), typeG, zdepthG) - -end subroutine ajsub_mars - - - -subroutine outputmoduleparameters - use miscparameters - use thermalmodelparam_mars - implicit none -! print *,'Parameters stored in modules' -! print *,' Ice bulk density',icedensity,'kg/m^3' -! print *,' dt=',dt,'Mars solar days' -! print *,' Fgeotherm=',Fgeotherm,'W/m^2' -! write(6,'(2x,a27,1x,f5.3)') 'Emissivity of dry surface=',emiss0 -! write(6,'(2x,a12,1x,f5.3,2x,a16,1x,f5.3)') 'CO2 albedo=',co2albedo,'CO2 emissivity=',co2emiss -! print *,' Thermal model equilibration time',EQUILTIME,'Mars years' -end subroutine outputmoduleparameters - - - Index: trunk/LMDZ.COMMON/libf/evolution/fast_subs_univ.f90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/fast_subs_univ.f90 (revision 3492) +++ (revision ) @@ -1,383 +1,0 @@ -!***************************************************** -! Commonly used subroutines for fast method -! written by Norbert Schorghofer 2007-2011 -!***************************************************** - -pure function zint(y1,y2,z1,z2) - ! interpolate between two values, y1*y2<0 - implicit none - real(8), intent(IN) :: y1,y2,z1,z2 - real(8) zint - zint = (y1*z2 - y2*z1)/(y1-y2) -end function zint - - - -pure function equildepth(nz, z, rhosatav, rhosatav0, avrho1) -!*********************************************************************** -! returns equilibrium depth for given ice content -! this is not the true (final) equilibrium depth -!*********************************************************************** - use allinterfaces, except_this_one => equildepth - implicit none - integer, intent(IN) :: nz - real(8), intent(IN) :: z(nz), rhosatav(nz) - real(8), intent(IN) :: rhosatav0, avrho1 - integer i, typeE - real(8) equildepth ! =zdepthE - !real(8), external :: zint ! defined in allinterfaces.mod - - typeE = -9; equildepth = -9999. - do i=1,nz - if (rhosatav(i) <= avrho1) then - typeE=i - exit - endif - enddo - if (typeE>1) then ! interpolate - equildepth=zint(avrho1-rhosatav(typeE-1),avrho1-rhosatav(typeE),z(typeE-1),z(typeE)) - end if - if (typeE==1) equildepth=zint(avrho1-rhosatav0,avrho1-rhosatav(1),0.d0,z(1)) - if (equildepth>z(nz)) equildepth=-9999. ! happens very rarely -end function equildepth - - - -subroutine depths_avmeth(typeT, nz, z, rhosatav, rhosatav0, rlow, avrho1, & - & porefill, typeF, zdepthF, B, ypp, typeG, zdepthG) -!*********************************************************************** -! returns interface depth and ypp -! also returns lower geothermally limited boundary, if applicable -! -! this is a crucial subroutine -! -! B = Diff*bigstep/(porosity*icedensity) [SI units] -!*********************************************************************** - use allinterfaces, except_this_one => depths_avmeth - implicit none - integer, intent(IN) :: nz, typeT - real(8), intent(IN), dimension(nz) :: z, rhosatav, porefill - real(8), intent(IN) :: rhosatav0, rlow, avrho1 - integer, intent(OUT) :: typeF ! index of depth below which filling occurs - real(8), intent(INOUT) :: zdepthF - real(8), intent(IN) :: B - real(8), intent(OUT) :: ypp(nz), zdepthG - integer, intent(INOUT) :: typeG ! positive on input when Fgeotherm>0 - - integer i, typeP, nlb, newtypeG - real(8) eta(nz), Jpump1, help(nz), yp(nz), zdepthFold, ap_one, ap(nz) - real(8) leak, cumfill, cumfillabove - - if (typeT<0) then - nlb = nz - do i=1,nz - eta(i) = constriction(porefill(i)) - enddo - else - !nlb = typeT-1 - nlb = typeT ! changed 2010-09-29 - do i=1,typeT-1 - eta(i) = constriction(porefill(i)) - enddo - do i=typeT,nz - eta(i)=0. - enddo - end if - -!-fill depth - zdepthFold = zdepthF - typeF = -9; zdepthF = -9999. - call deriv1(z,nz,rhosatav,rhosatav0,rlow,yp) ! yp also used below - do i=1,nlb - Jpump1 = (rhosatav(i)-avrho1)/z(i) ! <0 when stable - ! yp is always <0 - help(i) = Jpump1 - eta(i)*yp(i) - leak = porefill(i)/B*(z(i)-zdepthFold)/(18./8314.) - !help(i) = help(i)-leak ! optional - if (help(i) <= 0.) then - typeF=i - !print *,'#',typeF,Jpump1,eta(typeF)*yp(typeF),leak - exit - endif - enddo - if (typeF>1) zdepthF = zint(help(typeF-1),help(typeF),z(typeF-1),z(typeF)) - if (typeF==1) zdepthF=z(1) - - -!-depth to shallowest perennial ice - typeP = -9 - do i=1,nz - if (porefill(i)>0.) then - typeP = i ! first point with ice - exit - endif - enddo - -!-calculate ypp - !call deriv1(z,nz,rhosatav,rhosatav0,rlow,yp) - call deriv1(z,nz,eta(:),1.d0,eta(nz-1),ap) - if (typeP>0 .and. typeP0) typeG=-9 - if (typeG<0) zdepthG=-9999. - if (typeG>0 .and. typeT<0) then - typeG=-9 - do i=2,nz - if (yp(i)>0.) then ! first point with reversed flux - typeG=i - zdepthG=zint(yp(i-1),yp(i),z(i-1),z(i)) - !zdepthG=z(typeG) - exit - endif - enddo - else - typeG = -9 - endif - if (typeG>0 .and. typeT<0) then - cumfillabove = colint(porefill(:)/eta(:),z,nz,typeG-1,nz) - newtypeG = -9 - do i=typeG,nz - if (minval(eta(i:nz))<=0.) cycle ! eta=0 means completely full - cumfill=colint(porefill(:)/eta(:),z,nz,i,nz) - if (cumfilltypeG) then - write(34,*) '# adjustment to geotherm depth by',i-typeG - zdepthG = zint(yp(i-1)*18./8314.*B-cumfillabove, & ! no good - & yp(i)*18./8314.*B-cumfill,z(i-1),z(i)) - if (zdepthG>z(i) .or. zdepthG0) typeG=newtypeG - end if - ! if typeG>0, then all ice at and below typeG should be erased -end subroutine depths_avmeth - - - -pure function constriction(porefill) -! specify constriction function here, 0<=eta<=1 - implicit none - real(8), intent(IN) :: porefill - real(8) eta, constriction - if (porefill<=0.) eta = 1. - if (porefill>0. .and. porefill<1.) then - ! eta = 1. - ! eta = 1-porefill - eta = (1-porefill)**2 ! Hudson et al., JGR, 2009 - endif - if (porefill>=1.) eta = 0. - constriction = eta -end function constriction - - - -pure subroutine icechanges_poreonly(nz,z,typeF,typeG,avdrhoP,ypp,B,porefill) - use allinterfaces, except_this_one => icechanges_poreonly - implicit none - integer, intent(IN) :: nz, typeF, typeG - real(8), intent(IN) :: z(nz), ypp(nz), avdrhoP, B - real(8), intent(INOUT) :: porefill(nz) - integer j, erase, newtypeP, ub - real(8) integ - - !----retreat - ! avdrhoP>0 is outward loss from zdepthP - ! avdrhoP<0 means gain at zdepthP or no ice anywhere - if (avdrhoP>0.) then - erase=0 - do j=1,nz - if (typeF>0 .and. j>=typeF) exit ! don't retreat beyond typeF - integ = colint(porefill(1:nz)*z(1:nz),z(1:nz),nz,1,j) - erase = j - if (integ>B*avdrhoP*18./8314.) exit - end do - if (erase>0) porefill(1:erase)=0. - endif - - ! new depth - newtypeP = -9 - do j=1,nz - if (porefill(j)>0.) then - newtypeP = j ! first point with ice - exit - endif - enddo - - !----diffusive filling - ub = typeF - if (newtypeP>0 .and. typeF>0 .and. newtypeP0) then - do j=ub,nz - ! B=Diff/(porosity*icedensity)*86400*365.24*bigstep - porefill(j) = porefill(j) + B*ypp(j) - if (porefill(j)<0.) porefill(j)=0. - if (porefill(j)>1.) porefill(j)=1. - enddo - end if - - !----enact bottom boundary - if (typeG>0) porefill(typeG:nz)=0. - -end subroutine icechanges_poreonly - - - -pure subroutine icechanges(nz,z,typeF,avdrho,avdrhoP,ypp, & - & Diff,porosity,icefrac,bigstep,zdepthT,porefill,typeG) -!*********************************************************** -! advances ice table, advances interface, grows pore ice -! -! a crucial subroutine -!*********************************************************** - use miscparameters, only : icedensity - use allinterfaces, except_this_one => icechanges - implicit none - integer, intent(IN) :: nz, typeF, typeG - real(8), intent(IN) :: z(nz), ypp(nz), avdrho, avdrhoP - real(8), intent(IN) :: Diff, porosity, icefrac, bigstep - real(8), intent(INOUT) :: zdepthT, porefill(nz) - integer j, erase, newtypeP, ub, typeP, typeT - real(8) B, beta, integ - - B = Diff*bigstep*86400.*365.24/(porosity*icedensity) - - ! advance ice table, avdrho>0 is retreat - if (zdepthT>=0. .and. avdrho>0.) then - typeP=-9999; typeT=-9999 - do j=1,nz - if (z(j)>zdepthT) then - typeT=j - exit - endif - enddo - do j=1,nz - if (porefill(j)>0.) then - typeP=j - exit - endif - enddo - if (typeP==typeT) then ! new 2011-09-01 - beta = (1-icefrac)/(1-porosity)/icefrac - beta = Diff*bigstep*beta*86400*365.24/icedensity - zdepthT = sqrt(2*beta*avdrho*18./8314. + zdepthT**2) - endif - endif - if (zdepthT>z(nz)) zdepthT=-9999. - - ! advance interface, avdrhoP>0 is loss from zdepthP - if (avdrhoP>0.) then - erase=0 - do j=1,nz - if (typeF>0 .and. j>=typeF) exit ! don't retreat beyond typeF - if (zdepthT>=0. .and. z(j)>zdepthT) exit - integ = colint(porefill(1:nz)*z(1:nz),z(1:nz),nz,1,j) - erase = j - if (integ>B*avdrhoP*18./8314.) exit - end do - if (erase>0) porefill(1:erase)=0. - endif - - ! new depth - newtypeP = -9 - do j=1,nz - if (zdepthT>=0. .and. z(j)>zdepthT) exit - if (porefill(j)>0.) then - newtypeP = j ! first point with pore ice - exit - endif - enddo - - ! diffusive filling - ub = typeF - if (newtypeP>0 .and. typeF>0 .and. newtypeP0) then - do j=ub,nz - porefill(j) = porefill(j) + B*ypp(j) - if (porefill(j)<0.) porefill(j)=0. - if (porefill(j)>1.) porefill(j)=1. - if (zdepthT>=0. .and. z(j)>zdepthT) exit - enddo - end if - - ! below icetable - if (zdepthT>=0.) then - do j=1,nz - if (z(j)>zdepthT) porefill(j) = icefrac/porosity - enddo - else - ! geothermal lower boundary - if (typeG>0) porefill(typeG:nz)=0. - end if -end subroutine icechanges - - -subroutine assignthermalproperties(nz,thIn,rhoc, & - & ti,rhocv,typeT,icefrac,porosity,porefill) -!********************************************************* -! assign thermal properties of soil -!********************************************************* - implicit none - integer, intent(IN) :: nz - integer, intent(IN), optional :: typeT - real(8), intent(IN), optional :: icefrac - real(8), intent(IN) :: thIn, rhoc - real(8), intent(IN), optional :: porosity, porefill(nz) - real(8), intent(OUT) :: ti(nz), rhocv(nz) - integer j - real(8) newrhoc, newti, fill - real(8), parameter :: NULL=0. - - ti(1:nz) = thIn - rhocv(1:nz) = rhoc - if (typeT>0) then - call soilthprop(porosity,NULL,rhoc,thIn,2,newrhoc,newti,icefrac) - rhocv(typeT:nz) = newrhoc - ti(typeT:nz) = newti - endif - do j=1,nz - fill = porefill(j) ! off by half point - if (fill>0. .and. (typeT<0 .or. (typeT>0 .and. j0.d0) then - solarAttenuation = (1.- fracIR - fracScat)**(1./max(sinbeta,0.04)) - Qdir = Qo*sinbeta*solarAttenuation - Qscat = 0.5*Qo*fracScat - Q = (1.-albedo)*(Qdir + Qscat) - else - Q = 0. - endif - - ! atmospheric IR contribution based on Kieffer et al. (1977), JGR 82, 4249 - Fout = 1.*sigSB*150**4 ! matters only in polar region - Qlw = fracIR*max(Qo*sinbetaNoon,Fout) - - ! net flux = short-wavelength + long-wavelength irradiance - flux_mars77 = Q + Qlw -end function flux_mars77 - - - -pure subroutine flux_mars2(R,decl,latitude,HA,fracIR,fracScat, & - & SlopeAngle,azFac,emax,Qdir,Qscat,Qlw) -!***************************************************************************** -! flux_mars2: Insolation (incoming solar radiation) at Mars on sloped surface; -! returns several irradiances -! -! INPUTS: -! R: distance from sun [AU] -! decl: planetocentric solar declination [radians] -! latitude: [radians] -! HA: hour angle (radians from noon, clockwise) -! fracIR: fraction of absorption -! fracScat: fraction of scattering -! SlopeAngle: >0, [radians] -! azFac: azimuth of topographic gradient (radians east of north) -! azFac=0 is south-facing -! emax: maximum horizon elevation in direction of azimuth [radians] -! OUTPUTS: -! Qdir: direct incoming short-wavelength irradiance [W/m^2] -! Qscat: diffuse short-wavelength irradiance from atmosphere [W/m^2] -! Qlw: diffuse long-wavelength irradiance from atmosphere [W/m^2] -!***************************************************************************** - implicit none - real(8), parameter :: pi=3.1415926535897932, So=1365. - real(8), parameter :: sigSB=5.6704d-8 - real(8), intent(IN) :: R,decl,latitude,HA,SlopeAngle,azFac,emax - real(8), intent(IN) :: fracIR,fracScat - real(8), intent(OUT) :: Qdir,Qscat,Qlw - real(8) c1,s1,solarAttenuation,Qo - real(8) sinbeta,sinbetaNoon,cosbeta,sintheta,azSun,buf,Fout - - c1 = cos(latitude)*cos(decl) - s1 = sin(latitude)*sin(decl) - ! beta = 90 minus incidence angle for horizontal surface - ! beta = elevation of sun above (horizontal) horizon - sinbeta = c1*cos(HA) + s1 - sinbetaNoon = c1 + s1 - sinbetaNoon = max(sinbetaNoon,0.d0) ! horizon - - cosbeta = sqrt(1-sinbeta**2) - buf = (sin(decl)-sin(latitude)*sinbeta)/(cos(latitude)*cosbeta) - if (buf>+1.) buf=+1.d0 ! roundoff - if (buf<-1.) buf=-1.d0 ! roundoff - azSun = acos(buf) - if (sin(HA)>=0) azSun = 2*pi-azSun - -! theta = 90 minus incidence angle for sloped surface - sintheta = cos(SlopeAngle)*sinbeta - & - & sin(SlopeAngle)*cosbeta*cos(azSun-azFac) - if (cosbeta==0) sintheta=cos(SlopeAngle)*sinbeta ! zenith, where azSun=NaN - - if (sintheta<0.) sintheta=0. ! local horizon (self-shadowing) - if (sinbeta<0.) sintheta=0. ! horizontal horizon at infinity - if (sinbeta0.d0) then - Qscat = 0.5*fracScat*Qo - else - Qscat = 0. - endif - -! For a horizontal and unobstructed surface -! absorbed flux = (1-albedo)*(Qdir+Qscat) + emiss*Qlw -! -! For a tilted surface -! absorbed flux = (1-albedo)*(Qdir+Qscat*viewfactor) + emiss*Qlw*viewfactor -! then add irradiance from land in field of view -! in the case of a planar slope, viewfactor = cos(SlopeAngle/2.)**2 -end subroutine flux_mars2 Index: trunk/LMDZ.COMMON/libf/evolution/generalorbit.f =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/generalorbit.f (revision 3492) +++ (revision ) @@ -1,54 +1,0 @@ -C============================================================= -C Subroutine to return distance, longitude, and declination of -C the sun in planetocentric coordinates from orbital elements -C -C INPUTS: -C edays = time since perihelion (earth days) -C a = semimajor axis (AU) -C ecc = eccentricity -C omega = Ls of perihelion, relative to equinox (radians) -C eps = obliquity (radians) -C -C OUTPUTS: -C Ls = areocentric longitude (radians) -C dec = planetocentric solar declination (radians) -C r = heliocentric distance (AU) -C============================================================= - - SUBROUTINE generalorbit(edays,a,ecc,omega,eps,Ls,dec,r) - implicit none - real*8 edays,a,ecc,omega,eps ! input - real*8 Ls,dec,r ! output - real*8 pi,d2r - parameter (pi=3.1415926535897932,d2r=pi/180.d0) - integer j - real*8 M,E,nu,T,Eold - -c T = orbital period (days) - T = sqrt(4*pi**2/(6.674e-11*1.989e30)*(a*149.598e9)**3)/86400. - -c M = mean anomaly (radians) - M = 2.*pi*edays/T ! M=0 at perihelion - -c E = eccentric anomaly -c solve M = E - ecc*sin(E) by newton method - E = M - do j=1,10 - Eold = E - E = E - (E - ecc*sin(E) - M)/(1.-ecc*cos(E)) - if (abs(E-Eold)<1.e-8) exit - enddo - -c nu = true anomaly - !nu = acos(cos(E)-ecc/(1.-ecc*cos(E))) - !nu = sqrt(1-ecc^2)*sin(E)/(1.-ecc*cos(E)) - !nu = atan(sqrt(1-ecc^2)*sin(E)/(1-cos(E))) - nu = 2.*atan(sqrt((1.+ecc)/(1.-ecc))*tan(E/2.)) - - !r = a*(1.-ecc**2)/(1.+ecc*cos(nu)) - r = a*(1-ecc*cos(E)) - Ls = mod(nu + omega,2.*pi) - dec = asin(sin(eps)*sin(Ls)) - - END - Index: trunk/LMDZ.COMMON/libf/evolution/grids.f90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/grids.f90 (revision 3492) +++ (revision ) @@ -1,96 +1,0 @@ -subroutine setgrid(nz,z,zmax,zfac) - ! construct regularly or geometrically spaced 1D grid - ! z(n)-z(1) = 2*z(1)*(zfac**(n-1)-1)/(zfac-1) - ! choice of z(1) and z(2) is compatible with conductionQ - implicit none - integer, intent(IN) :: nz - real(8), intent(IN) :: zfac, zmax - real(8), intent(OUT) :: z(nz) - integer i - real(8) dz - - dz = zmax/nz - do i=1,nz - z(i) = (i-0.5)*dz - enddo - - if (zfac>1.) then - dz = zmax/(3.+2.*zfac*(zfac**(nz-2)-1.)/(zfac-1.)) - z(1) = dz - z(2) = 3*z(1) - do i=3,nz - z(i) = (1.+zfac)*z(i-1) - zfac*z(i-2) - ! z(i) = z(i-1) + zfac*(z(i-1)-z(i-2)) ! equivalent - enddo - endif - -end subroutine setgrid - - - -real(8) function smartzfac(nz_max,zmax,nz_delta,delta) - ! output can be used as input to setgrid - ! produces zfac with desired number of points within depth delta - implicit none - integer, intent(IN) :: nz_max, nz_delta - real(8), intent(IN) :: zmax, delta - integer j - real(8) f, g, gprime - - if (nz_max2.) then - print *,'zfac=',smartzfac - stop 'unwanted result in smartzfac' - endif -end function smartzfac - - - -!-grid-dependent utility functions - -function colint(y,z,nz,i1,i2) - ! column integrates y - implicit none - integer, intent(IN) :: nz, i1, i2 - real(8), intent(IN) :: y(nz), z(nz) - real(8) colint - integer i - real(8) dz(nz) - - dz(1) = (z(2)-0.)/2 - do i=2,nz-1 - dz(i) = (z(i+1)-z(i-1))/2. - enddo - dz(nz) = z(nz)-z(nz-1) - colint = sum(y(i1:i2)*dz(i1:i2)) -end function colint - - - -subroutine heatflux_from_temperature(nz,z,T,k,H) - ! calculates heat flux from temperature profile - ! like k, the heat flux H is defined mid-point - implicit none - integer, intent(IN) :: nz - real(8), intent(IN) :: z(nz), T(nz), k(nz) - real(8), intent(OUT) :: H(nz) - integer j - - ! H(1) = -k(1)*(T(1)-Tsurf)/z(1) - H(1) = 0. ! to avoid ill-defined value - do j=2,nz - H(j) = -k(j)*(T(j)-T(j-1))/(z(j)-z(j-1)) - enddo -end subroutine heatflux_from_temperature Index: trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90 (revision 3492) +++ trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90 (revision 3493) @@ -10,13 +10,14 @@ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -logical, save :: icetable_equilibrium ! Boolean to say if the PEM needs to recompute the icetable depth when at equilibrium -logical, save :: icetable_dynamic ! Boolean to say if the PEM needs to recompute the icetable depth with the dynamic method -real, allocatable, dimension(:,:), save :: porefillingice_depth ! ngrid x nslope: Depth of the ice table [m] -real, allocatable, dimension(:,:), save :: porefillingice_thickness ! ngrid x nslope: Thickness of the ice table [m] +logical, save :: icetable_equilibrium ! Boolean to say if the PEM needs to recompute the icetable depth when at equilibrium +logical, save :: icetable_dynamic ! Boolean to say if the PEM needs to recompute the icetable depth with the dynamic method +real, allocatable, dimension(:,:) :: icetable_depth ! ngrid x nslope: Depth of the ice table [m] +real, allocatable, dimension(:,:) :: icetable_thickness ! ngrid x nslope: Thickness of the ice table [m] +real, allocatable, dimension(:,:,:) :: ice_porefilling ! the amout of porefilling in each layer in each grid [m^3/m^3] !----------------------------------------------------------------------- contains !----------------------------------------------------------------------- -SUBROUTINE ini_ice_table_porefilling(ngrid,nslope) +SUBROUTINE ini_ice_table(ngrid,nslope,nsoil) implicit none @@ -24,19 +25,25 @@ integer, intent(in) :: ngrid ! number of atmospheric columns integer, intent(in) :: nslope ! number of slope within a mesh - -allocate(porefillingice_depth(ngrid,nslope)) -allocate(porefillingice_thickness(ngrid,nslope)) - -END SUBROUTINE ini_ice_table_porefilling - -!----------------------------------------------------------------------- -SUBROUTINE end_ice_table_porefilling - -implicit none - -if (allocated(porefillingice_depth)) deallocate(porefillingice_depth) -if (allocated(porefillingice_thickness)) deallocate(porefillingice_thickness) - -END SUBROUTINE end_ice_table_porefilling +integer, intent(in) :: nsoil ! number of soil layers + +allocate(icetable_depth(ngrid,nslope)) +if (icetable_equilibrium) then + allocate(icetable_thickness(ngrid,nslope)) +else if (icetable_dynamic) then + allocate(ice_porefilling(ngrid,nsoil,nslope)) +endif + +END SUBROUTINE ini_ice_table + +!----------------------------------------------------------------------- +SUBROUTINE end_ice_table + +implicit none + +if (allocated(icetable_depth)) deallocate(icetable_depth) +if (allocated(icetable_thickness)) deallocate(icetable_thickness) +if (allocated(ice_porefilling)) deallocate(ice_porefilling) + +END SUBROUTINE end_ice_table !------------------------------------------------------------------------------------------------------ Index: trunk/LMDZ.COMMON/libf/evolution/misc_params.f90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/misc_params.f90 (revision 3492) +++ (revision ) @@ -1,14 +1,0 @@ -module miscparameters - ! miscellaneous parameters that are very constant - real(8), parameter :: pi=3.1415926535897932, d2r=pi/180. - integer, parameter :: NMAX=1000 - real(8), parameter :: marsDay=88775.244, solsperyear=668.60 - real(8), parameter :: icedensity=927. - real(8), parameter :: earthDay=86400. - real(8), parameter :: sigSB=5.6704e-8 - ! for reference here are some parameters that are hard coded - ! mass of H2O molecule = 18 - ! universal gas constant = 8314 - ! length of Earth year in days = 365.24 -end module miscparameters - Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3492) +++ trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3493) @@ -54,6 +54,6 @@ use orbit_param_criterion_mod, only: orbit_param_criterion use recomp_orb_param_mod, only: recomp_orb_param -use ice_table_mod, only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, & - ini_ice_table_porefilling, icetable_equilibrium, icetable_dynamic, computeice_table_equilibrium,compute_massh2o_exchange_ssi +use ice_table_mod, only: icetable_depth, icetable_thickness, end_ice_table, ice_porefilling, & + ini_ice_table, icetable_equilibrium, icetable_dynamic, computeice_table_equilibrium, compute_massh2o_exchange_ssi use soil_thermalproperties_mod, only: update_soil_thermalproperties use time_phylmdz_mod, only: daysec, dtphys @@ -204,30 +204,28 @@ ! Variables for surface and soil -real, dimension(:,:), allocatable :: tsurf_avg ! Physic x SLOPE field: Averaged Surface Temperature [K] -real, dimension(:,:,:), allocatable :: tsoil_avg ! Physic x SOIL x SLOPE field: Averaged Soil Temperature [K] -real, dimension(:,:,:), allocatable :: tsurf_PCM_timeseries ! ngrid x SLOPE x TIMES field: Surface Temperature in timeseries [K] -real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries ! IG x SOIL x SLOPE x TIMES field: Non averaged Soil Temperature [K] -real, dimension(:,:,:,:), allocatable :: tsoil_anom ! IG x SOIL x SLOPE x TIMES field: Amplitude between instataneous and yearly average soil temperature at the beginning [K] -real, dimension(:,:), allocatable :: tsurf_avg_yr1 ! Physic x SLOPE field: Averaged Surface Temperature of first call of the PCM [K] -real, dimension(:,:), allocatable :: Tsoil_locslope ! Physic x Soil: Intermediate when computing Tsoil [K] -real, dimension(:), allocatable :: Tsurf_locslope ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K] -real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3] -real, dimension(:,:), allocatable :: watersurf_density_avg ! Physic x Slope, water surface density, yearly averaged [kg/m^3] -real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3] -real, dimension(:,:,:), allocatable :: watersoil_density_PEM_avg ! Physic x Soil x Slopes, water soil density, yearly averaged [kg/m^3] -real, dimension(:,:), allocatable :: Tsurfavg_before_saved ! Surface temperature saved from previous time step [K] -real, dimension(:), allocatable :: delta_co2_adsorbded ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2] -real, dimension(:), allocatable :: delta_h2o_adsorbded ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2] -real :: totmassco2_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg] -real :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg] -logical :: bool_sublim ! logical to check if there is sublimation or not -logical, dimension(:,:), allocatable :: co2ice_disappeared ! logical to check if a co2 ice reservoir already disappeared at a previous timestep -real, dimension(:,:), allocatable :: porefillingice_thickness_prev_iter ! ngrid x nslope: Thickness of the ice table at the previous iteration [m] -real, dimension(:), allocatable :: delta_h2o_icetablesublim ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg] - -! Dynamic ice table -real, dimension(:,:), allocatable :: ss_ice_pf ! the amout of porefilling in each layer in each grid [m^3/m^3] -real, dimension(:,:), allocatable :: ss_ice_pf_new ! the amout of porefilling in each layer in each grid after the ss module[m^3/m^3] -real, dimension(:,:), allocatable :: porefillingice_depth_new ! new pure ice table depth +real, dimension(:,:), allocatable :: tsurf_avg ! Physic x SLOPE field: Averaged Surface Temperature [K] +real, dimension(:,:,:), allocatable :: tsoil_avg ! Physic x SOIL x SLOPE field: Averaged Soil Temperature [K] +real, dimension(:,:,:), allocatable :: tsurf_PCM_timeseries ! ngrid x SLOPE x TIMES field: Surface Temperature in timeseries [K] +real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries ! IG x SOIL x SLOPE x TIMES field: Non averaged Soil Temperature [K] +real, dimension(:,:,:,:), allocatable :: tsoil_anom ! IG x SOIL x SLOPE x TIMES field: Amplitude between instataneous and yearly average soil temperature at the beginning [K] +real, dimension(:,:), allocatable :: tsurf_avg_yr1 ! Physic x SLOPE field: Averaged Surface Temperature of first call of the PCM [K] +real, dimension(:,:), allocatable :: Tsoil_locslope ! Physic x Soil: Intermediate when computing Tsoil [K] +real, dimension(:), allocatable :: Tsurf_locslope ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K] +real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3] +real, dimension(:,:), allocatable :: watersurf_density_avg ! Physic x Slope, water surface density, yearly averaged [kg/m^3] +real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3] +real, dimension(:,:,:), allocatable :: watersoil_density_PEM_avg ! Physic x Soil x Slopes, water soil density, yearly averaged [kg/m^3] +real, dimension(:,:), allocatable :: Tsurfavg_before_saved ! Surface temperature saved from previous time step [K] +real, dimension(:), allocatable :: delta_co2_adsorbded ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2] +real, dimension(:), allocatable :: delta_h2o_adsorbded ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2] +real :: totmassco2_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg] +real :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg] +logical :: bool_sublim ! logical to check if there is sublimation or not +logical, dimension(:,:), allocatable :: co2ice_disappeared ! logical to check if a co2 ice reservoir already disappeared at a previous timestep +real, dimension(:,:), allocatable :: icetable_thickness_prev_iter ! ngrid x nslope: Thickness of the ice table at the previous iteration [m] +real, dimension(:), allocatable :: delta_h2o_icetablesublim ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg] +real, dimension(:), allocatable :: porefill ! Pore filling (output) to compute the dynamic ice table +real :: ssi_depth ! Ice table depth (output) to compute the dynamic ice table + ! Some variables for the PEM run real, parameter :: year_step = 1 ! Timestep for the pem @@ -557,5 +555,5 @@ allocate(delta_co2_adsorbded(ngrid)) allocate(co2ice_disappeared(ngrid,nslope)) -allocate(porefillingice_thickness_prev_iter(ngrid,nslope)) +allocate(icetable_thickness_prev_iter(ngrid,nslope)) allocate(delta_h2o_icetablesublim(ngrid)) allocate(delta_h2o_adsorbded(ngrid)) @@ -583,6 +581,6 @@ call end_adsorption_h_PEM call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM) -call end_ice_table_porefilling -call ini_ice_table_porefilling(ngrid,nslope) +call end_ice_table +call ini_ice_table(ngrid,nslope,nsoilmx_PEM) if (soil_pem) then @@ -646,9 +644,9 @@ endif -call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, & - porefillingice_thickness,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,ps_timeseries, & - tsoil_phys_PEM_timeseries,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice,global_avg_press_PCM, & - watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbded_phys,delta_co2_adsorbded, & - h2o_adsorbded_phys,delta_h2o_adsorbded,stratif) +call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth, & + icetable_thickness,ice_porefilling,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys, & + ps_timeseries,tsoil_phys_PEM_timeseries,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice, & + global_avg_press_PCM,watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbded_phys, & + delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,stratif) ! We save the places where h2o ice is sublimating @@ -928,18 +926,24 @@ ! II_d.3 Update the ice table if (icetable_equilibrium) then - write(*,*) "Compute ice table" - porefillingice_thickness_prev_iter = porefillingice_thickness - call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness) - call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere + write(*,*) "Compute ice table (equilibrium method)" + icetable_thickness_prev_iter = icetable_thickness + call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness) + call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_prev_iter,icetable_thickness,icetable_depth,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere else if (icetable_dynamic) then - write(*,*) "Compute ice table" - porefillingice_thickness_prev_iter = porefillingice_thickness - call dyn_ss_ice_m_wrapper(ngrid,nsoilmx,TI_PEM,ps,mmol(igcm_h2o_vap),tsoil_PEM,porefillingice_depth,ss_ice_pf,ss_ice_pf_new,porefillingice_depth_new) - - call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere + write(*,*) "Compute ice table (dynamic method)" + allocate(porefill(nsoilmx_PEM)) + do ig = 1,ngrid + do islope = 1,nslope + call dyn_ss_ice_m(icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM(ig,:,islope),nsoilmx_PEM,TI_PEM(ig,1,nslope),ps(ig),sum(q_h2o_PEM_phys(ig,:))/size(q_h2o_PEM_phys,2),ice_porefilling(ig,:,islope),porefill,ssi_depth) + icetable_depth(ig,islope) = ssi_depth + ice_porefilling(ig,:,islope) = porefill + enddo + enddo + deallocate(porefill) + call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_prev_iter,icetable_thickness,icetable_depth,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere endif ! II_d.4 Update the soil thermal properties - call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tend_h2o_ice,h2o_ice,global_avg_press_new,porefillingice_depth,porefillingice_thickness,TI_PEM) + call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tend_h2o_ice,h2o_ice,global_avg_press_new,icetable_depth,icetable_thickness,TI_PEM) ! II_d.5 Update the mass of the regolith adsorbed @@ -952,7 +956,7 @@ totmassh2o_adsorbded = 0. do ig = 1,ngrid - do islope =1, nslope + do islope = 1,nslope do l = 1,nsoilmx_PEM - if (l==1) then + if (l == 1) then totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) @@ -960,7 +964,7 @@ subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) else - totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* & + totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) - totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* & + totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) endif @@ -987,9 +991,9 @@ call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope)) if (icetable_equilibrium) then - call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,porefillingice_depth(:,islope)) - call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,porefillingice_thickness(:,islope)) + call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope)) + call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,icetable_thickness(:,islope)) else if (icetable_dynamic) then - call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,porefillingice_depth(:,islope)) - call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,porefillingice_thickness(:,islope)) + call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope)) + call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,icetable_thickness(:,islope)) endif @@ -1205,5 +1209,5 @@ call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist) call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, & - TI_PEM,porefillingice_depth,porefillingice_thickness, & + TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, & co2_adsorbded_phys,h2o_adsorbded_phys,h2o_ice,stratif) write(*,*) "restartpem.nc has been written" @@ -1227,5 +1231,5 @@ deallocate(watersurf_density_avg,watersoil_density_timeseries,Tsurfavg_before_saved) deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_avg) -deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter) +deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim,icetable_thickness_prev_iter) deallocate(is_co2ice_ini,co2ice_disappeared,ini_co2ice_sublim,ini_h2oice_sublim,stratif) !----------------------------- END OUTPUT ------------------------------ Index: trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 (revision 3492) +++ trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 (revision 3493) @@ -7,7 +7,7 @@ !======================================================================= -SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table,ice_table_thickness, & - tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice, & - global_avg_pressure,watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys, & +SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table_depth,ice_table_thickness, & + ice_porefilling,tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice, & + global_avg_pressure,watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys, & m_h2o_regolith_phys,deltam_h2o_regolith_phys,stratif) @@ -60,6 +60,7 @@ real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM ! soil (mid-layer) thermal inertia in the PEM grid [SI] real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: tsoil_PEM ! soil (mid-layer) temperature [K] -real, dimension(ngrid,nslope), intent(inout) :: ice_table ! Ice table depth [m] +real, dimension(ngrid,nslope), intent(inout) :: ice_table_depth ! Ice table depth [m] real, dimension(ngrid,nslope), intent(inout) :: ice_table_thickness ! Ice table thickness [m] +real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: ice_porefilling ! Subsurface ice pore filling [m3/m3] real, dimension(ngrid,nsoil_PEM,nslope,timelen), intent(inout) :: tsoil_inst ! instantaneous soil (mid-layer) temperature [K] real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2] @@ -107,5 +108,5 @@ !0.2 Set to default values -ice_table = -1. ! by default, no ice table +ice_table_depth = -1. ! by default, no ice table ice_table_thickness = -1. !1. Run @@ -225,6 +226,6 @@ do ig = 1,ngrid if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then - inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ & - (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ & + inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ & + (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + & ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2)))) @@ -308,14 +309,31 @@ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !3. Ice Table - if (icetable_equilibrium .or. icetable_dynamic) then - call get_field("ice_table",ice_table,found) + if (icetable_equilibrium) then + call get_field("ice_table_depth",ice_table_depth,found) if (.not. found) then - write(*,*)'PEM settings: failed loading ' + write(*,*)'PEM settings: failed loading ' write(*,*)'will reconstruct the values of the ice table given the current state' - call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table,ice_table_thickness) - call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM) + call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table_depth,ice_table_thickness) + call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM) do islope = 1,nslope call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) enddo + endif + write(*,*) 'PEMETAT0: ICE TABLE done' + else if (icetable_dynamic) then + call get_field("ice_table_depth",ice_table_depth,found) + if (.not. found) then + write(*,*)'PEM settings: failed loading ' + write(*,*)'will reconstruct the values of the ice table given the current state' + call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table_depth,ice_table_thickness) + call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM) + do islope = 1,nslope + call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) + enddo + endif + call get_field("ice_porefilling",ice_porefilling,found) + if (.not. found) then + write(*,*)'PEM settings: failed loading ' + ice_porefilling = 0. endif write(*,*) 'PEMETAT0: ICE TABLE done' @@ -354,5 +372,5 @@ call close_startphy -else !No startfi, let's build all by hand +else !No startpem, let's build all by hand write(*,*)'There is no "startpem.nc"!' @@ -484,7 +502,15 @@ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !c) Ice table - if (icetable_equilibrium .or. icetable_dynamic) then - call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table,ice_table_thickness) - call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM) + if (icetable_equilibrium) then + call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table_depth,ice_table_thickness) + call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM) + do islope = 1,nslope + call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) + enddo + write(*,*) 'PEMETAT0: Ice table done' + else if (icetable_dynamic) then + ice_porefilling = 0. + ice_table_depth = 0. + call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table_depth,ice_table_thickness,TI_PEM) do islope = 1,nslope call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope)) Index: trunk/LMDZ.COMMON/libf/evolution/pemredem.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pemredem.F90 (revision 3492) +++ trunk/LMDZ.COMMON/libf/evolution/pemredem.F90 (revision 3493) @@ -58,5 +58,5 @@ SUBROUTINE pemdem1(filename,i_myear,nsoil_PEM,ngrid,nslope,tsoil_slope_PEM,inertiesoil_slope_PEM, & - ice_table_depth,ice_table_thickness,m_co2_regolith,m_h2o_regolith,h2o_ice,stratif) + ice_table_depth,ice_table_thickness,ice_porefilling,m_co2_regolith,m_h2o_regolith,h2o_ice,stratif) ! write time-dependent variable to restart file @@ -78,4 +78,5 @@ real, dimension(ngrid,nslope), intent(in) :: ice_table_depth ! under mesh bining according to slope real, dimension(ngrid,nslope), intent(in) :: ice_table_thickness ! under mesh bining according to slope +real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: ice_porefilling ! under mesh bining according to slope real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: m_co2_regolith, m_h2o_regolith real, dimension(ngrid,nslope), intent(in) :: h2o_ice @@ -122,4 +123,5 @@ call put_field("ice_table_depth","Depth of ice table",ice_table_depth,Year) call put_field("ice_table_thickness","Depth of ice table",ice_table_thickness,Year) + call put_field("ice_porefilling","Subsurface ice pore filling",ice_porefilling,Year) call put_field("inertiedat_PEM","Thermal inertie of PEM ",inertiedat_PEM,Year) endif ! soil_pem Index: trunk/LMDZ.COMMON/libf/evolution/psv.f =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/psv.f (revision 3492) +++ (revision ) @@ -1,54 +1,0 @@ - real*8 function psv(T) -C saturation vapor pressure of H2O ice [Pascal] -C input is temperature [Kelvin] - implicit none - real*8 T - -C-----parametrization 1 -c real*8 DHmelt,DHvap,DHsub,R,pt,Tt,C -c parameter (DHmelt=6008.,DHvap=45050.) -c parameter (DHsub=DHmelt+DHvap) ! sublimation enthalpy [J/mol] -c parameter (R=8.314,pt=6.11e2,Tt=273.16) -c C = (DHsub/R)*(1./T - 1./Tt) -c psv = pt*exp(-C) - -C-----parametrization 2 -C eq. (2) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005) -C differs from parametrization 1 by only 0.1% - real*8 A,B - parameter (A=-6143.7, B=28.9074) - psv = exp(A/T+B) ! Clapeyron - -C-----parametrization 3 -C eq. (7) in Murphy & Koop, Q. J. R. Meteor. Soc. 131, 1539 (2005) -c psv = exp(9.550426 - 5723.265/T + 3.53068*log(T) - 0.00728332*T) - - end - - - - real*8 function frostpoint(p) -C inverse of psv -C input is partial pressure [Pascal] -C output is temperature [Kelvin] - implicit none - real*8 p - -C-----inverse of parametrization 1 -c real*8 DHmelt,DHvap,DHsub,R,pt,Tt -c parameter (DHmelt=6008.,DHvap=45050.) -c parameter (DHsub=DHmelt+DHvap) -c parameter (R=8.314,pt=6.11e2,Tt=273.16) -c frostpoint = 1./(1./Tt-R/DHsub*log(p/pt)) - -C-----inverse of parametrization 2 -C inverse of eq. (2) in Murphy & Koop (2005) - real*8 A,B - parameter (A=-6143.7, B=28.9074) - frostpoint = A / (log(p) - B) - -C-----approximate inverse of parametrization 3 -C eq. (8) in Murphy & Koop (2005) -c frostpoint = (1.814625*log(p) + 6190.134)/(29.120 - log(p)) - - end Index: trunk/LMDZ.COMMON/libf/evolution/psvco2.f =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/psvco2.f (revision 3492) +++ (revision ) @@ -1,50 +1,0 @@ - real*8 function psvco2(T) -C solid-vapor transition for CO2 -C returns saturation pressure in Pascal, input is temperature in Kelvin - implicit none - real*8 T - - psvco2 = exp(23.3494 - 3182.48/T)*100. - end - - - real*8 function tfrostco2(p) -C the inverse of function psvco2 -C input is pressure in Pascal, output is temperature in Kelvin - implicit none - real*8 p - tfrostco2 = 3182.48/(23.3494+log(100./p)) - end - - - - -C------------------------------------------------------------------------------ -C Antoine equation parameters from NIST, 154K-196K -C based on Giauque and Egan (1937) -C A=6.81228, B=1301.679, C=-3.494 -c p = 10**(A-(B/(T+C)))*1.e5 - -C Expressions from Int. Crit. Tabl. 3.207, based on many references -C mm2Pa = 133.32 -C -135C to -56.7C (138K-216K) -c p = 10**(-0.05223*26179.3/T + 9.9082)*mm2Pa -C -183C to -135C (90K-138K) -c p = 10**(-1275.62/T + 0.006833*T + 8.3071)*mm2Pa -C Expressions from Int. Crit. Tabl. 3.208, based on Henning -C -110 to -80C (163K-193K) -c p = 10**(- 1279.11/T + 1.75*log10(T) - 0.0020757*T + 5.85242)*mm2Pa -c p = 10**(- 1352.6/T + 9.8318)*mm2Pa - -C Mars book (1992), p959, fit by chapter authors -c p = exp(23.3494 - 3182.48/T)*100. ! 120-160K -c p = exp(25.2194 - 3311.57/T - 6.71e-3*T)*100 ! 100-170K -C Mars book (1992), p960, based on Miller & Smythe (1970) -c p = exp(26.1228 - 3385.26/T - 9.4461e-3*T)*100 ! 123-173K - -C Fray & Schmitt, PSS 57, 2053 (2009) -C A0=1.476e1, A1=-2.571e3, A2=-7.781e4, A3=4.325e6, A4=-1.207e8, A5=1.350e9 -c p = exp(A0+A1/T+A2/T**2+A3/T**3+A4/T**4+A5/T**5)*1e5 ! 40K-194.7K -c A0=1.861e1, A1=-4.154e3, A2=1.041e5 -c p = exp(A0 + A1/T + A2/T**2)*1e5 ! 194.7K-216.58K -C------------------------------------------------------------------------------ Index: trunk/LMDZ.COMMON/libf/evolution/soilthprop_mars.f90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/soilthprop_mars.f90 (revision 3492) +++ (revision ) @@ -1,152 +1,0 @@ -subroutine soilthprop(porosity,fill,rhocobs,tiobs,layertype, & - & newrhoc,newti,icefrac) -!*********************************************************************** -! soilthprop: assign thermal properties of icy soil or dirty ice -! -! porositiy = void space / total volume -! rhof = density of free ice in space not occupied by regolith [kg/m^3] -! fill = rhof/icedensity <=1 (only relevant for layertype 1) -! rhocobs = heat capacity per volume of dry regolith [J/m^3] -! tiobs = thermal inertia of dry regolith [SI-units] -! layertype: 1=interstitial ice, 2=pure ice or ice with dirt -! 3=pure ice, 4=ice-cemented soil, 5=custom values -! icefrac: fraction of ice in icelayer (only relevant for layertype 2) -! output are newti and newrhoc -!*********************************************************************** - implicit none - integer, intent(IN) :: layertype - real(8), intent(IN) :: porosity, fill, rhocobs, tiobs - real(8), intent(OUT) :: newti, newrhoc - real(8), intent(IN) :: icefrac - real(8) kobs, cice, icedensity, kice - !parameter (cice=2000.d0, icedensity=926.d0, kice=2.4d0) ! unaffected by scaling - parameter (cice=1540.d0, icedensity=927.d0, kice=3.2d0) ! at 198 Kelvin - real(8) fA, ki0, ki, k - real(8), parameter :: kw=3. ! Mellon et al., JGR 102, 19357 (1997) - - kobs = tiobs**2/rhocobs - ! k, rhoc, and ti are defined in between grid points - ! rhof and T are defined on grid points - - newrhoc = -9999. - newti = -9999. - - select case (layertype) - case (1) ! interstitial ice - newrhoc = rhocobs + porosity*fill*icedensity*cice - if (fill>0.) then - !--linear addition (option A) - k = porosity*fill*kice + kobs - !--Mellon et al. 1997 (option B) - ki0 = porosity/(1/kobs-(1-porosity)/kw) - fA = sqrt(fill) - ki = (1-fA)*ki0 + fA*kice - !k = kw*ki/((1-porosity)*ki+porosity*kw) - else - k = kobs - endif - newti = sqrt(newrhoc*k) - - case (2) ! massive ice (pure or dirty ice) - newrhoc = rhocobs*(1.-icefrac)/(1.-porosity) + icefrac*icedensity*cice - k = icefrac*kice + (1.-icefrac)*kw - newti = sqrt(newrhoc*k) - - case (3) ! all ice, special case of layertype 2, which doesn't use porosity - newrhoc = icedensity*cice - k = kice - newti = sqrt(newrhoc*k) - - case (4) ! pores completely filled with ice, special case of layertype 1 - newrhoc = rhocobs + porosity*icedensity*cice - k = porosity*kice + kobs ! option A, end-member case of type 1, option A - !k = kw*kice/((1-porosity)*kice+porosity*kw) ! option B, harmonic average - newti = sqrt(newrhoc*k) - - case (5) ! custom values - ! values from Mellon et al. (2004) for ice-cemented soil - newrhoc = 2018.*1040. - k = 2.5 - newti = sqrt(newrhoc*k) - - case default - error stop 'invalid layer type' - - end select - -end subroutine soilthprop - - - -subroutine smartgrid(nz,z,zdepth,thIn,rhoc,porosity,ti,rhocv,layertype,icefrac) -!*********************************************************************** -! smartgrid: returns intelligently spaced grid and appropriate -! values of thermal inertia ti and rho*c in icy layer -! -! INPUTS: -! nz = number of grid points -! z = grid spacing for dry regolith -! (will be partially overwritten) -! zdepth = depth where ice table starts -! negative values indicate no ice -! rhoc = heat capacity per volume of dry regolith [J/m^3] -! thIn = thermal inertia of dry regolith [SI-units] -! porosity = void space / total volume -! layertypes are explained below -! icefrac = fraction of ice in icelayer -! -! OUTPUTS: z = grid coordinates -! vectors ti and rhocv -!*********************************************************************** - implicit none - integer, intent(IN) :: nz, layertype - real(8), intent(IN) :: zdepth, thIn, rhoc, porosity, icefrac - real(8), intent(INOUT) :: z(nz) - real(8), intent(OUT) :: ti(nz), rhocv(nz) - integer j, b - real(8) stretch, newrhoc, newti - real(8), parameter :: NULL=0. - - if (zdepth>0 .and. zdepth',n -c stop -c endif - bet=b(1) - u(1)=r(1)/bet - do 11 j=2,n - gam(j)=c(j-1)/bet - bet=b(j)-a(j)*gam(j) -c if(bet.eq.0.)pause 'tridag failed' - u(j)=(r(j)-a(j)*u(j-1))/bet -11 continue - do 12 j=n-1,1,-1 - u(j)=u(j)-gam(j+1)*u(j+1) -12 continue - return - END -C (C) Copr. 1986-92 Numerical Recipes Software 0(9p#31(+.