Changeset 4026 for trunk/LMDZ.PLUTO

Timestamp:

Jan 20, 2026, 9:50:44 AM (37 hours ago)

Author:

tbertrand

Message:

PLUTO PCM:
Adding Paleo options for Triton from Pluto.old
TB

Location:

trunk/LMDZ.PLUTO/libf/phypluto

Files:

• phyredem.F90 (modified) (5 diffs)
• physiq_mod.F90 (modified) (6 diffs)

trunk/LMDZ.PLUTO/libf/phypluto/phyredem.F90

@@ -96 +96 @@
   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 +108 @@
   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 +227 @@
   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 +234 @@
   endif ! of if (nq>0)
 
-! Non-orographic gavity waves !AF24: removed
-
 ! close file
       CALL close_restartphy(nid_restart)
@@ -246 +240 @@
 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

trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90

@@ -38 +38 @@
       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 +2174 @@
 
          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 +2186 @@
                            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 +2202 @@
             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 +2212 @@
                ENDDO
             ENDDO
+
             ! Finally ensure conservation of qsurf
             DO iq=1,nq
@@ -2255 +2254 @@
             ! 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, &