Index: trunk/LMDZ.PLUTO/libf/phypluto/phyredem.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/phyredem.F90 (revision 3992) +++ trunk/LMDZ.PLUTO/libf/phypluto/phyredem.F90 (revision 4026) @@ -96,6 +96,6 @@ tab_cntrl(4) = time -int(time) ! final time of day -! Informations about Mars, used by dynamics and physics - tab_cntrl(5) = rad ! radius of Mars (m) ~3397200 +! Informations, used by dynamics and physics + tab_cntrl(5) = rad ! radius tab_cntrl(6) = omeg ! rotation rate (rad.s-1) tab_cntrl(7) = g ! gravity (m.s-2) ~3.72 @@ -108,5 +108,5 @@ tab_cntrl(13) = 0. -! Informations about Mars, only for physics +! Informations, only for physics tab_cntrl(14) = year_day ! length of year (sols) ~668.6 tab_cntrl(15) = periastr ! min. star-planet distance (AU) @@ -227,8 +227,4 @@ call put_field(nid_restart,"q2","pbl wind variance",q2) -! cloud fraction and sea ice !AF24: removed - -! !Slab ocean !AF24: removed - ! tracers if (nq>0) then @@ -238,6 +234,4 @@ endif ! of if (nq>0) -! Non-orographic gavity waves !AF24: removed - ! close file CALL close_restartphy(nid_restart) @@ -246,3 +240,256 @@ end subroutine physdem1 +subroutine physdem0pal(filename,lonfi,latfi,nsoil,ngrid,nlay,nq, & + phystep,day_ini,time,airefi, & + alb,pzmea,pzstd,pzsig,pzgam,pzthe, & + oblipal,eccpal,tpalnew,adjustnew, & + phisfipal,peri_daypal) +! create physics restart file and write time-independent variables + use comsoil_h, only: volcapa, mlayer + use geometry_mod, only: cell_area + use surfdat_h, only: zmea, zstd, zsig, zgam, zthe, & + emisice, emissiv, & + iceradius, dtemisice, phisfi + use iostart, only : create_restartphy, open_restartphy, close_restartphy, & + put_var, put_field, length, ldscrpt, ndscrpt + use mod_grid_phy_lmdz, only : klon_glo + use planete_mod, only: year_day, periastr, apoastr, peri_day, & + obliquit, z0, lmixmin, emin_turb + use comcstfi_mod, only: rad, omeg, g, mugaz, rcp + use time_phylmdz_mod, only: daysec + + implicit none + + character(len=*), intent(in) :: filename + real,intent(in) :: lonfi(ngrid) + real,intent(in) :: latfi(ngrid) + integer,intent(in) :: nsoil + integer,intent(in) :: ngrid + integer,intent(in) :: nlay + integer,intent(in) :: nq + real,intent(in) :: phystep + real,intent(in) :: day_ini + real,intent(in) :: time + real,intent(in) :: airefi(ngrid) + real,intent(in) :: alb(ngrid) + real,intent(in) :: pzmea(ngrid) + real,intent(in) :: pzstd(ngrid) + real,intent(in) :: pzsig(ngrid) + real,intent(in) :: pzgam(ngrid) + real,intent(in) :: pzthe(ngrid) + real,intent(in) :: oblipal,peri_daypal,eccpal ! change of obliquity/periday/eccentricity + real,intent(in) :: tpalnew ! change of time for paleo run + real,intent(in) :: adjustnew ! change in N2 ice albedo + real,intent(in) :: phisfipal(ngrid) ! change of geopotential for paleo run + real :: peripal,aphepal ! change of perihelion/aphelion from eccentricity + + character(ndscrpt), dimension(ldscrpt), parameter :: dscrpt_tab_cntrl = (/ & + "(1) Number of atmospheric columns in physics ", & + "(2) Number of atmospheric layers ", & + "(3) Final day ", & + "(4) Final time of day ", & + "(5) Planet radius (m) ", & + "(6) Rotation rate (rad.s-1) ", & + "(7) Gravity (m.s-2) ", & + "(8) Molar mass of the atmosphere (g.mol-1) ", & + "(9) = r/Cp (=kappa in the dynamics) ", & + "(10) Length of a solar day (s) ", & + "(11) Physics time step (s) ", & + "(12) - ", & + "(13) - ", & + "(14) Length of year (in solar days) ", & + "(15) Minimum star-planet distance (AU) ", & + "(16) Maximum star-planet distance (AU) ", & + "(17) Date of periastro (sols since N. spring) ", & + "(18) Obliquity of the planet (deg) ", & + "(19) Default surface roughness (m) ", & + "(20) - ", & + "(21) - ", & + "(22) - ", & + "(23) - ", & + "(24) Emissivity of northern cap ~0.95 ", & + "(25) Emissivity of southern cap ~0.95 ", & + "(26) Emissivity of martian soil ~.95 ", & + "(27) - ", & + "(28) - ", & + "(29) - ", & + "(30) - ", & + "(31) Mean scat radius of CO2 snow (north) ", & + "(32) Mean scat radius of CO2 snow (south) ", & + "(33) Time scale for snow metamorphism (north) ", & + "(34) Time scale for snow metamorphism (south) ", & + "(35) Soil volumetric heat capacity "/) + real :: tab_cntrl(length) ! nb "length=100" defined in iostart module + + write(*,*) "change of eccentricity, ecc=",eccpal + peripal=(periastr+apoastr)/2.*(1-eccpal) + aphepal=(periastr+apoastr)/2.*(1+eccpal) + write(*,*) "check peri old/new=",periastr,peripal + write(*,*) "check aphe old/new=",apoastr,aphepal + write(*,*) "check a old/new=",apoastr+periastr,aphepal+peripal + + ! Create physics start file + call create_restartphy(filename,nid_restart) + +! tab_cntrl() contains run parameters + tab_cntrl(:)=0 ! initialization +!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +! Fill control array tab_cntrl(:) with paramleters for this run +!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +! Informations on the physics grid + tab_cntrl(1) = float(klon_glo) ! number of nodes on physics grid + tab_cntrl(2) = float(nlay) ! number of atmospheric layers + tab_cntrl(3) = day_ini + int(time) ! final day + tab_cntrl(4) = time -int(time) ! final time of day + +! Informations, used by dynamics and physics + tab_cntrl(5) = rad ! radius + tab_cntrl(6) = omeg ! rotation rate (rad.s-1) + tab_cntrl(7) = g ! gravity (m.s-2) ~3.72 + tab_cntrl(8) = mugaz ! Molar mass of the atmosphere (g.mol-1) ~43.49 + tab_cntrl(9) = rcp ! = r/cp ~0.256793 (=kappa dans dynamique) + tab_cntrl(10) = daysec ! length of a sol (s) ~88775 + + tab_cntrl(11) = phystep ! time step in the physics + tab_cntrl(12) = 0. + tab_cntrl(13) = 0. + +! Informations, only for physics + tab_cntrl(14) = year_day ! length of year (sols) ~668.6 + tab_cntrl(15) = peripal ! min. star-planet distance (AU) + tab_cntrl(16) = aphepal ! max. star-planet distance (AU) + tab_cntrl(17) = peri_daypal ! date of periastron (sols since N. spring) + tab_cntrl(18) = oblipal ! Obliquity of the planet (deg) ~23.98 + +! Boundary layer and turbulence + tab_cntrl(19) = z0 ! surface roughness (m) ~0.01 + tab_cntrl(20) = tpalnew + tab_cntrl(21) = adjustnew + +! Optical properties of polar caps and ground emissivity + tab_cntrl(24) = emisice(1) ! Emissivity of northern cap ~0.95 + tab_cntrl(25) = emisice(2) ! Emissivity of southern cap ~0.95 + tab_cntrl(26) = emissiv ! Emissivity of martian soil ~.95 + tab_cntrl(31) = iceradius(1) ! mean scat radius of N2 snow (north) + tab_cntrl(32) = iceradius(2) ! mean scat radius of N2 snow (south) + tab_cntrl(33) = dtemisice(1) ! time scale for snow metamorphism (north) + tab_cntrl(34) = dtemisice(2) ! time scale for snow metamorphism (south) + + tab_cntrl(28) = 0. + tab_cntrl(29) = 0. + tab_cntrl(30) = 0. +! Soil properties: + tab_cntrl(35) = volcapa ! soil volumetric heat capacity + + call put_var(nid_restart,"controle","Control parameters",tab_cntrl) + + ! Write the controle array descriptor + call put_var(nid_restart,"controle_descriptor",& + "Description of control parameters",dscrpt_tab_cntrl) + + ! Write the mid-layer depths + call put_var(nid_restart,"soildepth","Soil mid-layer depth",mlayer) + + ! Write longitudes + call put_field(nid_restart,"longitude","Longitudes of physics grid",lonfi) + + ! Write latitudes + call put_field(nid_restart,"latitude","Latitudes of physics grid",latfi) + + ! Write mesh areas + call put_field(nid_restart,"area","Mesh area",cell_area) + + ! Write surface geopotential + call put_field(nid_restart,"phisfi","Geopotential at the surface",phisfipal) + + ! Write surface albedo + !call put_field(nid_restart,"albedodat","Albedo of bare ground",alb) + + ! Subgrid topogaphy variables + call put_field(nid_restart,"ZMEA","Relief: mean relief",zmea) + call put_field(nid_restart,"ZSTD","Relief: standard deviation",zstd) + call put_field(nid_restart,"ZSIG","Relief: sigma parameter",zsig) + call put_field(nid_restart,"ZGAM","Relief: gamma parameter",zgam) + call put_field(nid_restart,"ZTHE","Relief: theta parameter",zthe) + + ! Close file + call close_restartphy(nid_restart) + +end subroutine physdem0pal + +subroutine physdem1pal(filename,nsoil,ngrid,nlay,nq, & + phystep,time,tsurf,tsoil,inertiesoil, & + emis,alb,q2,qsurf,n2frac, & + oblipal,eccpal,tpalnew,adjustnew,phisfipal,peri_daypal) + ! write time-dependent variable to restart file + use iostart, only : open_restartphy, close_restartphy, & + put_var, put_field + use tracer_h, only: noms + + implicit none + + character(len=*),intent(in) :: filename + integer,intent(in) :: nsoil + integer,intent(in) :: ngrid + integer,intent(in) :: nlay + integer,intent(in) :: nq + real,intent(in) :: phystep + real,intent(in) :: time + real,intent(in) :: tsurf(ngrid) + real,intent(in) :: tsoil(ngrid,nsoil) + real,intent(in) :: emis(ngrid) + real,intent(in) :: alb(ngrid) + real,intent(in) :: n2frac(ngrid) + real,intent(in) :: inertiesoil(ngrid,nsoil) + real,intent(in) :: q2(ngrid,nlay+1) + real,intent(in) :: qsurf(ngrid,nq) + real,intent(in) :: oblipal,peri_daypal,eccpal ! change of obliquity/periday/eccentricity + real,intent(in) :: tpalnew ! change of time for paleo run + real,intent(in) :: adjustnew ! change in N2 ice albedo + real,intent(in) :: phisfipal(ngrid) ! change of geopotential for paleo run + real :: peripal,aphepal ! change of perihelion/aphelion from eccentricity + + integer :: iq + + ! Open file + call open_restartphy(filename, nid_restart) + + ! First variable to write must be "Time", in order to correctly + ! set time counter in file + !call put_var("Time","Temps de simulation",time) + + ! Surface temperature + call put_field(nid_restart,"tsurf","Surface temperature",tsurf) + + ! Soil inertia + call put_field(nid_restart,"inertiedat","Soil thermal inertia",inertiesoil) + + ! Soil temperature + call put_field(nid_restart,"tsoil","Soil temperature",tsoil) + + ! Emissivity of the surface + call put_field(nid_restart,"emis","Surface emissivity",emis) + + ! Albedo of the surface + call put_field(nid_restart,"albedodat","Albedo of bare ground",alb) + + !n2 fraction on the surface + call put_field(nid_restart,"n2frac","N2 ice fraction on the surface",n2frac) + + ! Planetary Boundary Layer + call put_field(nid_restart,"q2","pbl wind variance",q2) + +! tracers + if (nq>0) then + do iq=1,nq + call put_field(nid_restart,noms(iq),"tracer on surface",qsurf(:,iq)) + enddo + endif ! of if (nq>0) + +! close file + CALL close_restartphy(nid_restart) +!$OMP BARRIER + +end subroutine physdem1pal + end module phyredem Index: trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (revision 3992) +++ trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (revision 4026) @@ -38,5 +38,5 @@ use tabfi_mod, only: tab_cntrl_mod use wstats_mod, only: callstats, wstats, mkstats - use phyredem, only: physdem0, physdem1 + use phyredem, only: physdem0, physdem1, physdem0pal, physdem1pal use planetwide_mod, only: planetwide_minval,planetwide_maxval,planetwide_sumval use mod_phys_lmdz_para, only : is_master @@ -2174,12 +2174,11 @@ if (paleo) then - ! time range for tendencies of ice flux qsurfyear + ! time range to be applied for tendencies of ice flux qsurfyear zdt_tot=year_day ! Last year of simulation + ! update new reservoir of ice on the surface masslost(:)=0. massacc(:)=0. - DO ig=1,ngrid - ! update new reservoir of ice on the surface DO iq=1,nq ! kg/m2 to be sublimed or condensed during paleoyears @@ -2187,18 +2186,14 @@ paleoyears*365.25/(zdt_tot*daysec/86400.) - ! special case if we sublime the entire reservoir - !! AF: TODO : fix following lines (real_area), using line below: - ! call planetwide_sumval((-qsurfyear(:,iq)-qsurf(:,iq))*cell_area(:),masslost) - - ! IF (-qsurfyear(ig,iq).gt.qsurf(ig,iq)) THEN - ! masslost(iq)=masslost(iq)+real_area(ig)* & - ! (-qsurfyear(ig,iq)-qsurf(ig,iq)) - ! qsurfyear(ig,iq)=-qsurf(ig,iq) - ! ENDIF - - ! IF (qsurfyear(ig,iq).gt.0.) THEN - ! massacc(iq)=massacc(iq)+real_area(ig)*qsurfyear(ig,iq) - ! ENDIF - + ! special case if we sublime the entire reservoir + IF (-qsurfyear(ig,iq).gt.qsurf(ig,iq)) THEN + masslost(iq)=masslost(iq)+cell_area(ig)* & + (-qsurfyear(ig,iq)-qsurf(ig,iq)) + qsurfyear(ig,iq)=-qsurf(ig,iq) + ENDIF + ! massacc needed to redistribute masslost + IF (qsurfyear(ig,iq).gt.0.) THEN + massacc(iq)=massacc(iq)+cell_area(ig)*qsurfyear(ig,iq) + ENDIF ENDDO @@ -2207,5 +2202,8 @@ DO ig=1,ngrid DO iq=1,nq + ! New reservoir paleo qsurfpal(ig,iq)=qsurf(ig,iq)+qsurfyear(ig,iq) + ! Redistribution of some of the mass lost when entire reservoir has sublimed + ! Accumulated mass is less than calculated because sublimed mass is less IF (qsurfyear(ig,iq).gt.0.) THEN qsurfpal(ig,iq)=qsurfpal(ig,iq)- & @@ -2214,4 +2212,5 @@ ENDDO ENDDO + ! Finally ensure conservation of qsurf DO iq=1,nq @@ -2255,8 +2254,10 @@ ! create restartfi if (ngrid.ne.1) then - print*, "physdem1pal not yet implemented" - stop - !TODO: import this routine from pluto.old - ! call physdem1pal("restartfi.nc",long,lati,nsoilmx,nq, & + call physdem0pal("restartfi.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, & + ptimestep,pdaypal,time_phys,cell_area, & + albedo_bareground,zmea,zstd,zsig,zgam,zthe, & + oblipal,eccpal,tpalnew,adjustnew,phisfipal,peri_daypal) + + !call physdem1pal("restartfi.nc",long,lati,nsoilmx,nq, & ! ptimestep,pdaypal, & ! ztime_restart,tsurf,tsoil,emis,q2,qsurfpal, &