Changeset 3203 for trunk/LMDZ.MARS/libf/phymars/dyn1d

Timestamp:

Feb 7, 2024, 12:14:38 PM (17 months ago)

Author:

jbclement

Message:

Mars PCM:

• Addition of the possibility to use subslopes parametization in 1D. To do so, put 'nslope=x' in "run.def" where x (1, 3, 5 or 7) is the number of slopes you want to, or modify it in an appropriate "startfi.nc". Then, default subslopes definition and distribution are used by the model. The already existing case of using one slope whose inclination and orientation are defined in "run.def" is still possible.
• Some cleanings throughout the code, in particular related to unused variables.

JBC

Location:

trunk/LMDZ.MARS/libf/phymars/dyn1d

Files:

• init_testphys1d_mod.F90 (modified) (22 diffs)
• testphys1d.F90 (modified) (4 diffs)

trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90 ¶

@@ -5 +5 @@
 contains
 
-SUBROUTINE init_testphys1d(start1Dname,startfiname,startfiles_1D,therestart1D,therestartfi,ngrid,nlayer,odpref, &
+SUBROUTINE init_testphys1d(start1Dname,startfiname,therestart1D,therestartfi,ngrid,nlayer,odpref, &
                            nq,q,time,psurf,u,v,temp,ndt,ptif,pks,dttestphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,  &
                            play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau,CO2cond_ps)
@@ -15 +15 @@
 use surfdat_h,                only: albedodat, z0_default, z0, emissiv, emisice, albedice, iceradius, dtemisice,      &
                                     zmea, zstd, zsig, zgam, zthe, hmons, summit, base, phisfi, watercaptag, watercap, &
-                                    tsurf, emis, qsurf, perennial_co2ice
+                                    tsurf, emis, qsurf, perennial_co2ice, ini_surfdat_h_slope_var, end_surfdat_h_slope_var
 use infotrac,                 only: nqtot, tname, nqperes, nqfils
 use read_profile_mod,         only: read_profile
 use iostart,                  only: open_startphy, get_var, close_startphy
 use physics_distribution_mod, only: init_physics_distribution
-use comsoil_h,                only: volcapa, nsoilmx, inertiesoil, inertiedat, layer, mlayer, flux_geo, tsoil, qsoil
+use comsoil_h,                only: volcapa, nsoilmx, inertiesoil, inertiedat, layer, mlayer, flux_geo, tsoil, qsoil, &
+                                    ini_comsoil_h_slope_var, end_comsoil_h_slope_var
 use comvert_mod,              only: ap, bp, aps, bps, pa, preff, presnivs, pseudoalt, scaleheight
-use dimradmars_mod,           only: tauvis, totcloudfrac, albedo
+use dimradmars_mod,           only: tauvis, totcloudfrac, albedo, ini_dimradmars_mod_slope_var, end_dimradmars_mod_slope_var
 use regular_lonlat_mod,       only: init_regular_lonlat
 use mod_interface_dyn_phys,   only: init_interface_dyn_phys
@@ -40 +41 @@
 use nonoro_gwd_ran_mod,       only: du_nonoro_gwd, dv_nonoro_gwd
 use conf_phys_mod,            only: conf_phys
-use paleoclimate_mod,         only: paleoclimate
+use paleoclimate_mod,         only: paleoclimate, ini_paleoclimate_h, end_paleoclimate_h
 use comslope_mod,             only: nslope, subslope_dist, ini_comslope_h, end_comslope_h
 ! Mostly for XIOS outputs:
@@ -54 +55 @@
 !=======================================================================
 integer,      intent(in) :: ngrid, nlayer
-real,         intent(in) :: odpref                                    ! DOD reference pressure (Pa)
-character(*), intent(in) :: start1Dname, startfiname                  ! Name of starting files for 1D
-logical,      intent(in) :: startfiles_1D, therestart1D, therestartfi ! Use of starting files for 1D
+real,         intent(in) :: odpref                     ! DOD reference pressure (Pa)
+character(*), intent(in) :: start1Dname, startfiname   ! Name of starting files for 1D
+logical,      intent(in) :: therestart1D, therestartfi ! Use of starting files for 1D
 
 integer, intent(inout) :: nq
@@ -96 +97 @@
 character(30)          :: header
 real, dimension(100)   :: tab_cntrl
-real, dimension(1,2,1) :: albedo_read ! surface albedo
 
 ! New flag to compute paleo orbital configurations + few variables JN
@@ -103 +103 @@
 
 ! MVals: isotopes as in the dynamics (CRisi)
-integer                                  :: ifils, ipere, generation, ierr0
+integer                                  :: ipere, ierr0
 character(30), dimension(:), allocatable :: tnom_transp ! transporting fluid short name
 character(80)                            :: line        ! to store a line of text
@@ -112 +112 @@
 real    :: inertieice = 2100. ! ice thermal inertia
 integer :: iref
-
-! LL: Subsurface geothermal flux
-real :: flux_geo_tmp
 
 !=======================================================================
@@ -122 +119 @@
 ! Prescribed constants to be set here
 !------------------------------------------------------
-
 ! Mars planetary constants
 ! ------------------------
@@ -147 +143 @@
 lmixmin = 30       ! mixing length ~100
 
-! cap properties and surface emissivities
+! Cap properties and surface emissivities
 ! ---------------------------------------
 emissiv = 0.95         ! Bare ground emissivity ~.95
@@ -159 +155 @@
 dtemisice(2) = 2.      ! time scale for snow metamorphism (south)
 
-! mesh surface (not a very usefull quantity in 1D)
+! Mesh surface (not a very usefull quantity in 1D)
 ! ------------------------------------------------
 cell_area(1) = 1.
 
-! check if there is a 'traceur.def' file and process it
-! load tracer names from file 'traceur.def'
+! Check if there is a 'traceur.def' file and process it
+! Load tracer names from file 'traceur.def'
+! -----------------------------------------------------
 open(90,file = 'traceur.def',status = 'old',form = 'formatted',iostat = ierr)
 if (ierr /= 0) then
@@ -188 +185 @@
 endif
 
-! allocate arrays:
+! Allocate arrays
 allocate(tname(nq),q(1,nlayer,nq),zqsat(nlayer))
 allocate(dq(1,nlayer,nq),dqdyn(1,nlayer,nq),tnom_transp(nq))
 
-! read tracer names from file traceur.def
+! Read tracer names from file traceur.def
 do iq = 1,nq
     read(90,'(80a)',iostat = ierr) line ! store the line from traceur.def
@@ -281 +278 @@
 nsoil = nsoilmx
 
-day_step = 48 ! default value for day_step
+day_step = 48 ! Default value for day_step
 write(*,*)'Number of time steps per sol?'
 call getin("day_step",day_step)
 write(*,*) " day_step = ",day_step
 
-ecritphy = day_step ! default value for ecritphy, output every time step
-
-ndt = 10 ! default value for ndt
+ecritphy = day_step ! Default value for ecritphy, output every time step
+
+ndt = 10 ! Default value for ndt
 write(*,*)'Number of sols to run?'
 call getin("ndt",ndt)
@@ -311 +308 @@
 ! Coefficient to control the surface pressure change
 ! To be defined in "callphys.def" so that the PEM can read it asa well
+! --------------------------------------------------------------------
 CO2cond_ps = 1. ! Default value
 write(*,*) 'Coefficient to control the surface pressure change?'
@@ -320 +318 @@
 
 ! Reference pressures
-pa = 20.     ! transition pressure (for hybrid coord.)
-preff = 610. ! reference surface pressure
+pa = 20.     ! Transition pressure (for hybrid coord.)
+preff = 610. ! Reference surface pressure
 
 ! Aerosol properties
 ! ------------------
-tauvis = 0.2 ! default value for tauvis (dust opacity)
+tauvis = 0.2 ! Default value for tauvis (dust opacity)
 write(*,'("Reference dust opacity at ",f4.0," Pa?")')odpref
 call getin("tauvis",tauvis)
@@ -408 +406 @@
 endif
 
-! Initialize tracers (surface and atmosphere) here:
+! Initialize local slope parameters
+! ---------------------------------
+nslope = 1 ! default: one slope
+if (.not. therestartfi) then
+    ! Number of subgrid scale slopes
+    write(*,*)'Number of slopes?'
+    call getin("nslope",nslope)
+    write(*,*) " nslope = ", nslope
+    call end_comslope_h
+    call ini_comslope_h(1,nslope)
+    select case (nslope)
+        case(1)
+            ! Sub-slopes parameters (minimum/maximun angles)
+            def_slope(1) = -90.
+            def_slope(2) = 90.
+            ! Fraction of subslopes in mesh
+            subslope_dist(1,1) = 1.
+        case(3)
+            ! Sub-slopes parameters (minimum/maximun angles)
+            def_slope(1) = -50.
+            def_slope(2) = -3.
+            def_slope(3) = 3.
+            def_slope(4) = 50.
+            ! Fraction of subslopes in mesh
+            subslope_dist(1,1) = 0.1
+            subslope_dist(1,2) = 0.8
+            subslope_dist(1,3) = 0.1
+        case(5)
+            ! Sub-slopes parameters (minimum/maximun angles)
+            def_slope(1) = -43.
+            def_slope(2) = -9.
+            def_slope(3) = -3.
+            def_slope(4) = 3.
+            def_slope(5) = 9.
+            def_slope(6) = 43.
+            ! Fraction of subslopes in mesh
+            subslope_dist(1,1) = 0.025
+            subslope_dist(1,2) = 0.075
+            subslope_dist(1,3) = 0.8
+            subslope_dist(1,4) = 0.075
+            subslope_dist(1,5) = 0.025
+        case(7)
+            ! Sub-slopes parameters (minimum/maximun angles)
+            def_slope(1) = -43.
+            def_slope(2) = -19.
+            def_slope(3) = -9.
+            def_slope(4) = -3.
+            def_slope(5) = 3.
+            def_slope(6) = 9.
+            def_slope(7) = 19.
+            def_slope(8) = 43.
+            ! Fraction of subslopes in mesh
+            subslope_dist(1,1) = 0.01
+            subslope_dist(1,2) = 0.02
+            subslope_dist(1,3) = 0.06
+            subslope_dist(1,4) = 0.8
+            subslope_dist(1,5) = 0.06
+            subslope_dist(1,6) = 0.02
+            subslope_dist(1,7) = 0.01
+        case default
+            error stop 'Number of slopes not possible: nslope should be 1, 3, 5 or 7!'
+    end select
+
+    call end_surfdat_h_slope_var
+    call ini_surfdat_h_slope_var(1,nq,nslope)
+    call end_comsoil_h_slope_var
+    call ini_comsoil_h_slope_var(1,nslope)
+    call end_dimradmars_mod_slope_var
+    call ini_dimradmars_mod_slope_var(1,nslope)
+    call end_paleoclimate_h
+    call ini_paleoclimate_h(1,nslope)
+
+    ! Fraction of subslopes in mesh
+    !write(*,*)'What subslope distribution do you want to use?'
+    !call getin("slope_distribution",subslope_dist)
+    !write(*,*) " subslope_dist = ", subslope_dist(1,:)
+endif
+
+! Slope inclination angle (deg) 0: horizontal, 90: vertical
+theta_sl = 0. ! default: no inclination
+call getin("slope_inclination",theta_sl)
+
+! Slope orientation (deg)
+! 0 == Northward, 90 == Eastward, 180 == Southward, 270 == Westward
+psi_sl = 0. ! default value
+call getin("slope_orientation",psi_sl)
+
+! Initialize tracers (surface and atmosphere)
+!--------------------------------------------
 write(*,*) "init_testphys1d: initializing tracers"
 if (.not. therestart1D) then
     call read_profile(nq,nlayer,qsurf(1,:,1),q)
+    do iq = 1,nq
+        qsurf(1,iq,:) = qsurf(1,iq,1)
+    enddo
 else
     do iq = 1,nq
@@ -430 +519 @@
     call getin("albedo",albedodat(1))
     write(*,*) " albedo = ",albedodat(1)
-    albedo(1,:,1) = albedodat(1)
+    albedo(1,:,:) = albedodat(1)
 
     inertiedat(1,1) = 400 ! default value for inertiedat
@@ -446 +535 @@
     write(*,*) " z0 = ",z0(1)
 endif !(.not. therestartfi)
-
-! Initialize local slope parameters (only matters if "callslope"
-! is .true. in callphys.def)
-! slope inclination angle (deg) 0: horizontal, 90: vertical
-theta_sl(1) = 0. ! default: no inclination
-call getin("slope_inclination",theta_sl(1))
-! slope orientation (deg)
-! 0 == Northward, 90 == Eastward, 180 == Southward, 270 == Westward
-psi_sl(1) = 0. ! default value
-call getin("slope_orientation",psi_sl(1))
-
-! Sub-slopes parameters (assuming no sub-slopes distribution for now).
-def_slope(1) = -90 ! minimum slope angle
-def_slope(2) = 90  ! maximum slope angle
-subslope_dist(1,1) = 1 ! fraction of subslopes in mesh
 
 ! For the gravity wave scheme
@@ -508 +582 @@
 
 ! Initialize winds for first time step
+!-------------------------------------
 if (.not. therestart1D) then
     u = gru
@@ -518 +593 @@
 
 ! Initialize turbulent kinetic energy
+!------------------------------------
 q2 = 0.
 
 ! CO2 ice on the surface
 ! ----------------------
-
 if (.not. therestartfi) then
-    qsurf(1,igcm_co2,1) = 0. ! default value for co2ice
+    qsurf(1,igcm_co2,:) = 0. ! default value for co2ice
     write(*,*)'Initial CO2 ice on the surface (kg.m-2)'
     call getin("co2ice",qsurf(1,igcm_co2,1))
+    qsurf(1,igcm_co2,:) = qsurf(1,igcm_co2,1)
     write(*,*) " co2ice = ",qsurf(1,igcm_co2,1)
 endif !(.not. therestartfi)
@@ -533 +609 @@
 ! ----------
 if (.not. therestartfi) then
-    emis(:,1) = emissiv
+    emis = emissiv
     if (qsurf(1,igcm_co2,1) == 1.) then
-        emis(:,1) = emisice(1) ! northern hemisphere
-        if (latitude(1) < 0) emis(:,1) = emisice(2) ! southern hemisphere
+        emis = emisice(1) ! northern hemisphere
+        if (latitude(1) < 0) emis = emisice(2) ! southern hemisphere
     endif
 endif !(.not. therestartfi)
@@ -568 +644 @@
 
 if (.not. therestart1D) then
-    tsurf(:,1) = tmp2(0)
+    tsurf = tmp2(0)
     temp = tmp2(1:)
 else
@@ -610 +686 @@
     endif ! ice_depth > 0
 
-    inertiesoil(1,:,1) = inertiedat(1,:)
-    tsoil(:,:,1) = tsurf(1,1) ! soil temperature
+    do isoil = 1,nsoil
+        inertiesoil(1,isoil,:) = inertiedat(1,isoil)
+        tsoil(1,isoil,:) = tsurf(1,:) ! soil temperature
+    enddo
 endif !(.not. therestartfi)
 
-flux_geo_tmp = 0.
-call getin("flux_geo",flux_geo_tmp)
-flux_geo(:,:) = flux_geo_tmp
+! Initialize subsurface geothermal flux
+! -------------------------------------
+flux_geo = 0.
+call getin("flux_geo",flux_geo(1,1))
+write(*,*) " flux_geo = ",flux_geo(1,1)
+flux_geo = flux_geo(1,1)
 
 ! Initialize soil content
-! -----------------
+! -----------------------
 if (.not. therestartfi) qsoil = 0.
 

trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90 ¶

@@ -7 +7 @@
 use watersat_mod,        only: watersat
 use tracer_mod,          only: igcm_h2o_vap, igcm_h2o_ice, igcm_co2, noms
-use comcstfi_h,          only: pi, rad, omeg, g, mugaz, rcp, r, cpp
+use comcstfi_h,          only: pi, g, rcp, cpp
 use time_phylmdz_mod,    only: daysec, day_step
 use dimradmars_mod,      only: tauvis, totcloudfrac, albedo
 use dust_param_mod,      only: tauscaling
-use comvert_mod,         only: ap, bp, aps, bps, pa, preff, sig
+use comvert_mod,         only: ap, bp, aps, bps
 use physiq_mod,          only: physiq
 use turb_mod,            only: q2
@@ -63 +63 @@
 integer, parameter                  :: nlayer = llm
 real, parameter                     :: odpref = 610. ! DOD reference pressure (Pa)
-integer                             :: unitstart     ! unite d'ecriture de "startfi"
-integer                             :: ndt, ilayer, ilevel, isoil, idt, iq
+integer                             :: ndt, ilayer, idt
 logical                             :: firstcall, lastcall
 integer                             :: day0          ! initial (sol ; =0 at Ls=0)
@@ -81 +80 @@
 
 ! Physical and dynamical tendencies (e.g. m.s-2, K/s, Pa/s)
-real, dimension(nlayer)             :: du, dv, dtemp, dudyn, dvdyn, dtempdyn
+real, dimension(nlayer)             :: du, dv, dtemp
 real, dimension(1)                  :: dpsurf
 real, dimension(:,:,:), allocatable :: dq, dqdyn
 
 ! Various intermediate variables
-integer                 :: aslun
-real                    :: zls, pks, ptif, qtotinit, qtot
+real                    :: zls, pks, ptif
 real, dimension(nlayer) :: phi, h, s, w
 integer                 :: nq = 1 ! number of tracers
 real, dimension(1)      :: latitude, longitude, cell_area
-character(2)            :: str2
-character(7)            :: str7
-character(44)           :: txt
 
 ! RV & JBC: Use of starting files for 1D
@@ -148 +143 @@
 endif
 
-call init_testphys1d('start1D.txt','startfi.nc',startfiles_1D,therestart1D,therestartfi,ngrid,nlayer,odpref, &
-                     nq,q,time,psurf,u,v,temp,ndt,ptif,pks,dttestphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,     &
+call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,odpref,nq,q, &
+                     time,psurf,u,v,temp,ndt,ptif,pks,dttestphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, &
                      play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau,CO2cond_ps)