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