Changeset 3067
- Timestamp:
- Oct 3, 2023, 11:21:28 AM (15 months ago)
- Location:
- trunk/LMDZ.MARS
- Files:
-
- 3 edited
-
changelog.txt (modified) (1 diff)
-
libf/phymars/dyn1d/init_testphys1d_mod.F90 (modified) (9 diffs)
-
libf/phymars/dyn1d/testphys1d.F90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/changelog.txt
r3066 r3067 4231 4231 Slight changes are expected for the PEM (TBC w/ JBC). 4232 4232 4233 == 03/10/2023 == JBC 4234 In 1D, 'q' has been converted from dimension (:,:) to (1,:,:) and 'q2' is now got through the module 'turb_mod'. It allows more generalization and to match dimension in the subroutines. -
trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90
r3066 r3067 5 5 contains 6 6 7 SUBROUTINE init_testphys1d(pem1d,ngrid,nlayer,odpref,nq,q,time,psurf,u,v,temp,startfiles_1D,therestart1D, therestartfi,&8 ndt,ptif,pks,dttestphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,q2,play,plev,&9 latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau)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 10 11 11 use ioipsl_getincom, only: getin ! To use 'getin' … … 13 13 use time_phylmdz_mod, only: daysec, day_step, ecritphy, iphysiq 14 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, &15 use surfdat_h, only: albedodat, z0_default, z0, emissiv, emisice, albedice, iceradius, dtemisice, & 16 16 zmea, zstd, zsig, zgam, zthe, hmons, summit, base, phisfi, watercaptag, watercap, & 17 17 tsurf, emis, qsurf … … 37 37 use mod_grid_phy_lmdz, only: regular_lonlat 38 38 use phys_state_var_init_mod, only: phys_state_var_init 39 use turb_mod, only: q2 39 40 ! Mostly for XIOS outputs: 40 41 use mod_const_mpi, only: COMM_LMDZ … … 48 49 ! Arguments 49 50 !======================================================================= 50 integer, 51 real, 52 logical, 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 53 54 54 55 integer, intent(inout) :: nq 55 56 56 real, dimension(:,:), allocatable, intent(out) :: q ! tracer mixing ratio (e.g. kg/kg) 57 real, intent(out) :: time ! time (0<time<1; time=0.5 at noon) 58 real, intent(out) :: psurf ! Surface pressure 59 real, dimension(nlayer), intent(out) :: u, v ! zonal, meridional wind 60 real, dimension(nlayer), intent(out) :: temp ! temperature at the middle of the layers 61 logical, intent(out) :: startfiles_1D, therestart1D, therestartfi ! Use of starting files for 1D 62 integer, intent(out) :: ndt 63 real, intent(out) :: ptif, pks 64 real, intent(out) :: dttestphys ! testphys1d timestep 65 real, dimension(:), allocatable, intent(out) :: zqsat ! useful for (atm_wat_profile=2) 66 real, dimension(:,:), allocatable, intent(out) :: dq, dqdyn ! Physical and dynamical tandencies 67 integer, intent(out) :: day0 ! initial (sol ; =0 at Ls=0) and final date 68 real, intent(out) :: day ! date during the run 69 real, intent(out) :: gru, grv ! prescribed "geostrophic" background wind 70 real, dimension(nlayer), intent(out) :: w ! "Dummy wind" in 1D 71 real, dimension(nlayer + 1), intent(out) :: q2 ! Turbulent Kinetic Energy 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 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 76 77 77 !======================================================================= … … 237 237 238 238 ! allocate arrays: 239 allocate(tname(nq),q( nlayer,nq),zqsat(nlayer))240 allocate(dq( nlayer,nq),dqdyn(nlayer,nq),tnom_transp(nq))239 allocate(tname(nq),q(1,nlayer,nq),zqsat(nlayer)) 240 allocate(dq(1,nlayer,nq),dqdyn(1,nlayer,nq),tnom_transp(nq)) 241 241 242 242 ! read tracer names from file traceur.def … … 289 289 write(*,*) 'nqfils=',nqfils 290 290 291 292 293 291 #ifdef CPP_XIOS 294 292 call init_physics_distribution(regular_lonlat,4,1,1,1,nlayer,COMM_LMDZ) … … 455 453 do iq = 1,nq 456 454 open(3,file = start1Dname,status = "old",action = "read") 457 read(3,*) header, qsurf(1,iq,1),(q( ilayer,iq), ilayer = 1,nlayer)455 read(3,*) header, qsurf(1,iq,1),(q(1,ilayer,iq), ilayer = 1,nlayer) 458 456 if (trim(tname(iq)) /= trim(header)) then 459 457 write(*,*) 'Tracer names not compatible for initialization with "'//trim(start1Dname)//'"!' … … 554 552 w = 0. ! default: no vertical wind 555 553 556 ! Initialize turbulent ekinetic energy554 ! Initialize turbulent kinetic energy 557 555 q2 = 0. 558 556 … … 690 688 ! q & psurf arrays are on the dynamics scalar grid 691 689 allocate(qdyn(2,1,llm,nq),psdyn(2,1)) 692 qdyn(1,1,1:llm,1:nq) = q(1 :llm,1:nq)690 qdyn(1,1,1:llm,1:nq) = q(1,1:llm,1:nq) 693 691 psdyn(1:2,1) = psurf 694 692 call inichim_newstart(ngrid,nq,qdyn,qsurf(1,:,1),psdyn,flagh2o,flagthermo) 695 q(1 :llm,1:nq) = qdyn(1,1,1:llm,1:nq)693 q(1,1:llm,1:nq) = qdyn(1,1,1:llm,1:nq) 696 694 endif 697 695 -
trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90
r3066 r3067 1 1 PROGRAM testphys1d 2 2 3 use comsoil_h, 4 use surfdat_h, 5 use comslope_mod, 6 use phyredem, 7 use watersat_mod, 8 use tracer_mod, 9 use comcstfi_h, 10 use time_phylmdz_mod, 11 use dimradmars_mod, 12 use dust_param_mod, 13 use comvert_mod, 14 use physiq_mod, 15 use phyetat0_mod, only: phyetat016 use write_output_mod, 17 use init_testphys1d_mod, 3 use comsoil_h, only: inertiedat, inertiesoil, nsoilmx, tsoil 4 use surfdat_h, only: albedodat, perenial_co2ice, watercap, tsurf, emis, qsurf 5 use comslope_mod, only: def_slope, subslope_dist 6 use phyredem, only: physdem0, physdem1 7 use watersat_mod, only: watersat 8 use tracer_mod, only: igcm_h2o_vap, igcm_h2o_ice, noms 9 use comcstfi_h, only: pi, rad, omeg, g, mugaz, rcp, r, cpp 10 use time_phylmdz_mod, only: daysec, day_step 11 use dimradmars_mod, only: tauvis, totcloudfrac, albedo 12 use dust_param_mod, only: tauscaling 13 use comvert_mod, only: ap, bp, aps, bps, pa, preff, sig 14 use physiq_mod, only: physiq 15 use turb_mod, only: q2 16 use write_output_mod, only: write_output 17 use init_testphys1d_mod, only: init_testphys1d 18 18 ! Mostly for XIOS outputs: 19 use mod_const_mpi, 20 use parallel_lmdz, 19 use mod_const_mpi, only: init_const_mpi 20 use parallel_lmdz, only: init_parallel 21 21 22 22 implicit none … … 58 58 ! Declarations 59 59 !-------------------------------------------------------------- 60 integer, parameter :: ngrid = 1 ! (2+(jjm-1)*iim - 1/jjm) 61 integer, parameter :: nlayer = llm 62 real, parameter :: odpref = 610. ! DOD reference pressure (Pa) 63 integer :: unitstart ! unite d'ecriture de "startfi" 64 integer :: ndt, ilayer, ilevel, isoil, idt, iq 65 logical :: firstcall, lastcall 66 integer :: day0 ! initial (sol ; =0 at Ls=0) 67 real :: day ! date during the run 68 real :: time ! time (0<time<1 ; time=0.5 a midi) 69 real :: dttestphys ! testphys1d timestep 70 real, dimension(nlayer) :: play ! Pressure at the middle of the layers (Pa) 71 real, dimension(nlayer + 1) :: plev ! intermediate pressure levels (pa) 72 real :: psurf ! Surface pressure 73 real, dimension(nlayer) :: u, v ! zonal, meridional wind 74 real :: gru, grv ! prescribed "geostrophic" background wind 75 real, dimension(nlayer) :: temp ! temperature at the middle of the layers 76 real, dimension(:,:), allocatable :: q ! tracer mixing ratio (e.g. kg/kg) 77 real, dimension(1) :: wstar = 0. ! Thermals vertical velocity 78 real, dimension(nlayer + 1) :: q2 ! Turbulent Kinetic Energy 60 integer, parameter :: ngrid = 1 ! (2+(jjm-1)*iim - 1/jjm) 61 integer, parameter :: nlayer = llm 62 real, parameter :: odpref = 610. ! DOD reference pressure (Pa) 63 integer :: unitstart ! unite d'ecriture de "startfi" 64 integer :: ndt, ilayer, ilevel, isoil, idt, iq 65 logical :: firstcall, lastcall 66 integer :: day0 ! initial (sol ; =0 at Ls=0) 67 real :: day ! date during the run 68 real :: time ! time (0<time<1 ; time=0.5 a midi) 69 real :: dttestphys ! testphys1d timestep 70 real, dimension(nlayer) :: play ! Pressure at the middle of the layers (Pa) 71 real, dimension(nlayer + 1) :: plev ! intermediate pressure levels (pa) 72 real :: psurf ! Surface pressure 73 real, dimension(nlayer) :: u, v ! zonal, meridional wind 74 real :: gru, grv ! prescribed "geostrophic" background wind 75 real, dimension(nlayer) :: temp ! temperature at the middle of the layers 76 real, dimension(:,:,:), allocatable :: q ! tracer mixing ratio (e.g. kg/kg) 77 real, dimension(1) :: wstar = 0. ! Thermals vertical velocity 79 78 80 79 ! Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s) 81 real, dimension(nlayer) :: du, dv, dtemp, dudyn, dvdyn, dtempdyn82 real, dimension(1) :: dpsurf83 real, dimension(:,: ), allocatable :: dq, dqdyn80 real, dimension(nlayer) :: du, dv, dtemp, dudyn, dvdyn, dtempdyn 81 real, dimension(1) :: dpsurf 82 real, dimension(:,:,:), allocatable :: dq, dqdyn 84 83 85 84 ! Various intermediate variables … … 115 114 !call initcomgeomphy 116 115 117 call init_testphys1d(.false.,ngrid,nlayer,odpref,nq,q,time,psurf,u,v,temp,startfiles_1D,therestart1D, therestartfi,&118 ndt,ptif,pks,dttestphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,q2,play,plev,&119 latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau)116 call init_testphys1d(.false.,ngrid,nlayer,odpref,nq,q,time,psurf,u,v,temp,startfiles_1D,therestart1D, & 117 therestartfi,ndt,ptif,pks,dttestphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, & 118 play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau) 120 119 121 120 ! Write a "startfi" file … … 171 170 ! If atmospheric water is monitored 172 171 if (atm_wat_tau < 0.) then ! Prescribed atm_wat_profile: wet if >0, dry if =0 173 q( :,igcm_h2o_vap) = min(zqsat(:),atm_wat_profile*g/psurf)174 q( :,igcm_h2o_ice) = 0. ! reset h2o ice172 q(1,:,igcm_h2o_vap) = min(zqsat(:),atm_wat_profile*g/psurf) 173 q(1,:,igcm_h2o_ice) = 0. ! reset h2o ice 175 174 else ! Relaxation towards the value atm_wat_profile with relaxation time atm_wat_tau 176 q( :,igcm_h2o_vap) = atm_wat_profile*g/psurf + (q(:,igcm_h2o_vap) - atm_wat_profile*g/psurf)*dexp(-dttestphys/atm_wat_tau)177 q( :,igcm_h2o_vap) = min(zqsat(:),q(:,igcm_h2o_vap))178 q( :,igcm_h2o_ice) = 0. ! reset h2o ice175 q(1,:,igcm_h2o_vap) = atm_wat_profile*g/psurf + (q(1,:,igcm_h2o_vap) - atm_wat_profile*g/psurf)*dexp(-dttestphys/atm_wat_tau) 176 q(1,:,igcm_h2o_vap) = min(zqsat(:),q(1,:,igcm_h2o_vap)) 177 q(1,:,igcm_h2o_ice) = 0. ! reset h2o ice 179 178 endif 180 179 endif … … 195 194 !enddo 196 195 ! For some tests: No coriolis force at equator 197 !if (latitude(1) == 0.) then196 !if (latitude(1) == 0.) then 198 197 du(:) = du(:) + (gru - u(:))/1.e4 199 198 dv(:) = dv(:) + (grv - v(:))/1.e4 … … 222 221 223 222 ! Increment tracers 224 q( :,:) = q(:,:) + dttestphys*dq(:,:)223 q(1,:,:) = q(1,:,:) + dttestphys*dq(1,:,:) 225 224 enddo ! End of time stepping loop (idt=1,ndt) 226 225
Note: See TracChangeset
for help on using the changeset viewer.
