Changeset 3053 for trunk/LMDZ.MARS/libf
- Timestamp:
- Sep 27, 2023, 10:45:21 AM (17 months ago)
- File:
-
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90
r3052 r3053 1 2 PROGRAM testphys1d 3 ! to use 'getin' 4 USE ioipsl_getincom, only: getin 5 use dimphy, only : init_dimphy 6 use mod_grid_phy_lmdz, only : regular_lonlat 7 use infotrac, only: nqtot, tname, nqperes,nqfils 8 use comsoil_h, only: volcapa, layer, mlayer, inertiedat, 9 & inertiesoil,nsoilmx,flux_geo 10 use comgeomfi_h, only: sinlat, ini_fillgeom 11 use surfdat_h, only: albedodat, z0_default, emissiv, emisice, 12 & albedice, iceradius, dtemisice, z0, 13 & zmea, zstd, zsig, zgam, zthe, phisfi, 14 & watercaptag, watercap, hmons, summit, base, 15 & perenial_co2ice 16 use slope_mod, only: theta_sl, psi_sl 17 use comslope_mod, only: def_slope,subslope_dist,def_slope_mean 18 use phyredem, only: physdem0,physdem1 19 use geometry_mod, only: init_geometry 20 use watersat_mod, only: watersat 21 use tracer_mod, only: igcm_h2o_vap,igcm_h2o_ice,igcm_co2,noms 22 use planete_h, only: year_day, periheli, aphelie, peri_day, 23 & obliquit, emin_turb, lmixmin 24 use comcstfi_h, only: pi, rad, omeg, g, mugaz, rcp, r, cpp 25 use time_phylmdz_mod, only: daysec, day_step, ecritphy, iphysiq 26 use dimradmars_mod, only: tauvis,totcloudfrac 27 use dust_param_mod, only: tauscaling 28 USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff,sig, 29 & presnivs,pseudoalt,scaleheight 30 USE vertical_layers_mod, ONLY: init_vertical_layers 31 USE logic_mod, ONLY: hybrid 32 use physics_distribution_mod, only: init_physics_distribution 33 use regular_lonlat_mod, only: init_regular_lonlat 34 use mod_interface_dyn_phys, only: init_interface_dyn_phys 35 USE phys_state_var_init_mod, ONLY: phys_state_var_init 36 USE physiq_mod, ONLY: physiq 37 USE read_profile_mod, ONLY: read_profile 38 use inichim_newstart_mod, only: inichim_newstart 39 use phyetat0_mod, only: phyetat0 40 USE netcdf, only: NF90_OPEN, NF90_NOERR, NF90_NOWRITE, 41 & nf90_strerror 42 use iostart, only: open_startphy,get_var, close_startphy 43 use write_output_mod, only: write_output 44 ! mostly for XIOS outputs: 45 use mod_const_mpi, only: init_const_mpi, COMM_LMDZ 46 use parallel_lmdz, only: init_parallel 47 48 IMPLICIT NONE 49 50 c======================================================================= 51 c subject: 52 c -------- 53 c PROGRAM useful to run physical part of the martian GCM in a 1D column 54 c 55 c Can be compiled with a command like (e.g. for 25 layers) 56 c "makegcm -p mars -d 25 testphys1d" 57 c It requires the files "testphys1d.def" "callphys.def" 58 c and a 'run.def' file (containing a "INCLUDEDEF=callphys.def" line) 59 c and a file describing the sigma layers (e.g. "z2sig.def") 60 c 61 c author: Frederic Hourdin, R.Fournier,F.Forget 62 c ------- 63 c 64 c update: 12/06/2003 including chemistry (S. Lebonnois) 65 c and water ice (F. Montmessin) 66 c 67 c======================================================================= 68 69 include "dimensions.h" 70 integer, parameter :: ngrid = 1 !(2+(jjm-1)*iim - 1/jjm) 71 integer, parameter :: nlayer = llm 1 PROGRAM testphys1d 2 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 37 ! Mostly for XIOS outputs: 38 use mod_const_mpi, only: init_const_mpi, COMM_LMDZ 39 use parallel_lmdz, only: init_parallel 40 41 implicit none 42 43 !======================================================================= 44 ! subject: 45 ! -------- 46 ! PROGRAM useful to run physical part of the martian GCM in a 1D column 47 ! 48 ! Can be compiled with a command like (e.g. for 25 layers) 49 ! "makegcm -p mars -d 25 testphys1d" 50 ! It requires the files "testphys1d.def" "callphys.def" 51 ! and a 'run.def' file (containing a "INCLUDEDEF=callphys.def" line) 52 ! and a file describing the sigma layers (e.g. "z2sig.def") 53 ! 54 ! author: Frederic Hourdin, R.Fournier,F.Forget 55 ! ------- 56 ! 57 ! update: 12/06/2003, including chemistry (S. Lebonnois) 58 ! and water ice (F. Montmessin) 59 ! 26/09/2023, conversion in F90 (JB Clément) 60 ! 61 !======================================================================= 62 63 include "dimensions.h" 72 64 !#include "dimradmars.h" 73 65 !#include "comgeomfi.h" … … 76 68 !#include "comsoil.h" 77 69 !#include "comdiurn.h" 78 70 include "callkeys.h" 79 71 !#include "comsaison.h" 80 72 !#include "control.h" 81 73 include "netcdf.inc" 82 74 !#include "advtrac.h" 83 75 84 c -------------------------------------------------------------- 85 c Declarations 86 c -------------------------------------------------------------- 87 c 88 INTEGER unitstart ! unite d'ecriture de "startfi" 89 INTEGER nlevel,nsoil,ndt 90 INTEGER ilayer,ilevel,isoil,idt,iq 91 LOGICAl firstcall,lastcall 92 c 93 real,parameter :: odpref=610. ! DOD reference pressure (Pa) 94 c 95 INTEGER day0,dayn ! initial (sol ; =0 at Ls=0) and final date 96 REAL day ! date during the run 97 REAL time ! time (0<time<1 ; time=0.5 a midi) 98 REAL dttestphys ! testphys1d timestep 99 REAL play(nlayer) ! Pressure at the middle of the layers (Pa) 100 REAL plev(nlayer+1) ! intermediate pressure levels (pa) 101 REAL psurf,tsurf(1) 102 REAL u(nlayer),v(nlayer) ! zonal, meridional wind 103 REAL gru,grv ! prescribed "geostrophic" background wind 104 REAL temp(nlayer) ! temperature at the middle of the layers 105 REAL,ALLOCATABLE :: q(:,:) ! tracer mixing ratio (e.g. kg/kg) 106 REAL,ALLOCATABLE :: qsurf(:) ! tracer surface budget (e.g. kg.m-2) 107 REAL tsoil(nsoilmx) ! subsurface soik temperature (K) 108 REAL emis(1) ! surface layer 109 REAL albedo(1,1) ! surface albedo 110 REAL :: wstar(1)=0. ! Thermals vertical velocity 111 REAL q2(nlayer+1) ! Turbulent Kinetic Energy 112 REAL zlay(nlayer) ! altitude estimee dans les couches (km) 113 114 c Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s) 115 REAL du(nlayer),dv(nlayer),dtemp(nlayer) 116 REAL dudyn(nlayer),dvdyn(nlayer),dtempdyn(nlayer) 117 REAL dpsurf(1) 118 REAL,ALLOCATABLE :: dq(:,:) 119 REAL,ALLOCATABLE :: dqdyn(:,:) 120 121 c Various intermediate variables 122 INTEGER flagthermo,flagh2o 123 REAL zls 124 REAL phi(nlayer),h(nlayer),s(nlayer) 125 REAL pks, ptif, w(nlayer) 126 REAL qtotinit,qtot 127 INTEGER ierr, aslun 128 REAL tmp1(0:nlayer),tmp2(0:nlayer) 129 integer :: nq=1 ! number of tracers 130 real :: latitude(1), longitude(1), cell_area(1) 131 ! dummy variables along "dynamics scalar grid" 132 real,allocatable :: qdyn(:,:,:,:),psdyn(:,:) 133 134 character*2 str2 135 character (len=7) :: str7 136 character(len=44) :: txt 137 138 c New flag to compute paleo orbital configurations + few variables JN 139 REAL halfaxe, excentric, Lsperi 140 Logical paleomars 141 142 c JN & JBC: Force atmospheric water profiles 143 real :: atm_wat_profile, atm_wat_tau 144 real, dimension(:), allocatable :: zqsat ! useful for (atm_wat_profile=2) 145 146 c MVals: isotopes as in the dynamics (CRisi) 147 INTEGER :: ifils,ipere,generation 148 CHARACTER(len=30), ALLOCATABLE, DIMENSION(:) :: tnom_transp ! transporting fluid short name 149 CHARACTER(len=80) :: line ! to store a line of text 150 INTEGER ierr0 151 LOGICAL :: continu 152 logical :: present 153 154 c LL: Subsurface geothermal flux 155 real :: flux_geo_tmp 156 157 c RV & JBC: Use of "start1D.txt" and "startfi.nc" files 158 logical :: startfiles_1D, found 159 logical :: therestart1D, therestartfi 160 character(len = 30) :: header 161 real, dimension(100) :: tab_cntrl 162 real, dimension(1,2,1) :: albedo_read(1,2,1) ! surface albedo 163 164 c LL: Possibility to add subsurface ice 165 REAL :: ice_depth ! depth of the ice table, ice_depth < 0. means no ice 166 REAL :: inertieice = 2100. ! ice thermal inertia 167 integer :: iref 168 169 c======================================================================= 170 171 c======================================================================= 172 c INITIALISATION 173 c======================================================================= 76 !-------------------------------------------------------------- 77 ! Declarations 78 !-------------------------------------------------------------- 79 integer, parameter :: ngrid = 1 ! (2+(jjm-1)*iim - 1/jjm) 80 integer, parameter :: nlayer = llm 81 real, 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 85 logical :: 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) 107 real, dimension(nlayer) :: du, dv, dtemp, dudyn, dvdyn, dtempdyn 108 real, dimension(1) :: dpsurf 109 real, dimension(:,:), allocatable :: dq 110 real, dimension(:,:), allocatable :: dqdyn 111 112 ! 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 126 127 character(len = 2) :: str2 128 character(len = 7) :: str7 129 character(len = 44) :: txt 130 131 ! New flag to compute paleo orbital configurations + few variables JN 132 real :: halfaxe, excentric, Lsperi 133 logical :: paleomars 134 135 ! JN & JBC: Force atmospheric water profiles 136 real :: atm_wat_profile, atm_wat_tau 137 real, 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 158 !======================================================================= 159 160 !======================================================================= 161 ! INITIALISATION 162 !======================================================================= 174 163 #ifdef CPP_XIOS 175 CALLinit_const_mpi176 CALLinit_parallel164 call init_const_mpi 165 call init_parallel 177 166 #endif 178 167 179 ! initialize "serial/parallel" related stuff180 ! CALLinit_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))181 ! CALLinit_phys_lmdz(1,1,llm,1,(/1/))182 ! 183 184 c------------------------------------------------------185 cLoading run parameters from "run.def" file186 c------------------------------------------------------168 ! Initialize "serial/parallel" related stuff 169 !call init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/)) 170 !call init_phys_lmdz(1,1,llm,1,(/1/)) 171 !call initcomgeomphy 172 173 !------------------------------------------------------ 174 ! Loading run parameters from "run.def" file 175 !------------------------------------------------------ 187 176 ! check if 'run.def' file is around (otherwise reading parameters 188 177 ! from callphys.def via getin() routine won't work. 189 open(99,file='run.def',status='old',form='formatted', 190 & iostat=ierr) 191 if (ierr.ne.0) then 192 write(*,*) 'Cannot find required file "run.def"' 193 write(*,*) ' (which should contain some input parameters' 194 write(*,*) ' along with the following line:' 195 write(*,*) ' INCLUDEDEF=callphys.def' 196 write(*,*) ' )' 197 write(*,*) ' ... might as well stop here ...' 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) " 198 265 stop 199 else 200 close(99) 201 endif 202 203 write(*,*)'Do you want to use "start1D.txt" and "startfi.nc"', 204 &' files?' 205 startfiles_1D = .false. 206 call getin("startfiles_1D",startfiles_1D) 207 write(*,*) " startfiles_1D = ", startfiles_1D 208 209 if (startfiles_1D) then 210 inquire(file ='start1D.txt',exist = therestart1D) 211 if (.not. therestart1D) then 212 write(*,*) 'There is no "start1D.txt" file!' 213 write(*,*) 'Initialization is done with default values.' 214 endif 215 inquire(file ='startfi.nc',exist = therestartfi) 216 if (.not. therestartfi) then 217 write(*,*) 'There is no "startfi.nc" file!' 218 write(*,*) 'Initialization is done with default values.' 219 endif 220 endif 221 222 c ------------------------------------------------------ 223 c Prescribed constants to be set here 224 c ------------------------------------------------------ 225 pi=2.E+0*asin(1.E+0) 226 227 c Mars planetary constants 228 c ---------------------------- 229 rad=3397200. ! mars radius (m) ~3397200 m 230 daysec=88775. ! length of a sol (s) ~88775 s 231 omeg=4.*asin(1.)/(daysec) ! rotation rate (rad.s-1) 232 g=3.72 ! gravity (m.s-2) ~3.72 233 mugaz=43.49 ! atmosphere mola mass (g.mol-1) ~43.49 234 rcp=.256793 ! = r/cp ~0.256793 235 r= 8.314511E+0 *1000.E+0/mugaz 236 cpp= r/rcp 237 year_day = 669 ! lenght of year (sols) ~668.6 238 periheli = 206.66 ! minimum sun-mars distance (Mkm) ~206.66 239 aphelie = 249.22 ! maximum sun-mars distance (Mkm) ~249.22 240 halfaxe = 227.94 ! demi-grand axe de l'ellipse 241 peri_day = 485. ! perihelion date (sols since N. Spring) 242 obliquit = 25.2 ! Obliquity (deg) ~25.2 243 excentric = 0.0934 ! Eccentricity (0.0934) 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 244 408 245 c Planetary Boundary Layer and Turbulence parameters 246 c -------------------------------------------------- 247 z0_default = 1.e-2 ! surface roughness (m) ~0.01 248 emin_turb = 1.e-6 ! minimal turbulent energy ~1.e-8 249 lmixmin = 30 ! mixing length ~100 250 251 c cap properties and surface emissivities 252 c ---------------------------------------------------- 253 emissiv= 0.95 ! Bare ground emissivity ~.95 254 emisice(1)=0.95 ! Northern cap emissivity 255 emisice(2)=0.95 ! Southern cap emisssivity 256 albedice(1)=0.5 ! Northern cap albedo 257 albedice(2)=0.5 ! Southern cap albedo 258 iceradius(1) = 100.e-6 ! mean scat radius of CO2 snow (north) 259 iceradius(2) = 100.e-6 ! mean scat radius of CO2 snow (south) 260 dtemisice(1) = 2. ! time scale for snow metamorphism (north) 261 dtemisice(2) = 2. ! time scale for snow metamorphism (south) 262 263 c mesh surface (not a very usefull quantity in 1D) 264 c ---------------------------------------------------- 265 cell_area(1)=1.E+0 266 267 ! check if there is a 'traceur.def' file 268 ! and process it. 269 ! load tracer names from file 'traceur.def' 270 open(90,file='traceur.def',status='old',form='formatted', 271 & iostat=ierr) 272 if (ierr.ne.0) then 273 write(*,*) 'Cannot find required file "traceur.def"' 274 write(*,*) ' If you want to run with tracers, I need it' 275 write(*,*) ' ... might as well stop here ...' 276 stop 277 else 278 write(*,*) "testphys1d: Reading file traceur.def" 279 ! read number of tracers: 280 read(90,*,iostat=ierr) nq 281 nqtot=nq ! set value of nqtot (in infotrac module) as nq 282 if (ierr.ne.0) then 283 write(*,*) "testphys1d: error reading number of tracers" 284 write(*,*) " (first line of traceur.def) " 285 stop 286 endif 287 if (nq<1) then 288 write(*,*) "testphys1d: error number of tracers" 289 write(*,*) "is nq=",nq," but must be >=1!" 290 stop 291 endif 292 endif 293 ! allocate arrays: 294 allocate(tname(nq)) 295 allocate(q(nlayer,nq)) 296 allocate(zqsat(nlayer)) 297 allocate(qsurf(nq)) 298 allocate(dq(nlayer,nq)) 299 allocate(dqdyn(nlayer,nq)) 300 allocate(tnom_transp(nq)) 301 302 ! read tracer names from file traceur.def 303 do iq=1,nq 304 read(90,'(80a)',iostat=ierr) line ! store the line from traceur.def 305 if (ierr.ne.0) then 306 write(*,*) 'testphys1d: error reading tracer names...' 307 stop 308 endif 309 ! if format is tnom_0, tnom_transp (isotopes) 310 read(line,*,iostat=ierr0) tname(iq),tnom_transp(iq) 311 if (ierr0.ne.0) then 312 read(line,*) tname(iq) 313 tnom_transp(iq)='air' 314 endif 315 316 enddo 317 close(90) 318 319 ! Isotopes: as in the 3D case we have to determine father/son relations for isotopes and carrying fluid 320 ALLOCATE(nqfils(nqtot)) 321 nqperes=0 322 nqfils(:)=0 323 DO iq=1,nqtot 324 if (tnom_transp(iq) == 'air') then 325 ! ceci est un traceur père 326 WRITE(*,*) 'Le traceur',iq,', appele ', 327 & trim(tname(iq)),', est un pere' 328 nqperes=nqperes+1 329 else !if (tnom_transp(iq) == 'air') then 330 ! ceci est un fils. Qui est son père? 331 WRITE(*,*) 'Le traceur',iq,', appele ', 332 & trim(tname(iq)),', est un fils' 333 continu=.true. 334 ipere=1 335 do while (continu) 336 if (tnom_transp(iq) .eq. tname(ipere)) then 337 ! Son père est ipere 338 WRITE(*,*) 'Le traceur',iq,'appele ', 339 & trim(tname(iq)),' est le fils de ', 340 & ipere,'appele ',trim(tname(ipere)) 341 nqfils(ipere)=nqfils(ipere)+1 342 continu=.false. 343 else !if (tnom_transp(iq) == tnom_0(ipere)) then 344 ipere=ipere+1 345 if (ipere.gt.nqtot) then 346 WRITE(*,*) 'Le traceur',iq,'appele ', 347 & trim(tname(iq)),', est orpelin.' 348 CALL abort_gcm('infotrac_init', 349 & 'Un traceur est orphelin',1) 350 endif !if (ipere.gt.nqtot) then 351 endif !if (tnom_transp(iq) == tnom_0(ipere)) then 352 enddo !do while (continu) 353 endif !if (tnom_transp(iq) == 'air') then 354 enddo !DO iq=1,nqtot 355 WRITE(*,*) 'nqperes=',nqperes 356 WRITE(*,*) 'nqfils=',nqfils 357 358 ! Initialize tracers here: 359 write(*,*) "testphys1d: initializing tracers" 360 if (.not. therestart1D) then 361 call read_profile(nq,nlayer,qsurf,q) 362 else 363 do iq = 1,nq 364 open(3,file = 'start1D.txt',status = "old", 365 &action = "read") 366 read(3,*) header, qsurf(iq), 367 & (q(ilayer,iq), ilayer = 1,nlayer) 368 if (tname(iq) /= trim(header)) then 369 write(*,*) 'Tracer names not compatible for', 370 &' initialization with "start1D.txt"!' 371 stop 372 endif 373 enddo 374 endif 375 376 call init_physics_distribution(regular_lonlat,4, 377 & 1,1,1,nlayer,COMM_LMDZ) 378 379 c Date and local time at beginning of run 380 c --------------------------------------- 381 if (.not. startfiles_1D) then 382 c Date (in sols since spring solstice) at beginning of run 383 day0 = 0 ! default value for day0 384 write(*,*) 'Initial date (in martian sols ; =0 at Ls=0)?' 385 call getin("day0",day0) 386 day=float(day0) 387 write(*,*) " day0 = ",day0 388 c Local time at beginning of run 389 time=0 ! default value for time 390 write(*,*)'Initial local time (in hours, between 0 and 24)?' 391 call getin("time",time) 392 write(*,*)" time = ",time 393 time=time/24.E+0 ! convert time (hours) to fraction of sol 394 395 else 396 call open_startphy("startfi.nc") 397 call get_var("controle",tab_cntrl,found) 398 if (.not.found) then 399 call abort_physic("open_startphy", 400 & "tabfi: Failed reading <controle> array",1) 401 else 402 write(*,*)'tabfi: tab_cntrl',tab_cntrl 403 endif 404 day0 = tab_cntrl(3) 405 day=float(day0) 406 407 call get_var("Time",time,found) 408 409 call close_startphy 410 411 endif !startfiles_1D 412 413 c Discretization (Definition of grid and time steps) 414 c -------------- 415 nlevel=nlayer+1 416 nsoil=nsoilmx 417 418 day_step=48 ! default value for day_step 419 write(*,*)'Number of time steps per sol ?' 420 call getin("day_step",day_step) 421 write(*,*) " day_step = ",day_step 422 423 ecritphy=day_step ! default value for ecritphy, output every time step 424 425 ndt=10 ! default value for ndt 426 write(*,*)'Number of sols to run ?' 427 call getin("ndt",ndt) 428 write(*,*) " ndt = ",ndt 429 430 dayn=day0+ndt 431 ndt=ndt*day_step 432 dttestphys=daysec/day_step 433 434 c Imposed surface pressure 435 c ------------------------------------ 436 psurf = 610. ! Default value for psurf 437 write(*,*) 'Surface pressure (Pa)?' 438 if (.not. therestart1D) then 439 call getin("psurf",psurf) 440 else 441 read(3,*) header, psurf 442 endif 443 write(*,*) " psurf = ",psurf 444 445 c Reference pressures 446 pa=20. ! transition pressure (for hybrid coord.) 447 preff=610. ! reference surface pressure 448 449 c Aerosol properties 450 c -------------------------------- 451 tauvis=0.2 ! default value for tauvis (dust opacity) 452 write(*,'("Reference dust opacity at ",f4.0," Pa ?")')odpref 453 call getin("tauvis",tauvis) 454 write(*,*) " tauvis = ",tauvis 455 456 c Orbital parameters 457 c ------------------ 458 if (.not. startfiles_1D) then 459 paleomars=.false. ! Default: no water ice reservoir 460 call getin("paleomars",paleomars) 461 if (paleomars) then 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 462 422 write(*,*) "paleomars=", paleomars 463 423 write(*,*) "Orbital parameters from callphys.def" 464 424 write(*,*) "Enter eccentricity & Lsperi" 465 write(*,*) 'Martian eccentricity (0<e<1) 425 write(*,*) 'Martian eccentricity (0<e<1)?' 466 426 call getin('excentric ',excentric) 467 427 write(*,*)"excentric =",excentric 468 write(*,*) 'Solar longitude of perihelion (0<Ls<360) 428 write(*,*) 'Solar longitude of perihelion (0<Ls<360)?' 469 429 call getin('Lsperi',Lsperi ) 470 430 write(*,*)"Lsperi=",Lsperi 471 431 Lsperi = Lsperi*pi/180.0 ! Put it in rad for peri_day 472 periheli = halfaxe*(1 -excentric)473 aphelie = halfaxe*(1 +excentric)432 periheli = halfaxe*(1 - excentric) 433 aphelie = halfaxe*(1 + excentric) 474 434 call call_dayperi(Lsperi,excentric,peri_day,year_day) 475 435 write(*,*) "Corresponding orbital params for GCM" … … 477 437 write(*,*) " aphelie = ",aphelie 478 438 write(*,*) "date of perihelion (sol)",peri_day 479 439 else 480 440 write(*,*) "paleomars=", paleomars 481 441 write(*,*) "Default present-day orbital parameters" … … 496 456 call getin("obliquit",obliquit) 497 457 write(*,*) " obliquit = ",obliquit 498 endif 499 500 endif !(.not. startfiles_1D ) 458 endif 459 endif !(.not. startfiles_1D ) 501 460 502 c latitude/longitude503 c------------------504 latitude(1)=0! default value for latitude505 write(*,*)'latitude (in degrees)?'506 507 508 latitude=latitude*pi/180.E+0 509 longitude=0.E+0 510 longitude=longitude*pi/180.E+0 511 512 ! some initializations (some of which have already been513 ! 514 ! 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 515 474 ! Mars possible matter with dttestphys in input and include!!! 516 475 ! Initializations below should mimick what is done in iniphysiq for 3D GCM 517 call init_interface_dyn_phys 518 call init_regular_lonlat(1,1,longitude,latitude, 519 & (/0.,0./),(/0.,0./)) 520 call init_geometry(1,longitude,latitude, 521 & (/0.,0.,0.,0./),(/0.,0.,0.,0./), 522 & cell_area) 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) 523 479 ! Ehouarn: init_vertial_layers called later (because disvert not called yet) 524 480 ! call init_vertical_layers(nlayer,preff,scaleheight, 525 481 ! & ap,bp,aps,bps,presnivs,pseudoalt) 526 call init_dimphy(1,nlayer) ! Initialize dimphy module 527 call phys_state_var_init(1,llm,nq,tname, 528 . day0,dayn,time, 529 . daysec,dttestphys, 530 . rad,g,r,cpp, 531 . nqperes,nqfils)! MVals: variables isotopes 532 call ini_fillgeom(1,latitude,longitude,(/1.0/)) 533 call conf_phys(1,llm,nq) 534 535 ! in 1D model physics are called every time step 536 ! ovverride iphysiq value that has been set by conf_phys 537 if (iphysiq/=1) then 538 write(*,*) "testphys1d: setting iphysiq=1" 539 iphysiq=1 540 endif 541 542 c Initialize albedo / soil thermal inertia 543 c ---------------------------------------- 544 if (.not. startfiles_1D) then 545 albedodat(1)=0.2 ! default value for albedodat 546 write(*,*)'Albedo of bare ground ?' 547 call getin("albedo",albedodat(1)) 548 write(*,*) " albedo = ",albedodat(1) 549 albedo(1,1)=albedodat(1) 550 551 inertiedat(1,1)=400 ! default value for inertiedat 552 write(*,*)'Soil thermal inertia (SI) ?' 553 call getin("inertia",inertiedat(1,1)) 554 write(*,*) " inertia = ",inertiedat(1,1) 555 556 ice_depth = -1 ! default value: no ice 557 CALL getin("subsurface_ice_depth",ice_depth) 558 559 z0(1)=z0_default ! default value for roughness 560 write(*,*) 'Surface roughness length z0 (m)?' 561 call getin("z0",z0(1)) 562 write(*,*) " z0 = ",z0(1) 563 564 endif !(.not. startfiles_1D ) 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 ) 565 516 566 517 ! Initialize local slope parameters (only matters if "callslope" 567 518 ! is .true. in callphys.def) 568 ! slope inclination angle (deg) 0: horizontal, 90: vertical 569 theta_sl(1)=0.0 ! default: no inclination 570 call getin("slope_inclination",theta_sl(1)) 571 ! slope orientation (deg) 572 ! 0 == Northward, 90 == Eastward, 180 == Southward, 270 == Westward 573 psi_sl(1)=0.0 ! default value 574 call getin("slope_orientation",psi_sl(1)) 575 576 ! sub-slopes parameters (assuming no sub-slopes distribution for now). 577 def_slope(1)=-90 ! minimum slope angle 578 def_slope(2)=90 ! maximum slope angle 579 subslope_dist(1,1)=1 ! fraction of subslopes in mesh 580 c 581 c for the gravity wave scheme 582 c --------------------------------- 583 zmea(1)=0.E+0 584 zstd(1)=0.E+0 585 zsig(1)=0.E+0 586 zgam(1)=0.E+0 587 zthe(1)=0.E+0 588 c 589 c for the slope wind scheme 590 c --------------------------------- 591 hmons(1)=0.E+0 592 write(*,*)'hmons is initialized to ',hmons(1) 593 summit(1)=0.E+0 594 write(*,*)'summit is initialized to ',summit(1) 595 base(1)=0.E+0 596 c 597 c Default values initializing the coefficients calculated later 598 c --------------------------------- 599 tauscaling(1)=1. ! calculated in aeropacity_mod.F 600 totcloudfrac(1)=1. ! calculated in watercloud_mod.F 601 602 c Specific initializations for "physiq" 603 c ------------------------------------- 604 c surface geopotential is not used (or useful) since in 1D 605 c everything is controled by surface pressure 606 phisfi(1)=0.E+0 607 608 c Initialization to take into account prescribed winds 609 c ------------------------------------------------------ 610 ptif=2.E+0*omeg*sinlat(1) 611 612 c geostrophic wind 613 gru=10. ! default value for gru 614 write(*,*)'zonal eastward component of the geostrophic', 615 &' wind (m/s) ?' 616 call getin("u",gru) 617 write(*,*) " u = ",gru 618 grv=0. !default value for grv 619 write(*,*)'meridional northward component of the geostrophic', 620 &' wind (m/s) ?' 621 call getin("v",grv) 622 write(*,*) " v = ",grv 623 624 c Initialize winds for first time step 625 if (.not. therestart1D) then 626 do ilayer = 1,nlayer 627 u(ilayer) = gru 628 v(ilayer) = grv 629 enddo 630 else 631 read(3,*) header, (u(ilayer), ilayer = 1,nlayer) 632 read(3,*) header, (v(ilayer), ilayer = 1,nlayer) 633 endif 634 w = 0. ! default: no vertical wind 635 636 c Initialize turbulente kinetic energy 637 DO ilevel=1,nlevel 638 q2(ilevel)=0.E+0 639 ENDDO 640 641 c CO2 ice on the surface 642 c ------------------- 643 ! get the index of co2 tracer (not known at this stage) 644 igcm_co2=0 645 do iq=1,nq 646 if (trim(tname(iq))=="co2") then 647 igcm_co2=iq 648 endif 649 enddo 650 if (igcm_co2==0) then 651 write(*,*) "testphys1d error, missing co2 tracer!" 652 stop 653 endif 654 655 if (.not. startfiles_1D) then 656 qsurf(igcm_co2)=0.E+0 ! default value for co2ice 657 write(*,*)'Initial CO2 ice on the surface (kg.m-2)' 658 call getin("co2ice",qsurf(igcm_co2)) 659 write(*,*) " co2ice = ",qsurf(igcm_co2) 660 endif !(.not. startfiles_1D ) 661 662 c 663 c emissivity 664 c ---------- 665 if (.not. startfiles_1D) then 666 emis=emissiv 667 IF (qsurf(igcm_co2).eq.1.E+0) THEN 668 emis=emisice(1) ! northern hemisphere 669 IF(latitude(1).LT.0) emis=emisice(2) ! southern hemisphere 670 ENDIF 671 endif !(.not. startfiles_1D ) 672 673 c Compute pressures and altitudes of atmospheric levels 674 c ---------------------------------------------------------------- 675 c Vertical Coordinates 676 c """""""""""""""""""" 677 hybrid=.true. 678 write(*,*)'Hybrid coordinates ?' 679 call getin("hybrid",hybrid) 680 write(*,*) " hybrid = ", hybrid 681 682 CALL disvert_noterre 683 ! now that disvert has been called, initialize module vertical_layers_mod 684 call init_vertical_layers(nlayer,preff,scaleheight, 685 & ap,bp,aps,bps,presnivs,pseudoalt) 686 687 DO ilevel=1,nlevel 688 plev(ilevel)=ap(ilevel)+psurf*bp(ilevel) 689 ENDDO 690 691 DO ilayer=1,nlayer 692 play(ilayer)=aps(ilayer)+psurf*bps(ilayer) 693 ENDDO 694 695 DO ilayer=1,nlayer 696 zlay(ilayer)=-200.E+0 *r*log(play(ilayer)/plev(1)) 697 & /g 698 ENDDO 699 700 701 c Initialize temperature profile 702 c -------------------------------------- 703 pks=psurf**rcp 704 705 c altitude in km in profile: divide zlay by 1000 706 tmp1(0)=0.E+0 707 DO ilayer=1,nlayer 708 tmp1(ilayer)=zlay(ilayer)/1000.E+0 709 ENDDO 710 711 call profile(nlayer+1,tmp1,tmp2) 712 713 if (.not. therestart1D) then 714 tsurf = tmp2(0) 715 do ilayer = 1,nlayer 716 temp(ilayer) = tmp2(ilayer) 717 enddo 718 else 719 read(3,*) header, tsurf, (temp(ilayer), ilayer = 1,nlayer) 720 close(3) 721 endif 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 722 650 723 651 ! Initialize soil properties and temperature 724 652 ! ------------------------------------------ 725 volcapa=1.e6 ! volumetric heat capacity 726 727 if (.not. startfiles_1D) then 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,:) = 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 728 695 ! Initialize depths 729 696 ! ----------------- 730 do isoil=1,nsoil 731 layer(isoil)=2.e-4*(2.**(isoil-1)) ! layer depth 732 enddo 733 734 ! Creating the new soil inertia table if there is subsurface ice : 735 IF (ice_depth.gt.0) THEN 736 iref = 1 ! ice/regolith boundary index 737 IF (ice_depth.lt.layer(1)) THEN 738 inertiedat(1,1) = sqrt( layer(1) / 739 & ( (ice_depth/inertiedat(1,1)**2) + 740 & ((layer(1)-ice_depth)/inertieice**2) ) ) 741 DO isoil=1,nsoil 742 inertiedat(1,isoil) = inertieice 743 ENDDO 744 ELSE ! searching for the ice/regolith boundary : 745 DO isoil=1,nsoil 746 IF ((ice_depth.ge.layer(isoil)).and. 747 & (ice_depth.lt.layer(isoil+1))) THEN 748 iref=isoil+1 749 EXIT 750 ENDIF 751 ENDDO 752 ! We then change the soil inertia table : 753 DO isoil=1,iref-1 754 inertiedat(1,isoil) = inertiedat(1,1) 755 ENDDO 756 ! We compute the transition in layer(iref) 757 inertiedat(1,iref) = sqrt( (layer(iref)-layer(iref-1)) / 758 & ( ((ice_depth-layer(iref-1))/inertiedat(1,1)**2) + 759 & ((layer(iref)-ice_depth)/inertieice**2) ) ) 760 ! Finally, we compute the underlying ice : 761 DO isoil=iref+1,nsoil 762 inertiedat(1,isoil) = inertieice 763 ENDDO 764 ENDIF ! (ice_depth.lt.layer(1)) 765 ELSE ! ice_depth <0 all is set to surface thermal inertia 766 DO isoil=1,nsoil 767 inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia 768 ENDDO 769 ENDIF ! ice_depth.gt.0 770 771 inertiesoil(1,:,1) = inertiedat(1,:) 772 773 DO isoil = 1,nsoil 774 tsoil(isoil)=tsurf(1) ! soil temperature 775 ENDDO 776 777 endif !(.not. startfiles_1D) 778 779 flux_geo_tmp=0. 780 call getin("flux_geo",flux_geo_tmp) 781 flux_geo(:,:) = flux_geo_tmp 782 783 ! Initialize depths 784 ! ----------------- 785 do isoil=0,nsoil-1 786 mlayer(isoil)=2.e-4*(2.**(isoil-0.5)) ! mid-layer depth 787 enddo 788 do isoil=1,nsoil 789 layer(isoil)=2.e-4*(2.**(isoil-1)) ! layer depth 790 enddo 791 792 c Initialize traceurs 793 c --------------------------- 794 if (photochem.or.callthermos) then 795 write(*,*) 'Initializing chemical species' 796 ! flagthermo=0: initialize over all atmospheric layers 797 flagthermo=0 798 ! check if "h2o_vap" has already been initialized 799 ! (it has been if there is a "profile_h2o_vap" file around) 800 inquire(file="profile_h2o_vap",exist=present) 801 if (present) then 802 flagh2o=0 ! 0: do not initialize h2o_vap 803 else 804 flagh2o=1 ! 1: initialize h2o_vap in inichim_newstart 805 endif 806 807 ! hack to accomodate that inichim_newstart assumes that 808 ! q & psurf arrays are on the dynamics scalar grid 809 allocate(qdyn(2,1,llm,nq),psdyn(2,1)) 810 qdyn(1,1,1:llm,1:nq)=q(1:llm,1:nq) 811 psdyn(1:2,1)=psurf 812 call inichim_newstart(ngrid,nq,qdyn,qsurf,psdyn, 813 $ flagh2o,flagthermo) 814 q(1:llm,1:nq)=qdyn(1,1,1:llm,1:nq) 815 endif 816 817 c Check if the surface is a water ice reservoir 818 c -------------------------------------------------- 819 if (.not. startfiles_1D) watercap(1,:)=0 ! Initialize watercap 820 watercaptag(1)=.false. ! Default: no water ice reservoir 821 write(*,*)'Water ice cap on ground ?' 822 call getin("watercaptag",watercaptag) 823 write(*,*) " watercaptag = ",watercaptag 824 825 c Check if the atmospheric water profile is specified 826 c --------------------------------------------------- 827 ! Adding an option to force atmospheric water values JN 828 atm_wat_profile = -1. ! Default: free atm wat profile 829 if (water) then 830 write(*,*)'Force atmospheric water vapor profile?' 831 call getin('atm_wat_profile',atm_wat_profile) 832 write(*,*) 'atm_wat_profile = ', atm_wat_profile 833 if (abs(atm_wat_profile + 1.) < 1.e-15) then ! if == -1. 834 write(*,*) 'Free atmospheric water vapor profile' 835 write(*,*) 'Total water is conserved in the column' 836 else if (abs(atm_wat_profile) < 1.e-15) then ! if == 0. 837 write(*,*) 'Dry atmospheric water vapor profile' 838 else if (0. < atm_wat_profile .and. atm_wat_profile <= 1.) 839 & then 840 write(*,*) 'Prescribed atmospheric water vapor profile' 841 write(*,*) 'Unless it reaches saturation (maximal value)' 842 else 843 write(*,*) 'Water vapor profile value not correct!' 844 stop 845 endif 846 endif 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 847 755 848 756 ! Check if the atmospheric water profile relaxation is specified 849 757 ! -------------------------------------------------------------- 850 ! Adding an option to relax atmospheric water values JBC 851 atm_wat_tau = -1. ! Default: no time relaxation 852 if (water) then 853 write(*,*) 'Relax atmospheric water vapor profile?' 854 call getin('atm_wat_tau',atm_wat_tau) 855 write(*,*) 'atm_wat_tau = ', atm_wat_tau 856 if (atm_wat_tau < 0.) then 857 write(*,*) 'Atmospheric water vapor profile is not', 858 &' relaxed.' 859 else 860 if (0. <= atm_wat_profile .and. atm_wat_profile <= 1.) 861 & then 862 write(*,*) 'Relaxed atmospheric water vapor profile', 863 &' towards ', atm_wat_profile 864 write(*,*) 'Unless it reaches saturation (maximal', 865 &' value)' 866 else 867 write(*,*) 'Reference atmospheric water vapor', 868 &' profile not known!' 869 write(*,*) 'Please, specify atm_wat_profile' 870 stop 871 endif 872 endif 873 endif 874 875 c Write a "startfi" file 876 c -------------------- 877 c This file will be read during the first call to "physiq". 878 c It is needed to transfert physics variables to "physiq"... 879 880 if (.not. startfiles_1D) then 881 call physdem0("startfi.nc",longitude,latitude,nsoilmx,ngrid, 882 & llm,nq,dttestphys,float(day0),0.,cell_area, 883 & albedodat,inertiedat,def_slope,subslope_dist) 884 call physdem1("startfi.nc",nsoilmx,ngrid,llm,nq, 885 & dttestphys,time, 886 & tsurf,tsoil,inertiesoil,albedo,emis, 887 & q2,qsurf,tauscaling, 888 & totcloudfrac,wstar,watercap,perenial_co2ice) 889 endif !(.not. startfiles_1D ) 890 891 c======================================================================= 892 c 1D MODEL TIME STEPPING LOOP 893 c======================================================================= 894 c 895 firstcall=.true. 896 lastcall=.false. 897 DO idt=1,ndt 898 IF (idt.eq.ndt) lastcall=.true. 899 c IF (idt.eq.ndt-day_step-1) then !test 900 c lastcall=.true. 901 c call solarlong(day*1.0,zls) 902 c write(103,*) 'Ls=',zls*180./pi 903 c write(103,*) 'Lat=', latitude(1)*180./pi 904 c write(103,*) 'Tau=', tauvis/odpref*psurf 905 c write(103,*) 'RunEnd - Atmos. Temp. File' 906 c write(103,*) 'RunEnd - Atmos. Temp. File' 907 c write(104,*) 'Ls=',zls*180./pi 908 c write(104,*) 'Lat=', latitude(1) 909 c write(104,*) 'Tau=', tauvis/odpref*psurf 910 c write(104,*) 'RunEnd - Atmos. Temp. File' 911 c ENDIF 912 913 c compute geopotential 914 c ~~~~~~~~~~~~~~~~~~~~~ 915 DO ilayer=1,nlayer 916 s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp 917 h(ilayer)=cpp*temp(ilayer)/(pks*s(ilayer)) 918 ENDDO 919 phi(1)=pks*h(1)*(1.E+0-s(1)) 920 DO ilayer=2,nlayer 921 phi(ilayer)=phi(ilayer-1)+ 922 & pks*(h(ilayer-1)+h(ilayer))*.5E+0 923 & *(s(ilayer-1)-s(ilayer)) 924 925 ENDDO 926 927 ! Modify atmospheric water if asked 928 ! Added "atm_wat_profile" flag (JN + JBC) 929 ! --------------------------------------- 930 if (water) then 931 call watersat(nlayer,temp,play,zqsat) 932 if (0. <= atm_wat_profile .and. atm_wat_profile <= 1.) then 933 ! If atmospheric water is monitored 934 if (atm_wat_tau < 0.) then ! Prescribed atm_wat_profile: wet if >0, dry if =0 935 q(:,igcm_h2o_vap) = min(zqsat(:),atm_wat_profile*g/psurf) 936 q(:,igcm_h2o_ice) = 0. ! reset h2o ice 937 else ! Relaxation towards the value atm_wat_profile with relaxation time atm_wat_tau 938 q(:,igcm_h2o_vap) = atm_wat_profile*g/psurf + (q(:,igcm_h2o_vap) 939 & - atm_wat_profile*g/psurf)*dexp(-dttestphys/atm_wat_tau) 940 q(:,igcm_h2o_vap) = min(zqsat(:),q(:,igcm_h2o_vap)) 941 q(:,igcm_h2o_ice) = 0. ! reset h2o ice 942 endif 943 endif 944 endif 945 946 ! write(*,*) "testphys1d avant q", q(1,:) 947 c call physics 948 c -------------------- 949 CALL physiq (1,llm,nq,firstcall,lastcall, 950 . day,time,dttestphys,plev,play,phi, 951 . u,v,temp,q,w, 952 C - outputs 953 . du, dv, dtemp, dq,dpsurf) 954 ! write(*,*) "testphys1d apres q", q(1,:) 955 956 c wind increment : specific for 1D 957 c -------------------------------- 958 c The physics compute the tendencies on u and v, 959 c here we just add Coriolos effect 960 c 961 c DO ilayer=1,nlayer 962 c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv) 963 c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru) 964 c ENDDO 965 966 c For some tests : No coriolis force at equator 967 c if(latitude(1).eq.0.) then 968 DO ilayer=1,nlayer 969 du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4 970 dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4 971 ENDDO 972 c end if 973 974 c Compute time for next time step 975 c --------------------------------------- 976 firstcall=.false. 977 time=time+dttestphys/daysec 978 IF (time.gt.1.E+0) then 979 time=time-1.E+0 980 day=day+1 981 ENDIF 982 983 c compute winds and temperature for next time step 984 c ---------------------------------------------------------- 985 986 DO ilayer=1,nlayer 987 u(ilayer)=u(ilayer)+dttestphys*du(ilayer) 988 v(ilayer)=v(ilayer)+dttestphys*dv(ilayer) 989 temp(ilayer)=temp(ilayer)+dttestphys*dtemp(ilayer) 990 ENDDO 991 992 c compute pressure for next time step 993 c ---------------------------------------------------------- 994 psurf=psurf+dttestphys*dpsurf(1) ! surface pressure change 995 DO ilevel=1,nlevel 996 plev(ilevel)=ap(ilevel)+psurf*bp(ilevel) 997 ENDDO 998 DO ilayer=1,nlayer 999 play(ilayer)=aps(ilayer)+psurf*bps(ilayer) 1000 ENDDO 1001 1002 ! increment tracers 1003 DO iq = 1,nq 1004 DO ilayer=1,nlayer 1005 q(ilayer,iq)=q(ilayer,iq)+dttestphys*dq(ilayer,iq) 1006 ENDDO 1007 ENDDO 1008 ENDDO ! of idt=1,ndt ! end of time stepping loop 1009 1010 ! Writing the "restart1D.txt" file for the next run 1011 if (startfiles_1D) call writerestart1D(psurf,tsurf,nlayer,temp, 1012 &u,v,nq,noms,qsurf,q) 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 777 778 ! Write a "startfi" file 779 ! ---------------------- 780 ! This file will be read during the first call to "physiq". 781 ! It is needed to transfert physics variables to "physiq"... 782 if (.not. startfiles_1D) then 783 call physdem0("startfi.nc",longitude,latitude,nsoilmx,ngrid, & 784 llm,nq,dttestphys,float(day0),0.,cell_area, & 785 albedodat,inertiedat,def_slope,subslope_dist) 786 call physdem1("startfi.nc",nsoilmx,ngrid,llm,nq,dttestphys,time, & 787 tsurf,tsoil,inertiesoil,albedo,emis,q2,qsurf,tauscaling, & 788 totcloudfrac,wstar,watercap,perenial_co2ice) 789 endif !(.not. startfiles_1D ) 790 791 !======================================================================= 792 ! 1D MODEL TIME STEPPING LOOP 793 !======================================================================= 794 firstcall = .true. 795 lastcall = .false. 796 do idt = 1,ndt 797 if (idt == ndt) lastcall = .true. 798 ! if (idt == ndt - day_step - 1) then ! test 799 ! lastcall = .true. 800 ! call solarlong(day*1.0,zls) 801 ! write(103,*) 'Ls=',zls*180./pi 802 ! write(103,*) 'Lat=', latitude(1)*180./pi 803 ! write(103,*) 'Tau=', tauvis/odpref*psurf 804 ! write(103,*) 'RunEnd - Atmos. Temp. File' 805 ! write(103,*) 'RunEnd - Atmos. Temp. File' 806 ! write(104,*) 'Ls=',zls*180./pi 807 ! write(104,*) 'Lat=', latitude(1) 808 ! write(104,*) 'Tau=', tauvis/odpref*psurf 809 ! write(104,*) 'RunEnd - Atmos. Temp. File' 810 ! endif 811 812 ! Compute geopotential 813 ! ~~~~~~~~~~~~~~~~~~~~ 814 s(:) = (aps(:)/psurf + bps(:))**rcp 815 h(:) = cpp*temp(:)/(pks*s(:)) 816 817 phi(1) = pks*h(1)*(1. - s(1)) 818 do ilayer = 2,nlayer 819 phi(ilayer) = phi(ilayer - 1) + pks*(h(ilayer - 1) + h(ilayer))*.5*(s(ilayer - 1) - s(ilayer)) 820 enddo 821 822 ! Modify atmospheric water if asked 823 ! Added "atm_wat_profile" flag (JN + JBC) 824 ! --------------------------------------- 825 if (water) then 826 call watersat(nlayer,temp,play,zqsat) 827 if (0. <= atm_wat_profile .and. atm_wat_profile <= 1.) then 828 ! If atmospheric water is monitored 829 if (atm_wat_tau < 0.) then ! Prescribed atm_wat_profile: wet if >0, dry if =0 830 q(:,igcm_h2o_vap) = min(zqsat(:),atm_wat_profile*g/psurf) 831 q(:,igcm_h2o_ice) = 0. ! reset h2o ice 832 else ! Relaxation towards the value atm_wat_profile with relaxation time atm_wat_tau 833 q(:,igcm_h2o_vap) = atm_wat_profile*g/psurf + (q(:,igcm_h2o_vap) - atm_wat_profile*g/psurf)*dexp(-dttestphys/atm_wat_tau) 834 q(:,igcm_h2o_vap) = min(zqsat(:),q(:,igcm_h2o_vap)) 835 q(:,igcm_h2o_ice) = 0. ! reset h2o ice 836 endif 837 endif 838 endif 839 840 ! Call physics 841 ! -------------------- 842 call physiq(1,llm,nq,firstcall,lastcall,day,time,dttestphys,plev,play,phi,u,v,temp,q,w,du,dv,dtemp,dq,dpsurf) 843 ! ^----- outputs -----^ 844 845 ! Wind increment: specific for 1D 846 ! -------------------------------- 847 ! The physics compute the tendencies on u and v, 848 ! here we just add Coriolos effect 849 !do ilayer = 1,nlayer 850 ! du(ilayer) = du(ilayer) + ptif*(v(ilayer) - grv) 851 ! dv(ilayer) = dv(ilayer) + ptif*(-u(ilayer) + gru) 852 !enddo 853 ! For some tests: No coriolis force at equator 854 !if(latitude(1) == 0.) then 855 du(:) = du(:) + (gru - u(:))/1.e4 856 dv(:) = dv(:) + (grv - v(:))/1.e4 857 !endif 858 859 ! Compute time for next time step 860 ! ------------------------------- 861 firstcall = .false. 862 time = time + dttestphys/daysec 863 if (time > 1.) then 864 time = time - 1. 865 day = day + 1 866 endif 867 868 ! Compute winds and temperature for next time step 869 ! ------------------------------------------------ 870 u(:) = u(:) + dttestphys*du(:) 871 v(:) = v(:) + dttestphys*dv(:) 872 temp(:) = temp(:) + dttestphys*dtemp(:) 873 874 ! Compute pressure for next time step 875 ! ----------------------------------- 876 psurf = psurf + dttestphys*dpsurf(1) ! surface pressure change 877 plev(:) = ap(:) + psurf*bp(:) 878 play(:) = aps(:) + psurf*bps(:) 879 880 ! Increment tracers 881 q(:,:) = q(:,:) + dttestphys*dq(:,:) 882 enddo ! End of time stepping loop (idt=1,ndt) 883 884 ! Writing the "restart1D.txt" file for the next run 885 if (startfiles_1D) call writerestart1D(psurf,tsurf,nlayer,temp,u,v,nq,noms,qsurf,q) 886 887 write(*,*) "testphys1d: everything is cool." 888 889 END PROGRAM testphys1d 1013 890 1014 1015 write(*,*) "testphys1d: Everything is cool." 1016 1017 END 891 !*********************************************************************** 892 !*********************************************************************** 893 ! Dummy subroutine used only in 3D, but required to 894 ! compile testphys1d (to cleanly use writediagfi) 895 SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn) 896 897 implicit none 898 899 integer :: im, jm, ngrid, nfield 900 real, dimension(im,jm,nfield) :: pdyn 901 real, dimension(ngrid,nfield) :: pfi 902 903 if (ngrid /= 1) then 904 write(*,*) "gr_fi_dyn error: in 1D ngrid should be 1!!!" 905 stop 906 endif 907 908 pdyn(1,1,1:nfield) = pfi(1,1:nfield) 909 910 END SUBROUTINE gr_fi_dyn 1018 911 1019 c*********************************************************************** 1020 c*********************************************************************** 1021 c Dummy subroutines used only in 3D, but required to 1022 c compile testphys1d (to cleanly use writediagfi) 1023 1024 subroutine gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn) 1025 1026 IMPLICIT NONE 1027 1028 INTEGER im,jm,ngrid,nfield 1029 REAL pdyn(im,jm,nfield) 1030 REAL pfi(ngrid,nfield) 1031 1032 if (ngrid.ne.1) then 1033 write(*,*) "gr_fi_dyn error: in 1D ngrid should be 1!!!" 1034 stop 1035 endif 1036 1037 pdyn(1,1,1:nfield)=pfi(1,1:nfield) 1038 1039 end 1040 1041 c*********************************************************************** 1042 c*********************************************************************** 1043 912 !*********************************************************************** 913 !*********************************************************************** 914
Note: See TracChangeset
for help on using the changeset viewer.
