Changeset 3056 for trunk


Ignore:
Timestamp:
Sep 27, 2023, 3:36:36 PM (14 months ago)
Author:
jbclement
Message:

Mars PCM 1D:
From now on, in "testphys1d.F90", initialization is entirely done by a new subroutine called 'init_testphys1d' ("init_testphys1d_mod.F90").
JBC

Location:
trunk/LMDZ.MARS
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/changelog.txt

    r3055 r3056  
    42154215Upgrade of "testphys1d" to Fortran 90. Cleaning of the subroutine and minor optimizations of the code.
    42164216Correction of a bug: 'inertiedat(1,1)' was overwritten by 'inertieice'.
     4217From now on, in "testphys1d.F90", initialization is entirely done by a new subroutine called 'init_testphys1d' ("init_testphys1d_mod.F90").
    42174218
    42184219== 27/09/2023 == EM

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

    r3054 r3056  
    11PROGRAM testphys1d
    22
    3 use ioipsl_getincom,          only: getin ! To use 'getin'
    4 use dimphy,                   only: init_dimphy
    5 use mod_grid_phy_lmdz,        only: regular_lonlat
    6 use infotrac,                 only: nqtot, tname, nqperes, nqfils
    7 use comsoil_h,                only: volcapa, layer, mlayer, inertiedat, inertiesoil, nsoilmx, flux_geo
    8 use comgeomfi_h,              only: sinlat, ini_fillgeom
    9 use surfdat_h,                only: albedodat, z0_default, emissiv, emisice, albedice, iceradius,     &
    10                                     dtemisice, z0, zmea, zstd, zsig, zgam, zthe, phisfi, watercaptag, &
    11                                     watercap, hmons, summit, base, perenial_co2ice
    12 use slope_mod,                only: theta_sl, psi_sl
    13 use comslope_mod,             only: def_slope, subslope_dist, def_slope_mean
    14 use phyredem,                 only: physdem0, physdem1
    15 use geometry_mod,             only: init_geometry
    16 use watersat_mod,             only: watersat
    17 use tracer_mod,               only: igcm_h2o_vap, igcm_h2o_ice, igcm_co2,noms
    18 use planete_h,                only: year_day, periheli, aphelie, peri_day, obliquit, emin_turb, lmixmin
    19 use comcstfi_h,               only: pi, rad, omeg, g, mugaz, rcp, r, cpp
    20 use time_phylmdz_mod,         only: daysec, day_step, ecritphy, iphysiq
    21 use dimradmars_mod,           only: tauvis, totcloudfrac
    22 use dust_param_mod,           only: tauscaling
    23 use comvert_mod,              only: ap, bp, aps, bps, pa, preff, sig, presnivs, pseudoalt, scaleheight
    24 use vertical_layers_mod,      only: init_vertical_layers
    25 use logic_mod,                only: hybrid
    26 use physics_distribution_mod, only: init_physics_distribution
    27 use regular_lonlat_mod,       only: init_regular_lonlat
    28 use mod_interface_dyn_phys,   only: init_interface_dyn_phys
    29 use phys_state_var_init_mod,  only: phys_state_var_init
    30 use physiq_mod,               only: physiq
    31 use read_profile_mod,         only: read_profile
    32 use inichim_newstart_mod,     only: inichim_newstart
    33 use phyetat0_mod,             only: phyetat0
    34 use netcdf,                   only: NF90_OPEN, NF90_NOERR, NF90_NOWRITE, nf90_strerror
    35 use iostart,                  only: open_startphy, get_var, close_startphy
    36 use write_output_mod,         only: write_output
     3use comsoil_h,               only: inertiedat, inertiesoil, nsoilmx
     4use surfdat_h,               only: albedodat, perenial_co2ice, watercap
     5use comslope_mod,            only: def_slope, subslope_dist
     6use phyredem,                only: physdem0, physdem1
     7use watersat_mod,            only: watersat
     8use tracer_mod,              only: igcm_h2o_vap, igcm_h2o_ice, noms
     9use comcstfi_h,              only: pi, rad, omeg, g, mugaz, rcp, r, cpp
     10use time_phylmdz_mod,        only: daysec, day_step
     11use dimradmars_mod,          only: tauvis, totcloudfrac
     12use dust_param_mod,          only: tauscaling
     13use comvert_mod,             only: ap, bp, aps, bps, pa, preff, sig
     14use physiq_mod,              only: physiq
     15use phyetat0_mod,            only: phyetat0
     16use write_output_mod,        only: write_output
     17use init_testphys1d_mod,     only: init_testphys1d
    3718! Mostly for XIOS outputs:
    38 use mod_const_mpi,            only: init_const_mpi, COMM_LMDZ
    39 use parallel_lmdz,            only: init_parallel
     19use mod_const_mpi,           only: init_const_mpi
     20use parallel_lmdz,           only: init_parallel
    4021
    4122implicit none
    … …  
    5637!
    5738!   update: 12/06/2003, including chemistry (S. Lebonnois)
    58 !                            and water ice (F. Montmessin)
    59 !           26/09/2023, conversion in F90 (JB Clément)
     39!                       and water ice (F. Montmessin)
     40!           27/09/2023, upgrade to F90 (JB Clément)
    6041!
    6142!=======================================================================
    … …  
    7758! Declarations
    7859!--------------------------------------------------------------
    79 integer, parameter                :: ngrid = 1 ! (2+(jjm-1)*iim - 1/jjm)
     60integer, parameter                :: ngrid = 1     ! (2+(jjm-1)*iim - 1/jjm)
    8061integer, parameter                :: nlayer = llm
    8162real, parameter                   :: odpref = 610. ! DOD reference pressure (Pa)
    82 integer                           :: unitstart ! unite d'ecriture de "startfi"
    83 integer                           :: nlevel, nsoil, ndt
    84 integer                           :: ilayer, ilevel, isoil, idt, iq
     63integer                           :: unitstart     ! unite d'ecriture de "startfi"
     64integer                           :: ndt, ilayer, ilevel, isoil, idt, iq
    8565logical                           :: firstcall, lastcall
    86 integer                           :: day0, dayn ! initial (sol ; =0 at Ls=0) and final date
    87 real                              :: day        ! date during the run
    88 real                              :: time       ! time (0<time<1 ; time=0.5 a midi)
    89 real                              :: dttestphys ! testphys1d timestep
    90 real, dimension(nlayer)           :: play       ! Pressure at the middle of the layers (Pa)
    91 real, dimension(nlayer + 1)       :: plev       ! intermediate pressure levels (pa)
    92 real                              :: psurf      ! Surface pressure
    93 real, dimension(1)                :: tsurf      ! Surface temperature
    94 real, dimension(nlayer)           :: u, v       ! zonal, meridional wind
    95 real                              :: gru, grv   ! prescribed "geostrophic" background wind
    96 real, dimension(nlayer)           :: temp       ! temperature at the middle of the layers
    97 real, dimension(:,:), allocatable :: q          ! tracer mixing ratio (e.g. kg/kg)
    98 real, dimension(:),   allocatable :: qsurf      ! tracer surface budget (e.g. kg.m-2)
    99 real, dimension(nsoilmx)          :: tsoil      ! subsurface soik temperature (K)
    100 real, dimension(1)                :: emis       ! surface layer
    101 real, dimension(1,1)              :: albedo     ! surface albedo
    102 real, dimension(1)                :: wstar = 0. ! Thermals vertical velocity
    103 real, dimension(nlayer + 1)       :: q2         ! Turbulent Kinetic Energy
    104 real, dimension(nlayer)           :: zlay       ! altitude estimee dans les couches (km)
    105 
    106 ! Physical and dynamical tandencies (e.g.  m.s-2, K/s, Pa/s)
     66integer                           :: day0          ! initial (sol ; =0 at Ls=0)
     67real                              :: day           ! date during the run
     68real                              :: time          ! time (0<time<1 ; time=0.5 a midi)
     69real                              :: dttestphys    ! testphys1d timestep
     70real, dimension(nlayer)           :: play          ! Pressure at the middle of the layers (Pa)
     71real, dimension(nlayer + 1)       :: plev          ! intermediate pressure levels (pa)
     72real                              :: psurf         ! Surface pressure
     73real, dimension(1)                :: tsurf         ! Surface temperature
     74real, dimension(nlayer)           :: u, v          ! zonal, meridional wind
     75real                              :: gru, grv      ! prescribed "geostrophic" background wind
     76real, dimension(nlayer)           :: temp          ! temperature at the middle of the layers
     77real, dimension(:,:), allocatable :: q             ! tracer mixing ratio (e.g. kg/kg)
     78real, dimension(:),   allocatable :: qsurf         ! tracer surface budget (e.g. kg.m-2)
     79real, dimension(nsoilmx)          :: tsoil         ! subsurface soik temperature (K)
     80real, dimension(1)                :: emis          ! surface layer
     81real, dimension(1,1)              :: albedo        ! surface albedo
     82real, dimension(1)                :: wstar = 0.    ! Thermals vertical velocity
     83real, dimension(nlayer + 1)       :: q2            ! Turbulent Kinetic Energy
     84
     85! Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s)
    10786real, dimension(nlayer)           :: du, dv, dtemp, dudyn, dvdyn, dtempdyn
    10887real, dimension(1)                :: dpsurf
    109 real, dimension(:,:), allocatable :: dq
    110 real, dimension(:,:), allocatable :: dqdyn
     88real, dimension(:,:), allocatable :: dq, dqdyn
    11189
    11290! Various intermediate variables
    113 integer                   :: flagthermo, flagh2o
    114 real                      :: zls
    115 real, dimension(nlayer)   :: phi, h, s, w
    116 real                      :: pks, ptif
    117 real                      :: qtotinit,qtot
    118 integer                   :: ierr, aslun
    119 real, dimension(0:nlayer) :: tmp1, tmp2
    120 integer                   :: nq = 1 ! number of tracers
    121 real, dimension(1)        :: latitude, longitude, cell_area
    122 
    123 ! Dummy variables along "dynamics scalar grid"
    124 real, dimension(:,:,:,:), allocatable :: qdyn
    125 real, dimension(:,:),     allocatable :: psdyn
     91integer                 :: aslun
     92real                    :: zls, pks, ptif, qtotinit, qtot
     93real, dimension(nlayer) :: phi, h, s, w
     94integer                 :: nq = 1 ! number of tracers
     95real, dimension(1)      :: latitude, longitude, cell_area
    12696
    12797character(len = 2)  :: str2
    … …  
    12999character(len = 44) :: txt
    130100
    131 ! New flag to compute paleo orbital configurations + few variables JN
    132 real    :: halfaxe, excentric, Lsperi
    133 logical :: paleomars
     101! RV & JBC: Use of "start1D.txt" and "startfi.nc" files
     102logical :: startfiles_1D
    134103
    135104! JN & JBC: Force atmospheric water profiles
    136105real                            :: atm_wat_profile, atm_wat_tau
    137106real, dimension(:), allocatable :: zqsat ! useful for (atm_wat_profile=2)
    138 
    139 ! MVals: isotopes as in the dynamics (CRisi)
    140 integer                                        :: ifils, ipere, generation, ierr0
    141 character(len = 30), dimension(:), allocatable :: tnom_transp ! transporting fluid short name
    142 character(len = 80)                            :: line ! to store a line of text
    143 logical                                        :: continu, there
    144 
    145 ! LL: Subsurface geothermal flux
    146 real :: flux_geo_tmp
    147 
    148 ! RV & JBC: Use of "start1D.txt" and "startfi.nc" files
    149 logical                :: startfiles_1D, found, therestart1D, therestartfi
    150 character(len = 30)    :: header
    151 real, dimension(100)   :: tab_cntrl
    152 real, dimension(1,2,1) :: albedo_read ! surface albedo
    153 
    154 ! LL: Possibility to add subsurface ice
    155 real    :: ice_depth          ! depth of the ice table, ice_depth < 0. means no ice
    156 real    :: inertieice = 2100. ! ice thermal inertia
    157 integer :: iref
    158107!=======================================================================
    159108
    … …  
    171120!call initcomgeomphy
    172121
    173 !------------------------------------------------------
    174 ! Loading run parameters from "run.def" file
    175 !------------------------------------------------------
    176 ! check if 'run.def' file is around (otherwise reading parameters
    177 ! from callphys.def via getin() routine won't work.
    178 inquire(file = 'run.def',exist = there)
    179 if (.not. there) then
    180     write(*,*) 'Cannot find required file "run.def"'
    181     write(*,*) '  (which should contain some input parameters along with the following line: INCLUDEDEF=callphys.def)'
    182     write(*,*) ' ... might as well stop here ...'
    183     stop
    184 endif
    185 
    186 write(*,*)'Do you want to use "start1D.txt" and "startfi.nc" files?'
    187 startfiles_1D = .false.
    188 call getin("startfiles_1D",startfiles_1D)
    189 write(*,*) " startfiles_1D = ", startfiles_1D
    190 
    191 if (startfiles_1D) then
    192     inquire(file = 'start1D.txt',exist = therestart1D)
    193     if (.not. therestart1D) then
    194         write(*,*) 'There is no "start1D.txt" file!'
    195         write(*,*) 'Initialization is done with default values.'
    196     endif
    197     inquire(file = 'startfi.nc',exist = therestartfi)
    198     if (.not. therestartfi) then
    199         write(*,*) 'There is no "startfi.nc" file!'
    200         write(*,*) 'Initialization is done with default values.'
    201     endif
    202 endif
    203 
    204 !------------------------------------------------------
    205 ! Prescribed constants to be set here
    206 !------------------------------------------------------
    207 pi = 2.*asin(1.)
    208 
    209 ! Mars planetary constants
    210 ! ------------------------
    211 rad = 3397200.                   ! mars radius (m) ~3397200 m
    212 daysec = 88775.                  ! length of a sol (s) ~88775 s
    213 omeg = 4.*asin(1.)/(daysec)      ! rotation rate (rad.s-1)
    214 g = 3.72                         ! gravity (m.s-2) ~3.72
    215 mugaz = 43.49                    ! atmosphere mola mass (g.mol-1) ~43.49
    216 rcp = .256793                    ! = r/cp ~0.256793
    217 r = 8.314511*1000./mugaz
    218 cpp = r/rcp
    219 year_day = 669                   ! lenght of year (sols) ~668.6
    220 periheli = 206.66                ! minimum sun-mars distance (Mkm) ~206.66
    221 aphelie = 249.22                 ! maximum sun-mars distance (Mkm) ~249.22
    222 halfaxe = 227.94                 ! demi-grand axe de l'ellipse
    223 peri_day = 485.                  ! perihelion date (sols since N. Spring)
    224 obliquit = 25.2                  ! Obliquity (deg) ~25.2
    225 excentric = 0.0934               ! Eccentricity (0.0934)
    226 
    227 ! Planetary Boundary Layer and Turbulence parameters
    228 ! --------------------------------------------------
    229 z0_default = 1.e-2 ! surface roughness (m) ~0.01
    230 emin_turb = 1.e-6  ! minimal turbulent energy ~1.e-8
    231 lmixmin = 30       ! mixing length ~100
    232 
    233 ! cap properties and surface emissivities
    234 ! ---------------------------------------
    235 emissiv = 0.95         ! Bare ground emissivity ~.95
    236 emisice(1) = 0.95      ! Northern cap emissivity
    237 emisice(2) = 0.95      ! Southern cap emisssivity
    238 albedice(1) = 0.5      ! Northern cap albedo
    239 albedice(2) = 0.5      ! Southern cap albedo
    240 iceradius(1) = 100.e-6 ! mean scat radius of CO2 snow (north)
    241 iceradius(2) = 100.e-6 ! mean scat radius of CO2 snow (south)
    242 dtemisice(1) = 2.      ! time scale for snow metamorphism (north)
    243 dtemisice(2) = 2.      ! time scale for snow metamorphism (south)
    244 
    245 ! mesh surface (not a very usefull quantity in 1D)
    246 ! ------------------------------------------------
    247 cell_area(1) = 1.
    248 
    249 ! check if there is a 'traceur.def' file and process it
    250 ! load tracer names from file 'traceur.def'
    251 open(90,file = 'traceur.def',status = 'old',form = 'formatted',iostat = ierr)
    252 if (ierr /= 0) then
    253     write(*,*) 'Cannot find required file "traceur.def"'
    254     write(*,*) ' If you want to run with tracers, I need it'
    255     write(*,*) ' ... might as well stop here ...'
    256     stop
    257 else
    258     write(*,*) "testphys1d: Reading file traceur.def"
    259     ! read number of tracers:
    260     read(90,*,iostat = ierr) nq
    261     nqtot = nq ! set value of nqtot (in infotrac module) as nq
    262     if (ierr /= 0) then
    263         write(*,*) "testphys1d: error reading number of tracers"
    264         write(*,*) "   (first line of traceur.def) "
    265         stop
    266     endif
    267     if (nq < 1) then
    268         write(*,*) "testphys1d: error number of tracers"
    269         write(*,*) "is nq=",nq," but must be >=1!"
    270         stop
    271     endif
    272 endif
    273 ! allocate arrays:
    274 allocate(tname(nq),q(nlayer,nq),zqsat(nlayer),qsurf(nq))
    275 allocate(dq(nlayer,nq),dqdyn(nlayer,nq),tnom_transp(nq))
    276 
    277 ! read tracer names from file traceur.def
    278 do iq = 1,nq
    279     read(90,'(80a)',iostat = ierr) line ! store the line from traceur.def
    280     if (ierr /= 0) then
    281         write(*,*) 'testphys1d: error reading tracer names...'
    282         stop
    283     endif
    284     ! if format is tnom_0, tnom_transp (isotopes)
    285     read(line,*,iostat = ierr0) tname(iq),tnom_transp(iq)
    286     if (ierr0 /= 0) then
    287         read(line,*) tname(iq)
    288         tnom_transp(iq) = 'air'
    289     endif
    290 enddo
    291 close(90)
    292 
    293 ! Isotopes: as in the 3D case we have to determine father/son relations for isotopes and carrying fluid
    294 allocate(nqfils(nqtot))
    295 nqperes = 0
    296 nqfils(:) = 0
    297 do iq = 1,nqtot
    298     if (tnom_transp(iq) == 'air') then
    299         ! ceci est un traceur père
    300         write(*,*) 'Le traceur',iq,', appele ',trim(tname(iq)),', est un pere'
    301         nqperes = nqperes + 1
    302     else !if (tnom_transp(iq) == 'air') then
    303         ! ceci est un fils. Qui est son père?
    304         write(*,*) 'Le traceur',iq,', appele ',trim(tname(iq)),', est un fils'
    305         continu = .true.
    306         ipere = 1
    307         do while (continu)
    308             if (tnom_transp(iq) == tname(ipere)) then
    309                 ! Son père est ipere
    310                 write(*,*) 'Le traceur',iq,'appele ',trim(tname(iq)),' est le fils de ',ipere,'appele ',trim(tname(ipere))
    311                 nqfils(ipere) = nqfils(ipere) + 1
    312                 continu = .false.
    313             else !if (tnom_transp(iq) == tnom_0(ipere)) then
    314                 ipere = ipere + 1
    315                 if (ipere > nqtot) then
    316                     write(*,*) 'Le traceur',iq,'appele ',trim(tname(iq)),', est orpelin.'
    317                     call abort_gcm('infotrac_init','Un traceur est orphelin',1)
    318                 endif !if (ipere > nqtot) then
    319             endif !if (tnom_transp(iq) == tnom_0(ipere)) then
    320         enddo !do while (continu)
    321     endif !if (tnom_transp(iq) == 'air') then
    322 enddo !do iq=1,nqtot
    323 write(*,*) 'nqperes=',nqperes
    324 write(*,*) 'nqfils=',nqfils
    325 
    326 ! Initialize tracers here:
    327 write(*,*) "testphys1d: initializing tracers"
    328 if (.not. therestart1D) then
    329     call read_profile(nq,nlayer,qsurf,q)
    330 else
    331     do iq = 1,nq
    332         open(3,file = 'start1D.txt',status = "old",action = "read")
    333         read(3,*) header, qsurf(iq),(q(ilayer,iq), ilayer = 1,nlayer)
    334         if (trim(tname(iq)) /= trim(header)) then
    335             write(*,*) 'Tracer names not compatible for initialization with "start1D.txt"!'
    336             stop
    337         endif
    338     enddo
    339 endif
    340 
    341 call init_physics_distribution(regular_lonlat,4,1,1,1,nlayer,COMM_LMDZ)
    342 
    343 ! Date and local time at beginning of run
    344 ! ---------------------------------------
    345 if (.not. startfiles_1D) then
    346     ! Date (in sols since spring solstice) at beginning of run
    347     day0 = 0 ! default value for day0
    348     write(*,*) 'Initial date (in martian sols ; =0 at Ls=0)?'
    349     call getin("day0",day0)
    350     day = float(day0)
    351     write(*,*) " day0 = ",day0
    352     ! Local time at beginning of run
    353     time = 0 ! default value for time
    354     write(*,*)'Initial local time (in hours, between 0 and 24)?'
    355     call getin("time",time)
    356     write(*,*)" time = ",time
    357     time = time/24. ! convert time (hours) to fraction of sol
    358 else
    359     call open_startphy("startfi.nc")
    360     call get_var("controle",tab_cntrl,found)
    361     if (.not. found) then
    362         call abort_physic("open_startphy","tabfi: Failed reading <controle> array",1)
    363     else
    364         write(*,*)'tabfi: tab_cntrl',tab_cntrl
    365     endif
    366     day0 = tab_cntrl(3)
    367     day = float(day0)
    368 
    369     call get_var("Time",time,found)
    370     call close_startphy
    371 endif !startfiles_1D
    372 
    373 ! Discretization (Definition of grid and time steps)
    374 ! --------------------------------------------------
    375 nlevel = nlayer + 1
    376 nsoil = nsoilmx
    377 
    378 day_step = 48 ! default value for day_step
    379 write(*,*)'Number of time steps per sol?'
    380 call getin("day_step",day_step)
    381 write(*,*) " day_step = ",day_step
    382 
    383 ecritphy = day_step ! default value for ecritphy, output every time step
    384 
    385 ndt = 10 ! default value for ndt
    386 write(*,*)'Number of sols to run?'
    387 call getin("ndt",ndt)
    388 write(*,*) " ndt = ",ndt
    389 
    390 dayn = day0 + ndt
    391 ndt = ndt*day_step
    392 dttestphys = daysec/day_step
    393 
    394 ! Imposed surface pressure
    395 ! ------------------------
    396 psurf = 610. ! Default value for psurf
    397 write(*,*) 'Surface pressure (Pa)?'
    398 if (.not. therestart1D) then
    399     call getin("psurf",psurf)
    400 else
    401     read(3,*) header, psurf
    402 endif
    403 write(*,*) " psurf = ",psurf
    404 
    405 ! Reference pressures
    406 pa = 20.     ! transition pressure (for hybrid coord.)
    407 preff = 610. ! reference surface pressure
    408 
    409 ! Aerosol properties
    410 ! ------------------
    411 tauvis = 0.2 ! default value for tauvis (dust opacity)
    412 write(*,'("Reference dust opacity at ",f4.0," Pa?")')odpref
    413 call getin("tauvis",tauvis)
    414 write(*,*) " tauvis = ",tauvis
    415 
    416 ! Orbital parameters
    417 ! ------------------
    418 if (.not. startfiles_1D) then
    419     paleomars = .false. ! Default: no water ice reservoir
    420     call getin("paleomars",paleomars)
    421     if (paleomars) then
    422         write(*,*) "paleomars=", paleomars
    423         write(*,*) "Orbital parameters from callphys.def"
    424         write(*,*) "Enter eccentricity & Lsperi"
    425         write(*,*) 'Martian eccentricity (0<e<1)?'
    426         call getin('excentric ',excentric)
    427         write(*,*)"excentric =",excentric
    428         write(*,*) 'Solar longitude of perihelion (0<Ls<360)?'
    429         call getin('Lsperi',Lsperi )
    430         write(*,*)"Lsperi=",Lsperi
    431         Lsperi = Lsperi*pi/180.0 ! Put it in rad for peri_day
    432         periheli = halfaxe*(1 - excentric)
    433         aphelie = halfaxe*(1 + excentric)
    434         call call_dayperi(Lsperi,excentric,peri_day,year_day)
    435         write(*,*) "Corresponding orbital params for GCM"
    436         write(*,*) " periheli = ",periheli
    437         write(*,*) " aphelie = ",aphelie
    438         write(*,*) "date of perihelion (sol)",peri_day
    439     else
    440         write(*,*) "paleomars=", paleomars
    441         write(*,*) "Default present-day orbital parameters"
    442         write(*,*) "Unless specified otherwise"
    443         write(*,*)'Min. distance Sun-Mars (Mkm)?'
    444         call getin("periheli",periheli)
    445         write(*,*) " periheli = ",periheli
    446 
    447         write(*,*)'Max. distance Sun-Mars (Mkm)?'
    448         call getin("aphelie",aphelie)
    449         write(*,*) " aphelie = ",aphelie
    450 
    451         write(*,*)'Day of perihelion?'
    452         call getin("periday",peri_day)
    453         write(*,*) " periday = ",peri_day
    454 
    455         write(*,*)'Obliquity?'
    456         call getin("obliquit",obliquit)
    457         write(*,*) " obliquit = ",obliquit
    458      endif
    459 endif !(.not. startfiles_1D )
    460 
    461 ! Latitude/longitude
    462 ! ------------------
    463 latitude = 0. ! default value for latitude
    464 write(*,*)'latitude (in degrees)?'
    465 call getin("latitude",latitude(1))
    466 write(*,*) " latitude = ",latitude
    467 latitude = latitude*pi/180.
    468 longitude = 0.
    469 longitude = longitude*pi/180.
    470 
    471 ! Some initializations (some of which have already been
    472 ! done above!) and loads parameters set in callphys.def
    473 ! and allocates some arrays
    474 ! Mars possible matter with dttestphys in input and include!!!
    475 ! Initializations below should mimick what is done in iniphysiq for 3D GCM
    476 call init_interface_dyn_phys
    477 call init_regular_lonlat(1,1,longitude,latitude,(/0.,0./),(/0.,0./))
    478 call init_geometry(1,longitude,latitude,(/0.,0.,0.,0./),(/0.,0.,0.,0./),cell_area)
    479 ! Ehouarn: init_vertial_layers called later (because disvert not called yet)
    480 !      call init_vertical_layers(nlayer,preff,scaleheight,
    481 !     &                      ap,bp,aps,bps,presnivs,pseudoalt)
    482 call init_dimphy(1,nlayer) ! Initialize dimphy module
    483 call phys_state_var_init(1,llm,nq,tname,day0,dayn,time,daysec,dttestphys,rad,g,r,cpp,nqperes,nqfils)! MVals: variables isotopes
    484 call ini_fillgeom(1,latitude,longitude,(/1.0/))
    485 call conf_phys(1,llm,nq)
    486 
    487 ! In 1D model physics are called every time step
    488 ! ovverride iphysiq value that has been set by conf_phys
    489 if (iphysiq /= 1) then
    490     write(*,*) "testphys1d: setting iphysiq=1"
    491     iphysiq = 1
    492 endif
    493 
    494 ! Initialize albedo / soil thermal inertia
    495 ! ----------------------------------------
    496 if (.not. startfiles_1D) then
    497     albedodat(1) = 0.2 ! default value for albedodat
    498     write(*,*)'Albedo of bare ground?'
    499     call getin("albedo",albedodat(1))
    500     write(*,*) " albedo = ",albedodat(1)
    501     albedo(1,1) = albedodat(1)
    502 
    503     inertiedat(1,1) = 400 ! default value for inertiedat
    504     write(*,*)'Soil thermal inertia (SI)?'
    505     call getin("inertia",inertiedat(1,1))
    506     write(*,*) " inertia = ",inertiedat(1,1)
    507 
    508     ice_depth = -1 ! default value: no ice
    509     call getin("subsurface_ice_depth",ice_depth)
    510 
    511     z0(1) = z0_default ! default value for roughness
    512     write(*,*) 'Surface roughness length z0 (m)?'
    513     call getin("z0",z0(1))
    514     write(*,*) " z0 = ",z0(1)
    515 endif !(.not. startfiles_1D )
    516 
    517 ! Initialize local slope parameters (only matters if "callslope"
    518 ! is .true. in callphys.def)
    519 ! slope inclination angle (deg) 0: horizontal, 90: vertical
    520 theta_sl(1) = 0. ! default: no inclination
    521 call getin("slope_inclination",theta_sl(1))
    522 ! slope orientation (deg)
    523 ! 0 == Northward, 90 == Eastward, 180 == Southward, 270 == Westward
    524 psi_sl(1) = 0. ! default value
    525 call getin("slope_orientation",psi_sl(1))
    526 
    527 ! Sub-slopes parameters (assuming no sub-slopes distribution for now).
    528 def_slope(1) = -90 ! minimum slope angle
    529 def_slope(2) = 90  ! maximum slope angle
    530 subslope_dist(1,1) = 1 ! fraction of subslopes in mesh
    531 
    532 ! For the gravity wave scheme
    533 ! ---------------------------
    534 zmea(1) = 0.
    535 zstd(1) = 0.
    536 zsig(1) = 0.
    537 zgam(1) = 0.
    538 zthe(1) = 0.
    539 
    540 ! For the slope wind scheme
    541 ! -------------------------
    542 hmons(1) = 0.
    543 write(*,*)'hmons is initialized to ',hmons(1)
    544 summit(1) = 0.
    545 write(*,*)'summit is initialized to ',summit(1)
    546 base(1) = 0.
    547 
    548 ! Default values initializing the coefficients calculated later
    549 ! -------------------------------------------------------------
    550 tauscaling(1) = 1.   ! calculated in aeropacity_mod.F
    551 totcloudfrac(1) = 1. ! calculated in watercloud_mod.F
    552 
    553 ! Specific initializations for "physiq"
    554 ! -------------------------------------
    555 ! surface geopotential is not used (or useful) since in 1D
    556 ! everything is controled by surface pressure
    557 phisfi(1) = 0.
    558 
    559 ! Initialization to take into account prescribed winds
    560 ! ----------------------------------------------------
    561 ptif = 2.*omeg*sinlat(1)
    562 
    563 ! Geostrophic wind
    564 gru = 10. ! default value for gru
    565 write(*,*)'zonal eastward component of the geostrophic wind (m/s)?'
    566 call getin("u",gru)
    567 write(*,*) " u = ",gru
    568 grv = 0. !default value for grv
    569 write(*,*)'meridional northward component of the geostrophic wind (m/s)?'
    570 call getin("v",grv)
    571 write(*,*) " v = ",grv
    572 
    573 ! Initialize winds for first time step
    574 if (.not. therestart1D) then
    575     u(:) = gru
    576     v(:) = grv
    577 else
    578     read(3,*) header, (u(ilayer), ilayer = 1,nlayer)
    579     read(3,*) header, (v(ilayer), ilayer = 1,nlayer)
    580 endif
    581 w = 0. ! default: no vertical wind
    582 
    583 ! Initialize turbulente kinetic energy
    584 q2 = 0.
    585 
    586 ! CO2 ice on the surface
    587 ! ----------------------
    588 ! get the index of co2 tracer (not known at this stage)
    589 igcm_co2 = 0
    590 do iq = 1,nq
    591     if (trim(tname(iq)) == "co2") igcm_co2 = iq
    592 enddo
    593 if (igcm_co2 == 0) then
    594     write(*,*) "testphys1d error, missing co2 tracer!"
    595     stop
    596 endif
    597 
    598 if (.not. startfiles_1D) then
    599     qsurf(igcm_co2) = 0. ! default value for co2ice
    600     write(*,*)'Initial CO2 ice on the surface (kg.m-2)'
    601     call getin("co2ice",qsurf(igcm_co2))
    602     write(*,*) " co2ice = ",qsurf(igcm_co2)
    603 endif !(.not. startfiles_1D )
    604 
    605 ! emissivity
    606 ! ----------
    607 if (.not. startfiles_1D) then
    608     emis = emissiv
    609     if (qsurf(igcm_co2) == 1.) then
    610         emis = emisice(1) ! northern hemisphere
    611         if (latitude(1) < 0) emis = emisice(2) ! southern hemisphere
    612     endif
    613 endif !(.not. startfiles_1D )
    614 
    615 ! Compute pressures and altitudes of atmospheric levels
    616 ! -----------------------------------------------------
    617 ! Vertical Coordinates
    618 ! """"""""""""""""""""
    619 hybrid = .true.
    620 write(*,*)'Hybrid coordinates?'
    621 call getin("hybrid",hybrid)
    622 write(*,*) " hybrid = ", hybrid
    623 
    624 call disvert_noterre
    625 ! Now that disvert has been called, initialize module vertical_layers_mod
    626 call init_vertical_layers(nlayer,preff,scaleheight,ap,bp,aps,bps,presnivs,pseudoalt)
    627 
    628 plev(:) = ap(:) + psurf*bp(:)
    629 play(:) = aps(:) + psurf*bps(:)
    630 zlay(:) = -200.*r*log(play(:)/plev(1))/g
    631 
    632 
    633 ! Initialize temperature profile
    634 ! ------------------------------
    635 pks = psurf**rcp
    636 
    637 ! Altitude in km in profile: divide zlay by 1000
    638 tmp1(0) = 0.
    639 tmp1(1:) = zlay(:)/1000.
    640 
    641 call profile(nlayer + 1,tmp1,tmp2)
    642 
    643 if (.not. therestart1D) then
    644     tsurf = tmp2(0)
    645     temp(:) = tmp2(1:)
    646 else
    647     read(3,*) header, tsurf, (temp(ilayer), ilayer = 1,nlayer)
    648     close(3)
    649 endif
    650 
    651 ! Initialize soil properties and temperature
    652 ! ------------------------------------------
    653 volcapa = 1.e6 ! volumetric heat capacity
    654 
    655 if (.not. startfiles_1D) then
    656     ! Initialize depths
    657     ! -----------------
    658     do isoil = 1,nsoil
    659         layer(isoil) = 2.e-4*(2.**(isoil - 1)) ! layer depth
    660     enddo
    661 
    662     ! Creating the new soil inertia table if there is subsurface ice:
    663     if (ice_depth > 0) then
    664         iref = 1 ! ice/regolith boundary index
    665         if (ice_depth < layer(1)) then
    666             inertiedat(1,1) = sqrt(layer(1)/((ice_depth/inertiedat(1,1)**2) + ((layer(1) - ice_depth)/inertieice**2)))
    667             inertiedat(1,2:) = inertieice
    668         else ! searching for the ice/regolith boundary:
    669             do isoil = 1,nsoil
    670                 if ((ice_depth >= layer(isoil)) .and. (ice_depth < layer(isoil + 1))) then
    671                     iref = isoil + 1
    672                     exit
    673                 endif
    674             enddo
    675             ! We then change the soil inertia table:
    676             inertiedat(1,:iref - 1) = inertiedat(1,1)
    677             ! We compute the transition in layer(iref)
    678             inertiedat(1,iref) = sqrt((layer(iref) - layer(iref - 1))/(((ice_depth - layer(iref - 1))/inertiedat(1,1)**2) + ((layer(iref) - ice_depth)/inertieice**2)))
    679             ! Finally, we compute the underlying ice:
    680             inertiedat(1,iref + 1:) = inertieice
    681         endif ! (ice_depth < layer(1))
    682     else ! ice_depth < 0 all is set to surface thermal inertia
    683         inertiedat(1,:) = inertiedat(1,1) ! soil thermal inertia
    684     endif ! ice_depth > 0
    685 
    686     inertiesoil(1,:,1) = inertiedat(1,:)
    687 
    688     tsoil(:) = tsurf(1) ! soil temperature
    689 endif !(.not. startfiles_1D)
    690 
    691 flux_geo_tmp = 0.
    692 call getin("flux_geo",flux_geo_tmp)
    693 flux_geo(:,:) = flux_geo_tmp
    694 
    695 ! Initialize depths
    696 ! -----------------
    697 do isoil = 0,nsoil - 1
    698     mlayer(isoil) = 2.e-4*(2.**(isoil - 0.5)) ! mid-layer depth
    699     layer(isoil + 1) = 2.e-4*(2.**isoil)      ! layer depth
    700 enddo
    701 
    702 ! Initialize traceurs
    703 ! -------------------
    704 if (photochem .or. callthermos) then
    705     write(*,*) 'Initializing chemical species'
    706     ! flagthermo=0: initialize over all atmospheric layers
    707     flagthermo = 0
    708     ! check if "h2o_vap" has already been initialized
    709     ! (it has been if there is a "profile_h2o_vap" file around)
    710     inquire(file = "profile_h2o_vap",exist = there)
    711     if (there) then
    712         flagh2o = 0 ! 0: do not initialize h2o_vap
    713     else
    714         flagh2o = 1 ! 1: initialize h2o_vap in inichim_newstart
    715     endif
    716 
    717     ! hack to accomodate that inichim_newstart assumes that
    718     ! q & psurf arrays are on the dynamics scalar grid
    719     allocate(qdyn(2,1,llm,nq),psdyn(2,1))
    720     qdyn(1,1,1:llm,1:nq) = q(1:llm,1:nq)
    721     psdyn(1:2,1) = psurf
    722     call inichim_newstart(ngrid,nq,qdyn,qsurf,psdyn,flagh2o,flagthermo)
    723     q(1:llm,1:nq) = qdyn(1,1,1:llm,1:nq)
    724 endif
    725 
    726 ! Check if the surface is a water ice reservoir
    727 ! ---------------------------------------------
    728 if (.not. startfiles_1D) watercap(1,:) = 0 ! Initialize watercap
    729 watercaptag(1) = .false. ! Default: no water ice reservoir
    730 write(*,*)'Water ice cap on ground?'
    731 call getin("watercaptag",watercaptag)
    732 write(*,*) " watercaptag = ",watercaptag
    733 
    734 ! Check if the atmospheric water profile is specified
    735 ! ---------------------------------------------------
    736 ! Adding an option to force atmospheric water values JN
    737 atm_wat_profile = -1. ! Default: free atm wat profile
    738 if (water) then
    739     write(*,*)'Force atmospheric water vapor profile?'
    740     call getin('atm_wat_profile',atm_wat_profile)
    741     write(*,*) 'atm_wat_profile = ', atm_wat_profile
    742     if (abs(atm_wat_profile + 1.) < 1.e-15) then ! if == -1.
    743         write(*,*) 'Free atmospheric water vapor profile'
    744         write(*,*) 'Total water is conserved in the column'
    745     else if (abs(atm_wat_profile) < 1.e-15) then ! if == 0.
    746         write(*,*) 'Dry atmospheric water vapor profile'
    747     else if (0. < atm_wat_profile .and. atm_wat_profile <= 1.) then
    748         write(*,*) 'Prescribed atmospheric water vapor profile'
    749         write(*,*) 'Unless it reaches saturation (maximal value)'
    750     else
    751         write(*,*) 'Water vapor profile value not correct!'
    752         stop
    753     endif
    754 endif
    755 
    756 ! Check if the atmospheric water profile relaxation is specified
    757 ! --------------------------------------------------------------
    758 ! Adding an option to relax atmospheric water values JBC
    759 atm_wat_tau = -1. ! Default: no time relaxation
    760 if (water) then
    761     write(*,*) 'Relax atmospheric water vapor profile?'
    762     call getin('atm_wat_tau',atm_wat_tau)
    763     write(*,*) 'atm_wat_tau = ', atm_wat_tau
    764     if (atm_wat_tau < 0.) then
    765         write(*,*) 'Atmospheric water vapor profile is not relaxed.'
    766     else
    767         if (0. <= atm_wat_profile .and. atm_wat_profile <= 1.) then
    768             write(*,*) 'Relaxed atmospheric water vapor profile towards ', atm_wat_profile
    769             write(*,*) 'Unless it reaches saturation (maximal value)'
    770         else
    771             write(*,*) 'Reference atmospheric water vapor profile not known!'
    772             write(*,*) 'Please, specify atm_wat_profile'
    773             stop
    774         endif
    775     endif
    776 endif
     122call init_testphys1d(ngrid,nq,nlayer,odpref,ndt,ptif,pks,dttestphys,startfiles_1D,q,zqsat,qsurf,dq,dqdyn, &
     123                     day0,day,time,psurf,tsurf,gru,grv,u,v,w,q2,play,plev,tsoil,temp,albedo,emis,         &
     124                     latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau)
    777125
    778126! Write a "startfi" file
    … …  
    912260!***********************************************************************
    913261!***********************************************************************
     262
Note: See TracChangeset for help on using the changeset viewer.