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