Index: trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90 (revision 3411) +++ trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90 (revision 3412) @@ -203,4 +203,7 @@ call getin_p("convergeps",convergeps) if (is_master) write(*,*) trim(rname)//" convergeps = ",convergeps + conservn2=.false. ! default value + call getin_p("conservn2",conservn2) + if (is_master) write(*,*) trim(rname)//" conservn2 = ",conservn2 if (is_master) write(*,*) trim(rname)//"KBO runs (eris, makemake) ?" Index: trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (revision 3411) +++ trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (revision 3412) @@ -54,8 +54,9 @@ fast,fasthaze,haze,metcloud,monoxcloud,& n2cond,nearn2cond,noseason_day,conservn2, & - convergeps,kbo,triton,paleo,paleoyears, & + convergeps,kbo,triton,paleo,paleoyears,glaflow, & carbox, methane,& oldplutovdifc,oldplutocorrk,oldplutosedim, & - aerohaze,haze_proffix,source_haze,& + aerohaze,haze_proffix,source_haze, tsurfmax, & + albmin_ch4, & season, sedimentation,generic_condensation, & specOLR, & @@ -291,4 +292,7 @@ REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer) + REAL,SAVE :: ptime0 ! store the first time + REAL dstep + REAL,SAVE :: glastep=20 ! step in pluto day to spread glacier @@ -617,4 +621,8 @@ albedo_snow_SPECTV(:)=0.0 albedo_n2_ice_SPECTV(:)=0.0 + + ptime0=ptime + write (*,*) 'In physiq ptime0 =', ptime + call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_n2_ice_SPECTV) @@ -759,4 +767,9 @@ taux(1:ngrid) = 0.0 tauy(1:ngrid) = 0.0 + + if (conservn2) then + write(*,*) 'conservn2 iniloop' + call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1)) + endif zday=pday+ptime ! Compute time, in sols (and fraction thereof). @@ -848,4 +861,10 @@ enddo endif + + if (conservn2) then + write(*,*) 'conservn2 thermo' + call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1)) + endif + !--------------------------------- ! II. Compute radiative tendencies @@ -1094,4 +1113,8 @@ !! call writediagfi(ngrid,"vdifc_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas)) + if (conservn2) then + write(*,*) 'conservn2 radiat' + call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1)) + endif ! -------------------------------------------- @@ -1229,4 +1252,9 @@ endif ! end of 'calldifv' + if (conservn2) then + write(*,*) 'conservn2 diff' + call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1)+ & + dqsurf(:,1)*ptimestep) + endif !------------------- @@ -1322,4 +1350,9 @@ endif ! end of 'n2cond' + if (conservn2) then + write(*,*) 'conservn2 n2cond' + call testconservmass(ngrid,nlayer,pplev(:,1)+ & + pdpsrf(:)*ptimestep,qsurf(:,1)+dqsurf(:,1)*ptimestep) + endif !--------------------------------------------- @@ -1592,6 +1625,6 @@ if (conservn2) then write(*,*) 'conservn2 tracer' - ! call testconservmass(ngrid,nlayer,pplev(:,1)+ & - ! pdpsrf(:)*ptimestep,qsurf(:,1)) + call testconservmass(ngrid,nlayer,pplev(:,1)+ & + pdpsrf(:)*ptimestep,qsurf(:,1)) endif @@ -1646,9 +1679,21 @@ - ! ! Increment surface temperature - ! if(ok_slab_ocean)then !AF24: removed - +! VII.1 Increment surface temperature +! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid) - ! Compute soil temperatures and subsurface heat flux. + + ! Prevent surface (.e.g. non volatile ch4) to exceed max temperature + ! Lellouch et al., 2000,2011 + IF (tsurfmax) THEN + DO ig=1,ngrid + if (albedo_equivalent(ig).gt.albmin_ch4.and. & + qsurf(ig,igcm_n2).eq.0.) then + tsurf(ig)=min(tsurf(ig),54.) + endif + ENDDO + ENDIF + +! VII.2 Compute soil temperatures and subsurface heat flux. +! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (callsoil) then call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat, & @@ -1656,6 +1701,9 @@ endif - -! if (ok_slab_ocean) then !AF24: removed + ! ! For output : + ! tidat_out(:,:)=0. + ! DO l=1,min(nlayermx,nsoilmx) + ! tidat_out(:,l)=tidat(:,l) + ! ENDDO ! Test energy conservation @@ -1665,4 +1713,49 @@ endif + + +! VII.3 multiply tendencies of cond/subli for paleo loop only in the +! last Pluto year of the simulation +! Year day must be adapted in the startfi for each object +! Paleo uses year_day to calculate the annual mean tendancies +! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + IF (paleo) then + if (zday.gt.day_ini+ptime0+nday-year_day) then + DO iq=1,nq + DO ig=1,ngrid + qsurfyear(ig,iq)=qsurfyear(ig,iq)+ & + (qsurf(ig,iq)-qsurf1(ig,iq)) !kg m-2 !ptimestep + ENDDO + ENDDO + endif + endif + +! VII.4 Glacial flow at each timestep glastep or at lastcall +! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + IF (fast.and.glaflow) THEN + if ((mod(zday-day_ini-ptime0,glastep)).lt.1. & + .or.lastcall) then + IF (lastcall) then + dstep=mod(zday-day_ini-ptime0,glastep)*daysec + else + dstep=glastep*daysec + endif + zdqflow(:,:)=qsurf(:,:) + IF (paleo) then + call spreadglacier_paleo(ngrid,nq,qsurf, & + phisfinew,dstep,tsurf) + else + call spreadglacier_simple(ngrid,nq,qsurf,dstep) + endif + zdqflow(:,:)=(zdqflow(:,:)-qsurf(:,:))/dstep + + if (conservn2) then + write(*,*) 'conservn2 glaflow' + call testconservmass(ngrid,nlayer,pplev(:,1)+ & + pdpsrf(:)*ptimestep,qsurf(:,1)) + endif + + endif + ENDIF !--------------------------------------------------- @@ -1816,5 +1909,4 @@ - ! ! Test for water conservation. !AF24: removed ! Calculate RH_generic (Generic Relative Humidity) for diagnostic. @@ -2271,9 +2363,9 @@ do iq=1,nq - ! call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq)) + call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq)) ! call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', & ! 'kg m^-2',2,qsurf_hist(1,iq) ) - ! call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', & - ! 'kg m^-2',2,qcol(1,iq) ) + call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col', & + 'kg m^-2',2,qcol(1,iq) ) ! call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf', & ! 'kg m^-2',2,qsurf(1,iq) ) Index: trunk/LMDZ.PLUTO/libf/phypluto/spreadglacier_paleo.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/spreadglacier_paleo.F90 (revision 3412) +++ trunk/LMDZ.PLUTO/libf/phypluto/spreadglacier_paleo.F90 (revision 3412) @@ -0,0 +1,340 @@ + subroutine spreadglacier_paleo(ngrid,nq,pqsurf, & + phisfi0,timstep,ptsurf) + + use callkeys_mod, only: carbox, methane + use comcstfi_mod, only: g + use comgeomfi_h + use geometry_mod, only: latitude, longitude, cell_area + use radinc_h, only : naerkind + use tracer_h, only: igcm_n2,igcm_ch4_ice,igcm_co_ice + + implicit none + +!================================================================== +! Purpose +! ------- +! Spreading the glacier : N2, cf Umurhan et al 2017 +! +! Inputs +! ------ +! ngrid Number of vertical columns +! pqsurf(ngrid,nq) N2 ice tracer on surface kg/m2 +! phisfi0 = phisfinew the geopotential of the bedrock +! below the N2 ice +! ptsurf surface temperature +! timstep dstep = timestep in pluto day for the call of glacial flow + +! Outputs +! ------- +! pqsurf(ngrid,nq) new value for N2 ice tracer on surface kg/m2 +! +! Authors +! ------- +! Tanguy Bertrand (2015,2022) +! +!================================================================== + +#include "dimensions.h" + +!----------------------------------------------------------------------- +! Arguments + + INTEGER ngrid, nq, ig, i + REAL :: pqsurf(ngrid,nq) + REAL :: phisfi0(ngrid) + REAL :: ptsurf(ngrid) + REAL timstep +!----------------------------------------------------------------------- +! Local variables + REAL tgla(ngrid),tbase(ngrid) !temperature at the base of glacier different of ptsurf + REAL :: zdqsurf(ngrid) ! tendency of qsurf + REAL massflow ! function + REAL hlim,qlim,difflim,diff,stock,H0,totmasstmp + INTEGER compt !compteur + INTEGER slid ! option slid 0 or 1 + INTEGER :: edge(4) ! index of the adjecent points, N,S,E,W + INTEGER cas,inddeb,indfin + REAL secu,dqch4,dqco + REAL :: masstmp(ngrid) + + REAL :: masstot + +!---------------- INPUT ------------------------------------------------ + + difflim=0.5 ! limit height (m) between two adjacent point to start spreading the glacier + zdqsurf(:)=0. +!--------------- Dimensions ------------------------------------- + + ! distance between two adjacent latitude + !distlim_lat=pi*rad/jjm/1000. ! km + !distlim_diag=(distlim_lat*distlim_lat+distlim_lon*distlim_lon)**0.5 + + !! Threshold for ice thickness for which we activate the glacial flow + hlim=10. !m + qlim=hlim*1000. ! kg + !! Security for not depleted all ice too fast in one timestep + secu=4 + + !************************************* + ! Loop over all points + !************************************* + DO ig=1,ngrid + !if (ig.eq.ngrid) then + ! print*, 'qpole=',pqsurf(ig,igcm_n2),qlim + !endif + + !************************************* + ! if significant amount of ice, where qsurf > qlim + !************************************* + IF (pqsurf(ig,igcm_n2).gt.qlim) THEN + + !************************************* + ! temp glacier with gradient 15 K / km + !************************************* + ! surface temperature pts + tgla(ig)=ptsurf(ig) + ! temperature at the base (we neglect convection) + tbase(ig)=tgla(ig)+20.*pqsurf(ig,igcm_n2)/1.e6 ! Umurhan et al. + if (tbase(ig).gt.63.) then + slid=1 ! activate slide + else + slid=0 + endif + + !************************************* + ! Selection of the adjacent index depending on the grid point + !************************************* + !! poles + IF (ig.eq.1) then + cas=0 + inddeb=2 ! First adj grid point + indfin=iim+1 ! Last adj grid point + ELSEIF (ig.eq.ngrid) then + cas=10 + inddeb=ngrid-iim + indfin=ngrid-1 + !! 9 other cases: edges East, west, or in the center , at edges north south or in the center + ELSEIF (mod(ig,iim).eq.2) then ! Edge east : N,S,W,E + IF (ig.eq.2) then + cas=1 + edge = (/1, ig+iim,ig-1+iim,ig+1 /) + ELSEIF (ig.eq.ngrid-iim) then + cas=3 + edge = (/ig-iim, ngrid,ig-1+iim,ig+1 /) + ELSE + cas=2 + edge = (/ig-iim, ig+iim,ig-1+iim,ig+1 /) + ENDIF + ELSEIF (mod(ig,iim).eq.1) then ! Edge west + IF (ig.eq.iim+1) then + cas=7 + edge = (/1,ig+iim,ig-1,ig+1-iim /) + ELSEIF (ig.eq.ngrid-1) then + cas=9 + edge = (/ig-iim,ngrid,ig-1,ig+1-iim /) + ELSE + cas=8 + edge = (/ig-iim,ig+iim,ig-1,ig+1-iim /) + ENDIF + ELSE + IF ((ig.gt.2).and.(ig.lt.iim+1)) then + cas=4 + edge = (/1,ig+iim,ig-1,ig+1 /) + ELSEIF ((ig.gt.ngrid-iim).and.(ig.lt.ngrid-1)) then + cas=6 + edge = (/ig-iim,ngrid,ig-1,ig+1 /) + ELSE + cas=5 + edge = (/ig-iim,ig+iim,ig-1,ig+1 /) + ENDIF + ENDIF + + masstmp(:)=0. ! mass to be transferred + totmasstmp=0. ! total mass to be transferred + H0=phisfi0(ig)/g+pqsurf(ig,igcm_n2)/1000. ! height of ice at ig (m) + + !************************************* + !!!!!!!!!!!!!!!!! POLE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !************************************* + IF ((cas.eq.0).or.(cas.eq.10)) then ! Mean over all latitudes near the pole + + !call testconservmass2d(ngrid,pqsurf(:,1)) + !! 1) Compute Diff H0-H1 and mass to be transferred for each adjacent point + DO i=inddeb,indfin + diff=H0-(phisfi0(i)/g+pqsurf(i,igcm_n2)/1000.) ! Height difference between ig and adjacent points (m) + if (diff.gt.difflim) then + masstmp(i)=massflow(ig,i,tgla(ig),diff,pqsurf(ig,igcm_n2),timstep,slid) ! Mass to be transferred (kg/m2/s) + else + masstmp(i)=0. + endif + totmasstmp=totmasstmp+masstmp(i) ! mass totale to be transferred if assume only one adjecent point + ENDDO + if (totmasstmp.gt.0.) then + !! 2) Available mass to be transferred + stock=maxval(masstmp(:))/secu ! kg/m2/s + !! 3) Real amounts of ice to be transferred : + masstmp(:)=masstmp(:)/totmasstmp*stock*cell_area(ig)/cell_area(inddeb) !kg/m2/s all area around the pole are the same + !! 4) Tendencies + zdqsurf(ig)=-stock !kg/m2/s + !! 5) Update PQSURF + + ! move CH4 and CO too + if (methane) then + !! Variation of CH4 ice + dqch4=pqsurf(ig,igcm_ch4_ice)*(zdqsurf(ig)/pqsurf(ig,igcm_n2)*timstep) ! amout of ch4 (kg/m2) to be removed at grid ig (negative) + !! remove CH4 + pqsurf(ig,igcm_ch4_ice)=pqsurf(ig,igcm_ch4_ice)+dqch4 + !! add CH4 + do i=inddeb,indfin + pqsurf(i,igcm_ch4_ice)=pqsurf(i,igcm_ch4_ice)-masstmp(i)/stock*dqch4 + enddo + endif + if (carbox) then + !! Variation of CO ice + dqco=pqsurf(ig,igcm_co_ice)*(zdqsurf(ig)/pqsurf(ig,igcm_n2)*timstep) ! amout of co (kg/m2) to be removed at grid ig (negative) + !! remove CO + pqsurf(ig,igcm_co_ice)=pqsurf(ig,igcm_co_ice)+dqco + !! add CO + do i=inddeb,indfin + pqsurf(i,igcm_co_ice)=pqsurf(i,igcm_co_ice)-masstmp(i)/stock*dqco + enddo + endif + + ! Add N2 + do i=inddeb,indfin + pqsurf(i,igcm_n2)=pqsurf(i,igcm_n2)+masstmp(i)*timstep + enddo + ! Remove N2 + pqsurf(ig,igcm_n2)=pqsurf(ig,igcm_n2)+zdqsurf(ig)*timstep + + endif + !call testconservmass2d(ngrid,pqsurf(:,1)) + + !!!!!!!!!!!!!!!!!!!! NOT THE POLE !!!!!!!!!!!!!!!!!!!!!!!! + ! here: ig = grid point with large amount of ice + ! Loop on each adjacent point + ! looking for adjacent point + ELSE + + !! 1) Compute Diff H0-H1 and mass to be transferred for each adjacent point + DO compt=1,4 + i=edge(compt) + diff=H0-(phisfi0(i)/g+pqsurf(i,igcm_n2)/1000.) ! (m) + if (diff.gt.difflim) then + masstmp(i)=massflow(ig,i,tgla(ig),diff,pqsurf(ig,igcm_n2),timstep,slid) + else + masstmp(i)=0. + endif + totmasstmp=totmasstmp+masstmp(i) + ENDDO + if (totmasstmp.gt.0) then + !! 2) Available mass to be transferred + stock=maxval(masstmp(:))/secu ! kg/m2/s + !! 3) Real amounts of ice to be transferred : + DO compt=1,4 + i=edge(compt) + masstmp(i)=masstmp(i)/totmasstmp*stock*cell_area(ig)/cell_area(i) !kg/m2/s + ENDDO + !! 4) Tendencies + zdqsurf(ig)=-stock + !! 5) Update PQSURF + + ! move CH4 and CO too + if (methane) then + !! Variation of CH4 ice + dqch4=pqsurf(ig,igcm_ch4_ice)*(zdqsurf(ig)/pqsurf(ig,igcm_n2)*timstep) ! amout of ch4 (kg/m2) to be removed at grid ig (negative) + !! remove CH4 + pqsurf(ig,igcm_ch4_ice)=pqsurf(ig,igcm_ch4_ice)+dqch4 + !! add CH4 + DO compt=1,4 + i=edge(compt) + pqsurf(i,igcm_ch4_ice)=pqsurf(i,igcm_ch4_ice)-masstmp(i)/stock*dqch4 + enddo + endif + if (carbox) then + !! Variation of CO ice + dqco=pqsurf(ig,igcm_co_ice)*(zdqsurf(ig)/pqsurf(ig,igcm_n2)*timstep) ! amout of co (kg/m2) to be removed at grid ig (negative) + !! remove CO + pqsurf(ig,igcm_co_ice)=pqsurf(ig,igcm_co_ice)+dqco + !! add CO + DO compt=1,4 + i=edge(compt) + pqsurf(i,igcm_co_ice)=pqsurf(i,igcm_co_ice)-masstmp(i)/stock*dqco + enddo + endif + + ! Add N2 + totmasstmp=0. + DO compt=1,4 + i=edge(compt) + pqsurf(i,igcm_n2)=pqsurf(i,igcm_n2)+masstmp(i)*timstep + totmasstmp=totmasstmp+masstmp(i)*cell_area(i) !! TB22 tmp + enddo + ! Remove N2 + pqsurf(ig,igcm_n2)=pqsurf(ig,igcm_n2)+zdqsurf(ig)*timstep + + endif + + if (pqsurf(ig,igcm_n2).lt.0) then + print*, pqsurf(ig,igcm_n2) + write(*,*) 'bug in spreadglacier' + stop + endif + + ENDIF ! cas + ENDIF ! qlim + + ENDDO ! ig + + end subroutine spreadglacier_paleo + + real function dist_pluto(lat1,lon1,lat2,lon2) + use comcstfi_mod, only: rad + implicit none + real, intent(in) :: lat1,lon1,lat2,lon2 + real dlon,dlat + real a,c + real radi + radi=rad/1000. + + dlon = lon2 - lon1 + dlat = lat2 - lat1 + a = (sin(dlat/2))**2 + cos(lat1) * cos(lat2) * (sin(dlon/2))**2 + c = 2 * atan2( sqrt(a), sqrt(1-a) ) + dist_pluto = radi * c !! km + return + end function dist_pluto + + real function massflow(ig1,ig2,pts,dif,pqs,dt,slid) + !! Mass of ice to be transferred from one grid point to another + use comgeomfi_h + use geometry_mod, only: latitude, longitude, cell_area + use comcstfi_mod, only: g, rad + use tracer_h, only: rho_n2 + + implicit none +#include "dimensions.h" + + real, intent(in) :: pts,dif,pqs,dt + integer, intent(in) :: ig1,ig2,slid + real lref,ac,n,theta,hdelta,ha,tau,qdry,qsl,usl0 + REAL dist_pluto ! function + + lref=dist_pluto(latitude(ig2),longitude(ig2),latitude(ig1),longitude(ig1))*1000. ! m + usl0=6.3e-6 ! m/s + ac=0.005*exp(9.37-422./pts) !0.005*exp(422./45.-422./pts) + n=2.1+0.0155*(pts-45.) + theta=atan(dif/(lref)) + hdelta=pts/20. + ha=pts*pts/8440. !pts*pts/20./422. + qdry=ac*(rho_n2*g*1.e-6)**n*(pqs/1000.)**(n+2)/(n+2)*tan(theta)**(n-1)/(1+tan(theta)*tan(theta))**(n/2.) + qdry=qdry*0.5*exp( (pqs/1000./ha) / ( 1+ pqs /hdelta) ) + qsl=(pqs/1000.)/abs(tan(theta))*usl0*slid + tau=pqs/1000.*lref/abs((qdry+qsl)*tan(theta)) + massflow=pqs*(1.-exp(-dt/tau))/dt + return + end function massflow + + + + Index: trunk/LMDZ.PLUTO/libf/phypluto/spreadglacier_simple.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/spreadglacier_simple.F90 (revision 3412) +++ trunk/LMDZ.PLUTO/libf/phypluto/spreadglacier_simple.F90 (revision 3412) @@ -0,0 +1,122 @@ + subroutine spreadglacier_simple(ngrid,nq,pqsurf,timstep) + + use geometry_mod, only: latitude, longitude, cell_area + use radinc_h, only : naerkind + use surfdat_h, only: phisfi + use tracer_h, only: igcm_n2 + + implicit none + +!================================================================== +! Purpose +! ------- +! Spreading the glacier in a circular crater using a speed of flow +! +! Inputs +! ------ +! ngrid Number of vertical columns +! pqsurf(ngrid,nq) N2 ice tracer on surface kg/m2 +! +! Outputs +! ------- +! pqsurf(ngrid,nq) new value for N2 ice tracer on surface kg/m2 +! +! Authors +! ------- +! Tanguy Bertrand (2015) +! +!================================================================== + +#include "dimensions.h" + +!----------------------------------------------------------------------- +! Arguments + + INTEGER ngrid, nq, ig, ig2 + REAL :: pqsurf(ngrid,nq) + REAL timstep +!----------------------------------------------------------------------- +! Local variables + + REAL :: zdqsurf(ngrid) + REAL phiref + REAL dist + REAL locdist + REAL tau + REAL speed + REAL dist_pluto2 ! function + LOGICAL firstcall + DATA firstcall/.true./ + INTEGER indglac(ngrid) + INTEGER nglac + +!---------------- INPUT ------------------------------------------------ +! Defined for Tombaugh crater : + + dist=300. ! reference radius km + phiref=-1000. ! geopotential reference +! tau=1.e10 + speed=6.4e-7 ! glacier speed km/s +! tau=1. ! direct redistribution glacier + +!--------------- FIRSTCALL : define crater index ----------------------- + IF (firstcall) then + indglac(:)=0 + DO ig=1,ngrid + IF (phisfi(ig).lt.phiref) THEN + nglac=nglac+1 + indglac(nglac)=ig + ENDIF + ENDDO + !write(*,*) 'inglac=',indglac + write(*,*) 'nglac=',nglac + + firstcall=.false. + + ENDIF + +!--------------- redistribution ------------------------------------- + DO ig=1,nglac + ! redistribution if large amount of ice only : detection of reservoirs + IF (pqsurf(indglac(ig),igcm_n2).gt.50.) THEN + ! looking where to redistribute : if less ice, lower altitude + ! and close to the reservoir + DO ig2=1,nglac + IF ((pqsurf(indglac(ig2),igcm_n2).lt.pqsurf(indglac(ig),igcm_n2)*0.8)) THEN + !write(*,*) 'loop=',latitude(indglac(ig2))*180/3.14,longitude(indglac(ig2))*180/3.14 + !.and.(phisfi(indglac(ig2)).le.phisfi(indglac(ig)))) THEN + locdist=dist_pluto2(latitude(indglac(ig2)),longitude(indglac(ig2)),latitude(indglac(ig)),longitude(indglac(ig))) + !write(*,*) 'dist=',locdist + IF (locdist.lt.dist) THEN + !write(*,*) 'calcul=',latitude(indglac(ig2))*180/3.14,longitude(indglac(ig2))*180/3.14,latitude(indglac(ig))*180/3.14,longitude(indglac(ig))*180/3.14 + tau=locdist/speed ! characteristic time for redistribution + zdqsurf(indglac(ig2))=pqsurf(indglac(ig),igcm_n2)*(1.-exp(-timstep/tau))/timstep*cell_area(indglac(ig))/cell_area(indglac(ig2)) + zdqsurf(indglac(ig))=pqsurf(indglac(ig),igcm_n2)*(exp(-timstep/tau)-1.)/timstep + pqsurf(indglac(ig2),igcm_n2)=pqsurf(indglac(ig2),igcm_n2)+zdqsurf(indglac(ig2))*timstep + pqsurf(indglac(ig),igcm_n2)=pqsurf(indglac(ig),igcm_n2)+zdqsurf(indglac(ig))*timstep + ENDIF + ENDIF + ENDDO + ENDIF + ENDDO + + end subroutine spreadglacier_simple + + real function dist_pluto2(lat1,lon1,lat2,lon2) + implicit none + !#include "comcstfi.h" + real, intent(in) :: lat1,lon1,lat2,lon2 + real dlon,dlat + real a,c,d + real rad + rad=1187. ! planetary radius (km) 1187 km sur Pluton + + dlon = lon2 - lon1 + dlat = lat2 - lat1 + a = (sin(dlat/2))**2 + cos(lat1) * cos(lat2) * (sin(dlon/2))**2 + c = 2 * atan2( sqrt(a), sqrt(1-a) ) + dist_pluto2 = rad * c + return + end function dist_pluto2 + + Index: trunk/LMDZ.PLUTO/libf/phypluto/testconserv.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/testconserv.F90 (revision 3412) +++ trunk/LMDZ.PLUTO/libf/phypluto/testconserv.F90 (revision 3412) @@ -0,0 +1,118 @@ + subroutine testconserv(ngrid,nlayer,nq,igcm1,igcm2,ptimestep, & + pplev,zdq,zdqs,car1,car2) + + use comgeomfi_h + use comcstfi_mod, only: pi, g + use geometry_mod, only: latitude, longitude, cell_area + USE tracer_h, only: igcm_co_ice, igcm_ch4_ice + implicit none + +!================================================================== +! Purpose +! ------- +! Test conservation of tracers +! +! Inputs +! ------ +! ngrid Number of vertical columns +! nlayer Number of layers +! pplev(ngrid,nlayer+1) Pressure levels +! zdq +! +! Outputs +! ------- +! +! Both +! ---- +! +! Authors +! ------- +! Tanguy Bertrand 2015 +! +!================================================================== + +#include "dimensions.h" + +!----------------------------------------------------------------------- +! Arguments + + INTEGER ngrid, nlayer, nq + INTEGER igcm1, igcm2 + REAL zdq(ngrid,nlayer,nq) + REAL zdqs(ngrid,nq) + REAL ptimestep + REAL pplev(ngrid,nlayer+1) + character(len=3) :: car1 + character(len=7) :: car2 + + ! LOCAL VARIABLES + INTEGER l,ig,iq + REAL masse + + ! OUTPUT + REAL dWtot + REAL dWtots + REAL nconsMAX + INTEGER myig + REAL vdifcncons(ngrid) +!----------------------------------------------------------------------- + + write(*,*) 'TB igcm=',igcm1,igcm2 + write(*,*) 'TB zdq1=',zdq(1,1,igcm1) + dWtot=0.0 + dWtots=0.0 + nconsMAX=0.0 + do ig = 1, ngrid + vdifcncons(ig)=0.0 + do l = 1, nlayer + masse = (pplev(ig,l) - pplev(ig,l+1))/g + + iq = igcm1 + ! sum of atmospheric mass : kg + dWtot = dWtot + masse*zdq(ig,l,iq)*ptimestep*cell_area(ig) + ! for each column, total mass lost per sec : kg(tracer) / m2 /s + vdifcncons(ig)=vdifcncons(ig) + masse*zdq(ig,l,iq) + + ! if clouds : + IF ((igcm2.eq.igcm_co_ice).or.(igcm2.eq.igcm_ch4_ice)) THEN + iq = igcm2 + dWtot = dWtot + masse*zdq(ig,l,iq)*ptimestep*cell_area(ig) + vdifcncons(ig)=vdifcncons(ig) + masse*zdq(ig,l,iq) + ENDIF + + enddo + + iq = igcm1 + dWtots = dWtots + zdqs(ig,iq)*ptimestep*cell_area(ig) + vdifcncons(ig)=vdifcncons(ig)+zdqs(ig,iq) + + IF ((igcm2.eq.igcm_co_ice).or.(igcm2.eq.igcm_ch4_ice)) THEN + iq = igcm2 + dWtots = dWtots + zdqs(ig,iq)*ptimestep*cell_area(ig) + vdifcncons(ig)=vdifcncons(ig)+zdqs(ig,iq) + ENDIF + + ! vdifcncons is the total amount of material that appear or + ! disapear per second in the routine + ! it is the non conservative factor + + if(vdifcncons(ig).gt.nconsMAX)then + nconsMAX=vdifcncons(ig) + myig=ig + endif + + enddo + write(*,*) 'TB Dwtot=',dWtots,dWtot + + dWtot = dWtot/totarea + dWtots = dWtots/totarea + print*,'-------------------------------------------' + print*,'In ',car2,' atmospheric ',car1,' change=',dWtot,' kg m-2' + print*,'In ',car2,' surface ',car1,' change =',dWtots,' kg m-2' + print*,'--> non-cons factor =',dWtot+dWtots,' kg m-2' + print*,'--> MAX non-cons factor =',nconsMAX,' kg m-2 s-1' + IF (nconsMAX.gt.0.) then + print*,'--> obtained at lat/lon=',latitude(myig)*180./pi,longitude(myig)*180./pi + ENDIF + end subroutine testconserv + Index: trunk/LMDZ.PLUTO/libf/phypluto/testconservfast.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/testconservfast.F90 (revision 3412) +++ trunk/LMDZ.PLUTO/libf/phypluto/testconservfast.F90 (revision 3412) @@ -0,0 +1,95 @@ + subroutine testconservfast(ngrid,nlayer,nq,ptimestep, & + pplev,zdq,zdqs,car1,car2) + + use comgeomfi_h + use comcstfi_mod, only: pi, g + use geometry_mod, only: latitude, longitude, cell_area + implicit none + +!================================================================== +! Purpose +! ------- +! Test conservation of tracers +! +! Inputs +! ------ +! ngrid Number of vertical columns +! nlayer Number of layers +! pplev(ngrid,nlayer+1) Pressure levels +! zdq +! +! Outputs +! ------- +! +! Both +! ---- +! +! Authors +! ------- +! Tanguy Bertrand 2015 +! +!================================================================== + +#include "dimensions.h" + +!----------------------------------------------------------------------- +! Arguments + + INTEGER ngrid, nlayer, nq + REAL zdq(ngrid) + REAL zdqs(ngrid) + REAL ptimestep + REAL pplev(ngrid,nlayer+1) + character(len=3) :: car1 + character(len=7) :: car2 + + ! LOCAL VARIABLES + INTEGER l,ig,iq + REAL masse + + ! OUTPUT + REAL dWtot + REAL dWtots + REAL nconsMAX + INTEGER myig + REAL vdifcncons(ngrid) +!----------------------------------------------------------------------- + + dWtot=0.0 + dWtots=0.0 + nconsMAX=0.0 + do ig = 1, ngrid + vdifcncons(ig)=0.0 + + ! sum of atmospheric mass : kg + dWtot = dWtot + zdq(ig)*ptimestep*cell_area(ig)*pplev(ig,1)/g + ! for each column, total mass lost per sec : kg(tracer) / m2 /s + vdifcncons(ig)=vdifcncons(ig) + pplev(ig,1)/g*zdq(ig) + + dWtots = dWtots + zdqs(ig)*ptimestep*cell_area(ig) + vdifcncons(ig)=vdifcncons(ig)+zdqs(ig) + + + ! vdifcncons is the total amount of material that appear or + ! disapear per second in the routine + ! it is the non conservative factor + + if(vdifcncons(ig).gt.nconsMAX)then + nconsMAX=vdifcncons(ig) + myig=ig + endif + + enddo + + dWtot = dWtot/totarea + dWtots = dWtots/totarea + print*,'-------------------------------------------' + print*,'In ',car2,' atmospheric ',car1,' change=',dWtot,' kg m-2' + print*,'In ',car2,' surface ',car1,' change =',dWtots,' kg m-2' + print*,'--> non-cons factor =',dWtot+dWtots,' kg m-2' + print*,'--> MAX non-cons factor =',nconsMAX,' kg m-2 s-1' + IF (nconsMAX.gt.0.) then + print*,'--> obtained at lat/lon=',latitude(myig)*180./pi,longitude(myig)*180./pi + ENDIF + end subroutine testconservfast + Index: trunk/LMDZ.PLUTO/libf/phypluto/testconservmass.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/testconservmass.F90 (revision 3412) +++ trunk/LMDZ.PLUTO/libf/phypluto/testconservmass.F90 (revision 3412) @@ -0,0 +1,87 @@ + subroutine testconservmass(ngrid,nlayer, & + ps,surfmass) + + use comcstfi_mod, only: g + use geometry_mod, only: cell_area + USE comgeomfi_h, only: totarea + implicit none + +!================================================================== +! Purpose +! ------- +! Test conservation of n2 mass +! +! Inputs +! ------ +! ngrid Number of vertical columns +! nlayer Number of layers +! pplev(ngrid,nlayer+1) Pressure levels +! surfmass +! +! Authors +! ------- +! Tanguy Bertrand 2015 +! +!================================================================== + +#include "dimensions.h" + +!----------------------------------------------------------------------- +! Arguments + INTEGER ngrid, nlayer + REAL ps(ngrid) + REAL surfmass(ngrid) + + ! LOCAL VARIABLES + INTEGER l,ig,iq + REAL atmmass + REAL pstot + + ! OUTPUT + REAL atmmasstot + REAL surfmasstot + REAL totmass +!----------------------------------------------------------------------- + + atmmass=0.0 + atmmasstot=0.0 + surfmasstot=0.0 + pstot=0.0 + + do ig = 2, ngrid-1 + + ! mass atm kg/m2 + atmmass = ps(ig)/g + + pstot = pstot+ps(ig)*cell_area(ig) + + ! mass atm kg + atmmasstot= atmmasstot + atmmass*cell_area(ig) + + surfmasstot= surfmasstot + surfmass(ig)*cell_area(ig) + + enddo + atmmasstot= atmmasstot + ps(1)/g*cell_area(1)*iim + atmmasstot= atmmasstot + ps(ngrid)/g*cell_area(ngrid)*iim + + surfmasstot= surfmasstot + surfmass(1)*cell_area(1)*iim + surfmasstot= surfmasstot + surfmass(ngrid)*cell_area(ngrid)*iim + + pstot= pstot + ps(1)*cell_area(1)*iim + pstot= pstot + ps(ngrid)*cell_area(ngrid)*iim + + atmmasstot= atmmasstot/(totarea+(iim-1)*(cell_area(1)+cell_area(ngrid))) + pstot= pstot/(totarea+(iim-1)*(cell_area(1)+cell_area(ngrid))) + surfmasstot = surfmasstot/(totarea+(iim-1)*(cell_area(1)+cell_area(ngrid))) + + totmass=surfmasstot+atmmasstot + +! print*,'-------------------------------------------' +! print*,'Total mass surface : ',surfmasstot,' kg m-2' +! print*,'Total mass atmosphere :',atmmasstot,' kg m-2' +! print*,'Total mean pressure:',pstot,' Pa' + print*,'Total mass : ',surfmasstot+atmmasstot,' kg m-2' + + + end subroutine testconservmass +