Changeset 2583 for trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90
- Timestamp:
- Nov 10, 2021, 8:46:21 AM (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90
r2542 r2583 23 23 use callkeys_mod 24 24 use ioipsl_getin_p_mod, only : getin_p 25 use mod_phys_lmdz_para, only : is_parallel 25 use mod_phys_lmdz_para, only : is_parallel, is_master, bcast 26 26 27 27 !======================================================================= … … 71 71 REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid) 72 72 integer,intent(in) :: day_ini 73 INTEGER ig,ierr 73 74 INTEGER ig 75 CHARACTER(len=20) :: rname="inifis" ! routine name, for messages 74 76 75 77 EXTERNAL iniorbit,orbite … … 109 111 ! -------------------------------------------------------------- 110 112 111 ! check that 'callphys.def' file is around 112 OPEN(99,file='callphys.def',status='old',form='formatted',iostat=ierr) 113 CLOSE(99) 114 IF(ierr.EQ.0) iscallphys=.true. !iscallphys initialised as false in callkeys_mod module 115 116 !!! IF(ierr.EQ.0) THEN 113 IF (is_master) THEN 114 ! check if 'callphys.def' file is around 115 inquire(file="callphys.def",exist=iscallphys) 116 write(*,*) trim(rname)//": iscallphys=",iscallphys 117 ENDIF 118 call bcast(iscallphys) 119 117 120 IF(iscallphys) THEN 118 PRINT* 119 PRINT* 120 PRINT*,'--------------------------------------------' 121 PRINT*,' inifis: Parametres pour la physique (callphys.def)' 122 PRINT*,'--------------------------------------------' 123 124 write(*,*) "Directory where external input files are:" 121 122 if (is_master) then 123 write(*,*) 124 write(*,*)'--------------------------------------------' 125 write(*,*)trim(rname)//': Parameters for the physics (callphys.def)' 126 write(*,*)'--------------------------------------------' 127 endif 128 129 if (is_master) write(*,*) trim(rname)//& 130 ": Directory where external input files are:" 125 131 ! default 'datadir' is set in "datadir_mod" 126 132 call getin_p("datadir",datadir) ! default path 127 write(*,*) " datadir = ",trim(datadir) 128 129 write(*,*) "Run with or without tracer transport ?" 133 if (is_master) write(*,*) trim(rname)//": datadir = ",trim(datadir) 134 135 if (is_master) write(*,*) trim(rname)//& 136 ": Run with or without tracer transport ?" 130 137 tracer=.false. ! default value 131 138 call getin_p("tracer",tracer) 132 write(*,*) " tracer = ",tracer 133 134 write(*,*) "Run with or without atm mass update ", & 135 " due to tracer evaporation/condensation?" 139 if (is_master) write(*,*) trim(rname)//": tracer = ",tracer 140 141 if (is_master) write(*,*) trim(rname)//& 142 ": Run with or without atm mass update "//& 143 " due to tracer evaporation/condensation?" 136 144 mass_redistrib=.false. ! default value 137 145 call getin_p("mass_redistrib",mass_redistrib) 138 write(*,*) " mass_redistrib = ",mass_redistrib 139 140 write(*,*) "Diurnal cycle ?" 141 write(*,*) "(if diurnal=false, diurnal averaged solar heating)" 146 if (is_master) write(*,*) trim(rname)//& 147 ": mass_redistrib = ",mass_redistrib 148 149 if (is_master) then 150 write(*,*) trim(rname)//": Diurnal cycle ?" 151 write(*,*) trim(rname)//& 152 ": (if diurnal=false, diurnal averaged solar heating)" 153 endif 142 154 diurnal=.true. ! default value 143 155 call getin_p("diurnal",diurnal) 144 write(*,*) " diurnal = ",diurnal 145 146 write(*,*) "Seasonal cycle ?" 147 write(*,*) "(if season=false, Ls stays constant, to value ", & 148 "set in 'start'" 156 if (is_master) write(*,*) trim(rname)//": diurnal = ",diurnal 157 158 if (is_master) then 159 write(*,*) trim(rname)//": Seasonal cycle ?" 160 write(*,*) trim(rname)//& 161 ": (if season=false, Ls stays constant, to value ", & 162 "set in 'start'" 163 endif 149 164 season=.true. ! default value 150 165 call getin_p("season",season) 151 write(*,*) "season = ",season166 if (is_master) write(*,*) trim(rname)//": season = ",season 152 167 153 write(*,*) "No seasonal cycle: initial day to lock the run during restart" 168 if (is_master) write(*,*) trim(rname)//& 169 ": No seasonal cycle: initial day to lock the run during restart" 154 170 noseason_day=0.0 ! default value 155 171 call getin_p("noseason_day",noseason_day) 156 write(*,*) "noseason_day=", noseason_day157 158 write(*,*) "Tidally resonant rotation ?"172 if (is_master) write(*,*) trim(rname)//": noseason_day=", noseason_day 173 174 if (is_master) write(*,*) trim(rname)//": Tidally resonant rotation ?" 159 175 tlocked=.false. ! default value 160 176 call getin_p("tlocked",tlocked) 161 write(*,*) "tlocked = ",tlocked162 163 write(*,*) "Saturn ring shadowing ?"177 if (is_master) write(*,*) trim(rname)//": tlocked = ",tlocked 178 179 if (is_master) write(*,*) trim(rname)//": Saturn ring shadowing ?" 164 180 rings_shadow = .false. 165 181 call getin_p("rings_shadow", rings_shadow) 166 write(*,*) "rings_shadow = ", rings_shadow182 if (is_master) write(*,*) trim(rname)//": rings_shadow = ", rings_shadow 167 183 168 write(*,*) "Compute latitude-dependent gravity field?" 184 if (is_master) write(*,*) trim(rname)//& 185 ": Compute latitude-dependent gravity field?" 169 186 oblate = .false. 170 187 call getin_p("oblate", oblate) 171 write(*,*) "oblate = ", oblate172 173 write(*,*) "Flattening of the planet (a-b)/a "188 if (is_master) write(*,*) trim(rname)//": oblate = ", oblate 189 190 if (is_master) write(*,*) trim(rname)//": Flattening of the planet (a-b)/a " 174 191 flatten = 0.0 175 192 call getin_p("flatten", flatten) 176 write(*,*) "flatten = ", flatten193 if (is_master) write(*,*) trim(rname)//": flatten = ", flatten 177 194 178 195 179 write(*,*) "Needed if oblate=.true.: J2"196 if (is_master) write(*,*) trim(rname)//": Needed if oblate=.true.: J2" 180 197 J2 = 0.0 181 198 call getin_p("J2", J2) 182 write(*,*) "J2 = ", J2199 if (is_master) write(*,*) trim(rname)//": J2 = ", J2 183 200 184 write(*,*) "Needed if oblate=.true.: Planet mass (*1e24 kg)" 201 if (is_master) write(*,*) trim(rname)//& 202 ": Needed if oblate=.true.: Planet mass (*1e24 kg)" 185 203 MassPlanet = 0.0 186 204 call getin_p("MassPlanet", MassPlanet) 187 write(*,*) "MassPlanet = ", MassPlanet 188 189 write(*,*) "Needed if oblate=.true.: Planet mean radius (m)" 205 if (is_master) write(*,*) trim(rname)//": MassPlanet = ", MassPlanet 206 207 if (is_master) write(*,*) trim(rname)//& 208 ": Needed if oblate=.true.: Planet mean radius (m)" 190 209 Rmean = 0.0 191 210 call getin_p("Rmean", Rmean) 192 write(*,*) "Rmean = ", Rmean211 if (is_master) write(*,*) trim(rname)//": Rmean = ", Rmean 193 212 194 213 ! Test of incompatibility: 195 214 ! if tlocked, then diurnal should be false 196 215 if (tlocked.and.diurnal) then 197 print*,'If diurnal=true, we should turn off tlocked.' 198 stop 199 endif 200 201 write(*,*) "Tidal resonance ratio ?" 216 call abort_physic(rname,"If diurnal=true, we should turn off tlocked",1) 217 endif 218 219 if (is_master) write(*,*) trim(rname)//": Tidal resonance ratio ?" 202 220 nres=0 ! default value 203 221 call getin_p("nres",nres) 204 write(*,*) "nres = ",nres 205 206 write(*,*) "Write some extra output to the screen ?" 222 if (is_master) write(*,*) trim(rname)//": nres = ",nres 223 224 if (is_master) write(*,*) trim(rname)//& 225 ": Write some extra output to the screen ?" 207 226 lwrite=.false. ! default value 208 227 call getin_p("lwrite",lwrite) 209 write(*,*) " lwrite = ",lwrite 210 211 write(*,*) "Save statistics in file stats.nc ?" 228 if (is_master) write(*,*) trim(rname)//": lwrite = ",lwrite 229 230 if (is_master) write(*,*) trim(rname)//& 231 ": Save statistics in file stats.nc ?" 212 232 callstats=.true. ! default value 213 233 call getin_p("callstats",callstats) 214 write(*,*) " callstats = ",callstats 215 216 write(*,*) "Test energy conservation of model physics ?" 234 if (is_master) write(*,*) trim(rname)//": callstats = ",callstats 235 236 if (is_master) write(*,*) trim(rname)//& 237 ": Test energy conservation of model physics ?" 217 238 enertest=.false. ! default value 218 239 call getin_p("enertest",enertest) 219 write(*,*) " enertest = ",enertest 220 221 write(*,*) "Check to see if cpp values used match gases.def ?" 240 if (is_master) write(*,*) trim(rname)//": enertest = ",enertest 241 242 if (is_master) write(*,*) trim(rname)//& 243 ": Check to see if cpp values used match gases.def ?" 222 244 check_cpp_match=.true. ! default value 223 245 call getin_p("check_cpp_match",check_cpp_match) 224 write(*,*) " check_cpp_match = ",check_cpp_match 225 226 write(*,*) "call radiative transfer ?" 246 if (is_master) write(*,*) trim(rname)//& 247 ": check_cpp_match = ",check_cpp_match 248 249 if (is_master) write(*,*) trim(rname)//": call radiative transfer ?" 227 250 callrad=.true. ! default value 228 251 call getin_p("callrad",callrad) 229 write(*,*) " callrad = ",callrad 230 231 write(*,*) "call correlated-k radiative transfer ?" 252 if (is_master) write(*,*) trim(rname)//": callrad = ",callrad 253 254 if (is_master) write(*,*) trim(rname)//& 255 ": call correlated-k radiative transfer ?" 232 256 corrk=.true. ! default value 233 257 call getin_p("corrk",corrk) 234 write(*,*) " corrk = ",corrk 235 236 write(*,*) "prohibit calculations outside corrk T grid?" 258 if (is_master) write(*,*) trim(rname)//": corrk = ",corrk 259 260 if (is_master) write(*,*) trim(rname)//& 261 ": prohibit calculations outside corrk T grid?" 237 262 strictboundcorrk=.true. ! default value 238 263 call getin_p("strictboundcorrk",strictboundcorrk) 239 write(*,*) "strictboundcorrk = ",strictboundcorrk 240 241 write(*,*) "Minimum atmospheric temperature for Planck function integration ?" 264 if (is_master) write(*,*) trim(rname)//& 265 ": strictboundcorrk = ",strictboundcorrk 266 267 if (is_master) write(*,*) trim(rname)//& 268 ": Minimum atmospheric temperature for Planck function integration ?" 242 269 tplanckmin=30.0 ! default value 243 270 call getin_p("tplanckmin",tplanckmin) 244 write(*,*) "tplanckmin = ",tplanckmin271 if (is_master) write(*,*) trim(rname)//": tplanckmin = ",tplanckmin 245 272 246 write(*,*) "Maximum atmospheric temperature for Planck function integration ?" 273 if (is_master) write(*,*) trim(rname)//& 274 ": Maximum atmospheric temperature for Planck function integration ?" 247 275 tplanckmax=1500.0 ! default value 248 276 call getin_p("tplanckmax",tplanckmax) 249 write(*,*) "tplanckmax = ",tplanckmax277 if (is_master) write(*,*) trim(rname)//": tplanckmax = ",tplanckmax 250 278 251 write(*,*) "Temperature step for Planck function integration ?" 279 if (is_master) write(*,*) trim(rname)//& 280 ": Temperature step for Planck function integration ?" 252 281 dtplanck=0.1 ! default value 253 282 call getin_p("dtplanck",dtplanck) 254 write(*,*) "dtplanck = ",dtplanck283 if (is_master) write(*,*) trim(rname)//": dtplanck = ",dtplanck 255 284 256 write(*,*) "call gaseous absorption in the visible bands?", & 257 "(matters only if callrad=T)" 285 if (is_master) write(*,*) trim(rname)//& 286 ": call gaseous absorption in the visible bands?"//& 287 " (matters only if callrad=T)" 258 288 callgasvis=.false. ! default value 259 289 call getin_p("callgasvis",callgasvis) 260 write(*,*) "callgasvis = ",callgasvis290 if (is_master) write(*,*) trim(rname)//": callgasvis = ",callgasvis 261 291 262 write(*,*) "call continuum opacities in radiative transfer ?", & 263 "(matters only if callrad=T)" 292 if (is_master) write(*,*) trim(rname)//& 293 ": call continuum opacities in radiative transfer ?"//& 294 " (matters only if callrad=T)" 264 295 continuum=.true. ! default value 265 296 call getin_p("continuum",continuum) 266 write(*,*) "continuum = ",continuum297 if (is_master) write(*,*) trim(rname)//": continuum = ",continuum 267 298 268 write(*,*) "version for H2H2 CIA file ?"299 if (is_master) write(*,*) trim(rname)//": version for H2H2 CIA file ?" 269 300 versH2H2cia=2011 ! default value (should be 2018 but retrocompatibility first) 270 301 call getin_p("versH2H2cia",versH2H2cia) 271 write(*,*) "versH2H2cia = ",versH2H2cia302 if (is_master) write(*,*) trim(rname)//": versH2H2cia = ",versH2H2cia 272 303 ! Sanity check 273 304 if (versH2H2cia.ne.2011 .and. versH2H2cia.ne.2018) then 274 print*,'Error: Choose a correct value (2011 or 2018) for versH2H2cia !' 275 call abort 305 call abort_physic(rname,"Error: Choose a correct value (2011 or 2018) for versH2H2cia !",1) 276 306 endif 277 307 278 write(*,*) "call turbulent vertical diffusion ?" 308 if (is_master) write(*,*) trim(rname)//& 309 ": call turbulent vertical diffusion ?" 279 310 calldifv=.true. ! default value 280 311 call getin_p("calldifv",calldifv) 281 write(*,*) "calldifv = ",calldifv282 283 write(*,*) "use turbdiff instead of vdifc ?"312 if (is_master) write(*,*) trim(rname)//": calldifv = ",calldifv 313 314 if (is_master) write(*,*) trim(rname)//": use turbdiff instead of vdifc ?" 284 315 UseTurbDiff=.true. ! default value 285 316 call getin_p("UseTurbDiff",UseTurbDiff) 286 write(*,*) "UseTurbDiff = ",UseTurbDiff287 288 write(*,*) "call convective adjustment ?"317 if (is_master) write(*,*) trim(rname)//": UseTurbDiff = ",UseTurbDiff 318 319 if (is_master) write(*,*) trim(rname)//": call convective adjustment ?" 289 320 calladj=.true. ! default value 290 321 call getin_p("calladj",calladj) 291 write(*,*) " calladj = ",calladj 292 293 write(*,*)"call Lott's non-oro GWs parameterisation scheme ?" 322 if (is_master) write(*,*) trim(rname)//": calladj = ",calladj 323 324 if (is_master) write(*,*)trim(rname)//& 325 ": call Lott's non-oro GWs parameterisation scheme ?" 294 326 calllott_nonoro=.false. ! default value 295 327 call getin_p("calllott_nonoro",calllott_nonoro) 296 write(*,*)" calllott_nonoro = ",calllott_nonoro 328 if (is_master) write(*,*)trim(rname)//& 329 ": calllott_nonoro = ",calllott_nonoro 297 330 298 write(*,*)"Max Value of EP flux at launching altitude" 331 if (is_master) write(*,*)trim(rname)//& 332 ": Max Value of EP flux at launching altitude" 299 333 epflux_max=0.0 ! default value 300 334 call getin_p("epflux_max",epflux_max) 301 write(*,*)"epflux_max = ",epflux_max302 303 write(*,*) "call thermal plume model ?"335 if (is_master) write(*,*)trim(rname)//": epflux_max = ",epflux_max 336 337 if (is_master) write(*,*) trim(rname)//": call thermal plume model ?" 304 338 calltherm=.false. ! default value 305 339 call getin_p("calltherm",calltherm) 306 write(*,*) "calltherm = ",calltherm307 308 write(*,*) "call CO2 condensation ?"340 if (is_master) write(*,*) trim(rname)//": calltherm = ",calltherm 341 342 if (is_master) write(*,*) trim(rname)//": call CO2 condensation ?" 309 343 co2cond=.false. ! default value 310 344 call getin_p("co2cond",co2cond) 311 write(*,*) "co2cond = ",co2cond345 if (is_master) write(*,*) trim(rname)//": co2cond = ",co2cond 312 346 ! Test of incompatibility 313 347 if (co2cond.and.(.not.tracer)) then 314 print*,'We need a CO2 ice tracer to condense CO2' 315 call abort 348 call abort_physic(rname,"Error: We need a CO2 ice tracer to condense CO2!",1) 316 349 endif 317 350 318 write(*,*) "CO2 supersaturation level ?"351 if (is_master) write(*,*) trim(rname)//": CO2 supersaturation level ?" 319 352 co2supsat=1.0 ! default value 320 353 call getin_p("co2supsat",co2supsat) 321 write(*,*) " co2supsat = ",co2supsat 322 323 write(*,*) "Radiative timescale for Newtonian cooling ?" 354 if (is_master) write(*,*) trim(rname)//": co2supsat = ",co2supsat 355 356 if (is_master) write(*,*) trim(rname)//& 357 ": Radiative timescale for Newtonian cooling (in Earth days)?" 324 358 tau_relax=30. ! default value 325 359 call getin_p("tau_relax",tau_relax) 326 write(*,*) "tau_relax = ",tau_relax360 if (is_master) write(*,*) trim(rname)//": tau_relax = ",tau_relax 327 361 tau_relax=tau_relax*24*3600 ! convert Earth days --> seconds 328 362 329 write(*,*)"call thermal conduction in the soil ?" 363 if (is_master) write(*,*)trim(rname)//& 364 ": call thermal conduction in the soil ?" 330 365 callsoil=.true. ! default value 331 366 call getin_p("callsoil",callsoil) 332 write(*,*) "callsoil = ",callsoil367 if (is_master) write(*,*) trim(rname)//": callsoil = ",callsoil 333 368 334 write(*,*)"Rad transfer is computed every iradia", & 335 " physical timestep" 369 if (is_master) write(*,*)trim(rname)//& 370 ": Rad transfer is computed every iradia", & 371 " physical timestep" 336 372 iradia=1 ! default value 337 373 call getin_p("iradia",iradia) 338 write(*,*)"iradia = ",iradia374 if (is_master) write(*,*)trim(rname)//": iradia = ",iradia 339 375 340 write(*,*)"Rayleigh scattering ?"376 if (is_master) write(*,*)trim(rname)//": Rayleigh scattering ?" 341 377 rayleigh=.false. 342 378 call getin_p("rayleigh",rayleigh) 343 write(*,*)" rayleigh = ",rayleigh 344 345 write(*,*) "Use blackbody for stellar spectrum ?" 379 if (is_master) write(*,*)trim(rname)//": rayleigh = ",rayleigh 380 381 if (is_master) write(*,*) trim(rname)//& 382 ": Use blackbody for stellar spectrum ?" 346 383 stelbbody=.false. ! default value 347 384 call getin_p("stelbbody",stelbbody) 348 write(*,*) "stelbbody = ",stelbbody349 350 write(*,*) "Stellar blackbody temperature ?"385 if (is_master) write(*,*) trim(rname)//": stelbbody = ",stelbbody 386 387 if (is_master) write(*,*) trim(rname)//": Stellar blackbody temperature ?" 351 388 stelTbb=5800.0 ! default value 352 389 call getin_p("stelTbb",stelTbb) 353 write(*,*) "stelTbb = ",stelTbb354 355 write(*,*)"Output mean OLR in 1D?"390 if (is_master) write(*,*) trim(rname)//": stelTbb = ",stelTbb 391 392 if (is_master) write(*,*) trim(rname)//": Output mean OLR in 1D?" 356 393 meanOLR=.false. 357 394 call getin_p("meanOLR",meanOLR) 358 write(*,*)"meanOLR = ",meanOLR359 360 write(*,*)"Output spectral OLR in 3D?"395 if (is_master) write(*,*) trim(rname)//": meanOLR = ",meanOLR 396 397 if (is_master) write(*,*)trim(rname)//": Output spectral OLR in 3D?" 361 398 specOLR=.false. 362 399 call getin_p("specOLR",specOLR) 363 write(*,*)" specOLR = ",specOLR 364 365 write(*,*)"Output diagnostic optical thickness attenuation exp(-tau)?" 400 if (is_master) write(*,*)trim(rname)//": specOLR = ",specOLR 401 402 if (is_master) write(*,*)trim(rname)//& 403 ": Output diagnostic optical thickness attenuation exp(-tau)?" 366 404 diagdtau=.false. 367 405 call getin_p("diagdtau",diagdtau) 368 write(*,*)"diagdtau = ",diagdtau369 370 write(*,*)"Operate in kastprof mode?"406 if (is_master) write(*,*)trim(rname)//": diagdtau = ",diagdtau 407 408 if (is_master) write(*,*)trim(rname)//": Operate in kastprof mode?" 371 409 kastprof=.false. 372 410 call getin_p("kastprof",kastprof) 373 write(*,*)" kastprof = ",kastprof 374 375 write(*,*)"Uniform absorption in radiative transfer?" 411 if (is_master) write(*,*)trim(rname)//": kastprof = ",kastprof 412 413 if (is_master) write(*,*)trim(rname)//& 414 ": Uniform absorption in radiative transfer?" 376 415 graybody=.false. 377 416 call getin_p("graybody",graybody) 378 write(*,*)"graybody = ",graybody417 if (is_master) write(*,*)trim(rname)//": graybody = ",graybody 379 418 380 419 ! Soil model 381 write(*,*)"Number of sub-surface layers for soil scheme?" 420 if (is_master) write(*,*)trim(rname)//& 421 ": Number of sub-surface layers for soil scheme?" 382 422 ! default value of nsoilmx set in comsoil_h 383 423 call getin_p("nsoilmx",nsoilmx) 384 write(*,*)"nsoilmx=",nsoilmx424 if (is_master) write(*,*)trim(rname)//": nsoilmx=",nsoilmx 385 425 386 write(*,*)"Thickness of topmost soil layer (m)?" 426 if (is_master) write(*,*)trim(rname)//& 427 ": Thickness of topmost soil layer (m)?" 387 428 ! default value of lay1_soil set in comsoil_h 388 429 call getin_p("lay1_soil",lay1_soil) 389 write(*,*)"lay1_soil = ",lay1_soil430 if (is_master) write(*,*)trim(rname)//": lay1_soil = ",lay1_soil 390 431 391 write(*,*)"Coefficient for soil layer thickness distribution?" 432 if (is_master) write(*,*)trim(rname)//& 433 ": Coefficient for soil layer thickness distribution?" 392 434 ! default value of alpha_soil set in comsoil_h 393 435 call getin_p("alpha_soil",alpha_soil) 394 write(*,*)"alpha_soil = ",alpha_soil436 if (is_master) write(*,*)trim(rname)//": alpha_soil = ",alpha_soil 395 437 396 438 ! Slab Ocean 397 write(*,*) "Use slab-ocean ?"439 if (is_master) write(*,*) trim(rname)//": Use slab-ocean ?" 398 440 ok_slab_ocean=.false. ! default value 399 441 call getin_p("ok_slab_ocean",ok_slab_ocean) 400 write(*,*) "ok_slab_ocean = ",ok_slab_ocean442 if (is_master) write(*,*) trim(rname)//": ok_slab_ocean = ",ok_slab_ocean 401 443 ! Sanity check: for now slab oncean only works in serial mode 402 444 if (ok_slab_ocean.and.is_parallel) then 403 write(*,*) " Error: slab ocean should only be used in serial mode!" 404 call abort 405 endif 406 407 write(*,*) "Use slab-sea-ice ?" 445 call abort_physic(rname," Error: slab ocean should only be used in serial mode!",1) 446 endif 447 448 if (is_master) write(*,*) trim(rname)//": Use slab-sea-ice ?" 408 449 ok_slab_sic=.true. ! default value 409 450 call getin_p("ok_slab_sic",ok_slab_sic) 410 write(*,*) "ok_slab_sic = ",ok_slab_sic 411 412 write(*,*) "Use heat transport for the ocean ?" 451 if (is_master) write(*,*) trim(rname)//": ok_slab_sic = ",ok_slab_sic 452 453 if (is_master) write(*,*) trim(rname)//& 454 ": Use heat transport for the ocean ?" 413 455 ok_slab_heat_transp=.true. ! default value 414 456 call getin_p("ok_slab_heat_transp",ok_slab_heat_transp) 415 write(*,*) "ok_slab_heat_transp = ",ok_slab_heat_transp 457 if (is_master) write(*,*) trim(rname)//& 458 ": ok_slab_heat_transp = ",ok_slab_heat_transp 416 459 417 460 ! Photochemistry and chemistry in the thermosphere 418 461 419 write(*,*) "Use photochemistry ?"462 if (is_master) write(*,*) trim(rname)//": Use photochemistry ?" 420 463 photochem=.false. ! default value 421 464 call getin_p("photochem",photochem) 422 write(*,*) "photochem = ",photochem423 424 write(*,*) "Use photolysis heat table ?"465 if (is_master) write(*,*) trim(rname)//": photochem = ",photochem 466 467 if (is_master) write(*,*) trim(rname)//": Use photolysis heat table ?" 425 468 photoheat=.false. ! default value 426 469 call getin_p("photoheat",photoheat) 427 write(*,*) "photoheat = ",photoheat 428 429 write(*,*) "Use photolysis online calculation ?" 470 if (is_master) write(*,*) trim(rname)//": photoheat = ",photoheat 471 472 if (is_master) write(*,*) trim(rname)//& 473 ": Use photolysis online calculation ?" 430 474 jonline=.false. ! default value 431 475 call getin_p("jonline",jonline) 432 write(*,*) "jonline = ",jonline433 434 write(*,*) "Use deposition ?"476 if (is_master) write(*,*) trim(rname)//": jonline = ",jonline 477 478 if (is_master) write(*,*) trim(rname)//": Use deposition ?" 435 479 depos=.false. ! default value 436 480 call getin_p("depos",depos) 437 write(*,*) "depos = ",depos438 439 write(*,*)"Production of haze ?"481 if (is_master) write(*,*) trim(rname)//": depos = ",depos 482 483 if (is_master) write(*,*)trim(rname)//": Production of haze ?" 440 484 haze=.false. ! default value 441 485 call getin_p("haze",haze) 442 write(*,*)"haze = ",haze486 if (is_master) write(*,*)trim(rname)//": haze = ",haze 443 487 444 488 ! Global1D mean and solar zenith angle … … 452 496 ! Test of incompatibility : if global1d is true, there should not be any diurnal cycle. 453 497 if (global1d.and.diurnal) then 454 call abort_physic( "inifis",'if global1d is true, diurnal must be set to false',1)498 call abort_physic(rname,'if global1d is true, diurnal must be set to false',1) 455 499 endif 456 500 … … 462 506 write(*,*) "szangle = ",szangle 463 507 endif 464 endif 508 endif ! of if (ngrid.eq.1) 465 509 466 510 ! Test of incompatibility: 467 511 ! if kastprof used, we must be in 1D 468 512 if (kastprof.and.(ngrid.gt.1)) then 469 print*,'kastprof can only be used in 1D!'470 call abort471 endif 472 473 write(*,*)"Stratospheric temperature for kastprof mode?"513 call abort_physic(rname,'kastprof can only be used in 1D!',1) 514 endif 515 516 if (is_master) write(*,*)trim(rname)//& 517 ": Stratospheric temperature for kastprof mode?" 474 518 Tstrat=167.0 475 519 call getin_p("Tstrat",Tstrat) 476 write(*,*)"Tstrat = ",Tstrat477 478 write(*,*)"Remove lower boundary?"520 if (is_master) write(*,*)trim(rname)//": Tstrat = ",Tstrat 521 522 if (is_master) write(*,*)trim(rname)//": Remove lower boundary?" 479 523 nosurf=.false. 480 524 call getin_p("nosurf",nosurf) 481 write(*,*)"nosurf = ",nosurf525 if (is_master) write(*,*)trim(rname)//": nosurf = ",nosurf 482 526 483 527 ! Tests of incompatibility: 484 528 if (nosurf.and.callsoil) then 485 print*,'nosurf not compatible with soil scheme!' 486 print*,'... got to make a choice!' 487 call abort 488 endif 489 490 write(*,*)"Add an internal heat flux?", & 529 if (is_master) then 530 write(*,*)trim(rname)//'nosurf not compatible with soil scheme!' 531 write(*,*)trim(rname)//'... got to make a choice!' 532 endif 533 call abort_physic(rname,"nosurf not compatible with soil scheme!",1) 534 endif 535 536 if (is_master) write(*,*)trim(rname)//": Add an internal heat flux?", & 491 537 "... matters only if callsoil=F" 492 538 intheat=0. 493 539 call getin_p("intheat",intheat) 494 write(*,*)" intheat = ",intheat 495 496 write(*,*)"Use Newtonian cooling for radiative transfer?" 540 if (is_master) write(*,*)trim(rname)//": intheat = ",intheat 541 542 if (is_master) write(*,*)trim(rname)//& 543 ": Use Newtonian cooling for radiative transfer?" 497 544 newtonian=.false. 498 545 call getin_p("newtonian",newtonian) 499 write(*,*)"newtonian = ",newtonian546 if (is_master) write(*,*)trim(rname)//": newtonian = ",newtonian 500 547 501 548 ! Tests of incompatibility: 502 549 if (newtonian.and.corrk) then 503 print*,'newtonian not compatible with correlated-k!' 504 call abort 550 call abort_physic(rname,'newtonian not compatible with correlated-k!',1) 505 551 endif 506 552 if (newtonian.and.calladj) then 507 print*,'newtonian not compatible with adjustment!' 508 call abort 553 call abort_physic(rname,'newtonian not compatible with adjustment!',1) 509 554 endif 510 555 if (newtonian.and.calldifv) then 511 print*,'newtonian not compatible with a boundary layer!' 512 call abort 513 endif 514 515 write(*,*)"Test physics timescale in 1D?" 556 call abort_physic(rname,'newtonian not compatible with a boundary layer!',1) 557 endif 558 559 if (is_master) write(*,*)trim(rname)//": Test physics timescale in 1D?" 516 560 testradtimes=.false. 517 561 call getin_p("testradtimes",testradtimes) 518 write(*,*)"testradtimes = ",testradtimes562 if (is_master) write(*,*)trim(rname)//": testradtimes = ",testradtimes 519 563 520 564 ! Test of incompatibility: 521 565 ! if testradtimes used, we must be in 1D 522 566 if (testradtimes.and.(ngrid.gt.1)) then 523 print*,'testradtimes can only be used in 1D!' 524 call abort 525 endif 526 527 write(*,*)"Default planetary temperature?" 567 call abort_physic(rname,'testradtimes can only be used in 1D!',1) 568 endif 569 570 if (is_master) write(*,*)trim(rname)//": Default planetary temperature?" 528 571 tplanet=215.0 529 572 call getin_p("tplanet",tplanet) 530 write(*,*)"tplanet = ",tplanet531 532 write(*,*)"Which star?"573 if (is_master) write(*,*)trim(rname)//": tplanet = ",tplanet 574 575 if (is_master) write(*,*)trim(rname)//": Which star?" 533 576 startype=1 ! default value = Sol 534 577 call getin_p("startype",startype) 535 write(*,*)"startype = ",startype536 537 write(*,*)"Value of stellar flux at 1 AU?"578 if (is_master) write(*,*)trim(rname)//": startype = ",startype 579 580 if (is_master) write(*,*)trim(rname)//": Value of stellar flux at 1 AU?" 538 581 Fat1AU=1356.0 ! default value = Sol today 539 582 call getin_p("Fat1AU",Fat1AU) 540 write(*,*)"Fat1AU = ",Fat1AU583 if (is_master) write(*,*)trim(rname)//": Fat1AU = ",Fat1AU 541 584 542 585 543 586 ! TRACERS: 544 587 545 write(*,*)"Varying H2O cloud fraction?"588 if (is_master) write(*,*)trim(rname)//": Varying H2O cloud fraction?" 546 589 CLFvarying=.false. ! default value 547 590 call getin_p("CLFvarying",CLFvarying) 548 write(*,*)" CLFvarying = ",CLFvarying 549 550 write(*,*)"Value of fixed H2O cloud fraction?" 591 if (is_master) write(*,*)trim(rname)//": CLFvarying = ",CLFvarying 592 593 if (is_master) write(*,*)trim(rname)//& 594 ": Value of fixed H2O cloud fraction?" 551 595 CLFfixval=1.0 ! default value 552 596 call getin_p("CLFfixval",CLFfixval) 553 write(*,*)"CLFfixval = ",CLFfixval554 555 write(*,*)"fixed radii for Cloud particles?"597 if (is_master) write(*,*)trim(rname)//": CLFfixval = ",CLFfixval 598 599 if (is_master) write(*,*)trim(rname)//": fixed radii for Cloud particles?" 556 600 radfixed=.false. ! default value 557 601 call getin_p("radfixed",radfixed) 558 write(*,*)"radfixed = ",radfixed602 if (is_master) write(*,*)trim(rname)//": radfixed = ",radfixed 559 603 560 604 if(kastprof)then … … 562 606 endif 563 607 564 write(*,*)"Number mixing ratio of CO2 ice particles:" 608 if (is_master) write(*,*)trim(rname)//& 609 ": Number mixing ratio of CO2 ice particles:" 565 610 Nmix_co2=1.e6 ! default value 566 611 call getin_p("Nmix_co2",Nmix_co2) 567 write(*,*)"Nmix_co2 = ",Nmix_co2612 if (is_master) write(*,*)trim(rname)//": Nmix_co2 = ",Nmix_co2 568 613 569 614 ! write(*,*)"Number of radiatively active aerosols:" … … 572 617 ! write(*,*)" naerkind = ",naerkind 573 618 574 write(*,*)"Opacity of dust (if used):"619 if (is_master) write(*,*)trim(rname)//": Opacity of dust (if used):" 575 620 dusttau=0. ! default value 576 621 call getin_p("dusttau",dusttau) 577 write(*,*)"dusttau = ",dusttau578 579 write(*,*)"Radiatively active CO2 aerosols?"622 if (is_master) write(*,*)trim(rname)//": dusttau = ",dusttau 623 624 if (is_master) write(*,*)trim(rname)//": Radiatively active CO2 aerosols?" 580 625 aeroco2=.false. ! default value 581 626 call getin_p("aeroco2",aeroco2) 582 write(*,*)"aeroco2 = ",aeroco2583 584 write(*,*)"Fixed CO2 aerosol distribution?"627 if (is_master) write(*,*)trim(rname)//": aeroco2 = ",aeroco2 628 629 if (is_master) write(*,*)trim(rname)//": Fixed CO2 aerosol distribution?" 585 630 aerofixco2=.false. ! default value 586 631 call getin_p("aerofixco2",aerofixco2) 587 write(*,*)"aerofixco2 = ",aerofixco2588 589 write(*,*)"Radiatively active water ice?"632 if (is_master) write(*,*)trim(rname)//": aerofixco2 = ",aerofixco2 633 634 if (is_master) write(*,*)trim(rname)//": Radiatively active water ice?" 590 635 aeroh2o=.false. ! default value 591 636 call getin_p("aeroh2o",aeroh2o) 592 write(*,*)"aeroh2o = ",aeroh2o593 594 write(*,*)"Fixed H2O aerosol distribution?"637 if (is_master) write(*,*)trim(rname)//": aeroh2o = ",aeroh2o 638 639 if (is_master) write(*,*)trim(rname)//": Fixed H2O aerosol distribution?" 595 640 aerofixh2o=.false. ! default value 596 641 call getin_p("aerofixh2o",aerofixh2o) 597 write(*,*)" aerofixh2o = ",aerofixh2o 598 599 write(*,*)"Radiatively active sulfuric acid aerosols?" 642 if (is_master) write(*,*)trim(rname)//": aerofixh2o = ",aerofixh2o 643 644 if (is_master) write(*,*)trim(rname)//& 645 ": Radiatively active sulfuric acid aerosols?" 600 646 aeroh2so4=.false. ! default value 601 647 call getin_p("aeroh2so4",aeroh2so4) 602 write(*,*)"aeroh2so4 = ",aeroh2so4648 if (is_master) write(*,*)trim(rname)//": aeroh2so4 = ",aeroh2so4 603 649 604 write(*,*)"Radiatively active auroral aerosols?" 650 if (is_master) write(*,*)trim(rname)//& 651 ": Radiatively active auroral aerosols?" 605 652 aeroaurora=.false. ! default value 606 653 call getin_p("aeroaurora",aeroaurora) 607 write(*,*)"aeroaurora = ",aeroaurora654 if (is_master) write(*,*)trim(rname)//": aeroaurora = ",aeroaurora 608 655 609 656 !================================= … … 611 658 ! You should now use N-LAYER scheme (see below). 612 659 613 write(*,*)"Radiatively active two-layer aerosols?" 660 if (is_master) write(*,*)trim(rname)//& 661 ": Radiatively active two-layer aerosols?" 614 662 aeroback2lay=.false. ! default value 615 663 call getin_p("aeroback2lay",aeroback2lay) 616 write(*,*)"aeroback2lay = ",aeroback2lay617 618 if (aeroback2lay ) then664 if (is_master) write(*,*)trim(rname)//": aeroback2lay = ",aeroback2lay 665 666 if (aeroback2lay.and.is_master) then 619 667 print*,'Warning : The TWOLAY AEROSOL scheme is deprecated and buggy...' 620 668 print*,'You should use the generic n-layer scheme (see aeronlay).' 621 669 endif 622 670 623 write(*,*)"Radiatively active ammonia cloud?"671 if (is_master) write(*,*)trim(rname)//": Radiatively active ammonia cloud?" 624 672 aeronh3=.false. ! default value 625 673 call getin_p("aeronh3",aeronh3) 626 write(*,*)"aeronh3 = ",aeronh3627 628 if (aeronh3 ) then674 if (is_master) write(*,*)trim(rname)//": aeronh3 = ",aeronh3 675 676 if (aeronh3.and.is_master) then 629 677 print*,'Warning : You are using specific NH3 cloud scheme ...' 630 678 print*,'You should use the generic n-layer scheme (see aeronlay).' 631 679 endif 632 680 633 write(*,*)"TWOLAY AEROSOL: total optical depth ", & 634 "in the tropospheric layer (visible)" 681 if (is_master) write(*,*)trim(rname)//& 682 ": TWOLAY AEROSOL: total optical depth "//& 683 "in the tropospheric layer (visible)" 635 684 obs_tau_col_tropo=8.D0 636 685 call getin_p("obs_tau_col_tropo",obs_tau_col_tropo) 637 write(*,*)" obs_tau_col_tropo = ",obs_tau_col_tropo 638 639 write(*,*)"TWOLAY AEROSOL: total optical depth ", & 640 "in the stratospheric layer (visible)" 686 if (is_master) write(*,*)trim(rname)//& 687 ": obs_tau_col_tropo = ",obs_tau_col_tropo 688 689 if (is_master) write(*,*)trim(rname)//& 690 ": TWOLAY AEROSOL: total optical depth "//& 691 "in the stratospheric layer (visible)" 641 692 obs_tau_col_strato=0.08D0 642 693 call getin_p("obs_tau_col_strato",obs_tau_col_strato) 643 write(*,*)" obs_tau_col_strato = ",obs_tau_col_strato 644 645 write(*,*)"TWOLAY AEROSOL: optprop_back2lay_vis?" 694 if (is_master) write(*,*)trim(rname)//& 695 ": obs_tau_col_strato = ",obs_tau_col_strato 696 697 if (is_master) write(*,*)trim(rname)//& 698 ": TWOLAY AEROSOL: optprop_back2lay_vis?" 646 699 optprop_back2lay_vis = 'optprop_saturn_vis_n20.dat' 647 700 call getin_p("optprop_back2lay_vis",optprop_back2lay_vis) 648 write(*,*)" optprop_back2lay_vis = ",optprop_back2lay_vis 649 650 write(*,*)"TWOLAY AEROSOL: optprop_back2lay_ir?" 701 if (is_master) write(*,*)trim(rname)//& 702 ": optprop_back2lay_vis = ",trim(optprop_back2lay_vis) 703 704 if (is_master) write(*,*)trim(rname)//& 705 ": TWOLAY AEROSOL: optprop_back2lay_ir?" 651 706 optprop_back2lay_ir = 'optprop_saturn_ir_n20.dat' 652 707 call getin_p("optprop_back2lay_ir",optprop_back2lay_ir) 653 write(*,*)" optprop_back2lay_ir = ",optprop_back2lay_ir 708 if (is_master) write(*,*)trim(rname)//& 709 ": optprop_back2lay_ir = ",trim(optprop_back2lay_ir) 654 710 655 write(*,*)"TWOLAY AEROSOL: pres_bottom_tropo? in pa" 711 if (is_master) write(*,*)trim(rname)//& 712 ": TWOLAY AEROSOL: pres_bottom_tropo? in pa" 656 713 pres_bottom_tropo=66000.0 657 714 call getin_p("pres_bottom_tropo",pres_bottom_tropo) 658 write(*,*)" pres_bottom_tropo = ",pres_bottom_tropo 659 660 write(*,*)"TWOLAY AEROSOL: pres_top_tropo? in pa" 715 if (is_master) write(*,*)trim(rname)//& 716 ": pres_bottom_tropo = ",pres_bottom_tropo 717 718 if (is_master) write(*,*)trim(rname)//& 719 ": TWOLAY AEROSOL: pres_top_tropo? in pa" 661 720 pres_top_tropo=18000.0 662 721 call getin_p("pres_top_tropo",pres_top_tropo) 663 write(*,*)" pres_top_tropo = ",pres_top_tropo 664 665 write(*,*)"TWOLAY AEROSOL: pres_bottom_strato? in pa" 722 if (is_master) write(*,*)trim(rname)//& 723 ": pres_top_tropo = ",pres_top_tropo 724 725 if (is_master) write(*,*)trim(rname)//& 726 ": TWOLAY AEROSOL: pres_bottom_strato? in pa" 666 727 pres_bottom_strato=2000.0 667 728 call getin_p("pres_bottom_strato",pres_bottom_strato) 668 write(*,*)" pres_bottom_strato = ",pres_bottom_strato 729 if (is_master) write(*,*)trim(rname)//& 730 ": pres_bottom_strato = ",pres_bottom_strato 669 731 670 732 ! Sanity check 671 733 if (pres_bottom_strato .gt. pres_top_tropo) then 734 if(is_master) then 672 735 print*,'Error : TWOLAY AEROSOL, Please ensure that in callphys.def' 673 736 print*,'you have pres_top_tropo > pres_bottom_strato !' 674 call abort 675 endif 676 677 write(*,*)"TWOLAY AEROSOL: pres_top_strato? in pa" 737 endif 738 call abort_physic(rname," pres_top_tropo > pres_bottom_strato!",1) 739 endif 740 741 if (is_master) write(*,*)trim(rname)//& 742 ": TWOLAY AEROSOL: pres_top_strato? in pa" 678 743 pres_top_strato=100.0 679 744 call getin_p("pres_top_strato",pres_top_strato) 680 write(*,*)" pres_top_strato = ",pres_top_strato 681 682 write(*,*)"TWOLAY AEROSOL: particle size in the ", & 683 "tropospheric layer, in meters" 745 if (is_master) write(*,*)trim(rname)//& 746 ": pres_top_strato = ",pres_top_strato 747 748 if (is_master) write(*,*)trim(rname)//& 749 ": TWOLAY AEROSOL: particle size in the ", & 750 "tropospheric layer, in meters" 684 751 size_tropo=2.e-6 685 752 call getin_p("size_tropo",size_tropo) 686 write(*,*)" size_tropo = ",size_tropo 687 688 write(*,*)"TWOLAY AEROSOL: particle size in the ", & 689 "stratospheric layer, in meters" 753 if (is_master) write(*,*)trim(rname)//": size_tropo = ",size_tropo 754 755 if (is_master) write(*,*)trim(rname)//& 756 ": TWOLAY AEROSOL: particle size in the ", & 757 "stratospheric layer, in meters" 690 758 size_strato=1.e-7 691 759 call getin_p("size_strato",size_strato) 692 write(*,*)" size_strato = ",size_strato 693 694 write(*,*)"NH3 (thin) cloud: total optical depth" 760 if (is_master) write(*,*)trim(rname)//": size_strato = ",size_strato 761 762 if (is_master) write(*,*)trim(rname)//& 763 ": NH3 (thin) cloud: total optical depth" 695 764 tau_nh3_cloud=7.D0 696 765 call getin_p("tau_nh3_cloud",tau_nh3_cloud) 697 write(*,*)"tau_nh3_cloud = ",tau_nh3_cloud698 699 write(*,*)"NH3 (thin) cloud pressure level"766 if (is_master) write(*,*)trim(rname)//": tau_nh3_cloud = ",tau_nh3_cloud 767 768 if (is_master) write(*,*)trim(rname)//": NH3 (thin) cloud pressure level" 700 769 pres_nh3_cloud=7.D0 701 770 call getin_p("pres_nh3_cloud",pres_nh3_cloud) 702 write(*,*)"pres_nh3_cloud = ",pres_nh3_cloud703 704 write(*,*)"NH3 (thin) cloud: particle sizes"771 if (is_master) write(*,*)trim(rname)//": pres_nh3_cloud = ",pres_nh3_cloud 772 773 if (is_master) write(*,*)trim(rname)//": NH3 (thin) cloud: particle sizes" 705 774 size_nh3_cloud=3.e-6 706 775 call getin_p("size_nh3_cloud",size_nh3_cloud) 707 write(*,*)"size_nh3_cloud = ",size_nh3_cloud776 if (is_master) write(*,*)trim(rname)//": size_nh3_cloud = ",size_nh3_cloud 708 777 709 778 !================================= 710 779 ! Generic N-LAYER aerosol scheme 711 780 712 write(*,*)"Radiatively active generic n-layer aerosols?" 781 if (is_master) write(*,*)trim(rname)//& 782 ": Radiatively active generic n-layer aerosols?" 713 783 aeronlay=.false. ! default value 714 784 call getin_p("aeronlay",aeronlay) 715 write(*,*)" aeronlay = ",aeronlay 716 717 write(*,*)"Number of generic aerosols layers?" 785 if (is_master) write(*,*)trim(rname)//": aeronlay = ",aeronlay 786 787 if (is_master) write(*,*)trim(rname)//& 788 ": Number of generic aerosols layers?" 718 789 nlayaero=1 ! default value 719 790 call getin_p("nlayaero",nlayaero) 720 791 ! Avoid to allocate arrays of size 0 721 792 if (aeronlay .and. nlayaero.lt.1) then 793 if (is_master) then 722 794 print*, " You are trying to set no generic aerosols..." 723 795 print*, " Set aeronlay=.false. instead ! I abort." 724 call abort 796 endif 797 call abort_physic(rname,"no generic aerosols...",1) 725 798 endif 726 799 if (.not. aeronlay) nlayaero=1 727 write(*,*)"nlayaero = ",nlayaero800 if (is_master) write(*,*)trim(rname)//": nlayaero = ",nlayaero 728 801 729 802 ! This is necessary, we just set the number of aerosol layers … … 739 812 IF(.NOT.ALLOCATED(optprop_aeronlay_vis)) ALLOCATE(optprop_aeronlay_vis(nlayaero)) 740 813 741 write(*,*)"Generic n-layer aerosols: Optical depth at reference wavelenght" 814 if (is_master) write(*,*)trim(rname)//& 815 ": Generic n-layer aerosols: Optical depth at reference wavelenght" 742 816 aeronlay_tauref=1.0E-1 743 817 call getin_p("aeronlay_tauref",aeronlay_tauref) 744 write(*,*)" aeronlay_tauref = ",aeronlay_tauref 745 746 write(*,*)"Generic n-layer aerosols: Reference wavelenght for optical depths (m)" 818 if (is_master) write(*,*)trim(rname)//& 819 ": aeronlay_tauref = ",aeronlay_tauref 820 821 if (is_master) write(*,*)trim(rname)//& 822 ": Generic n-layer aerosols: Reference wavelenght for optical depths (m)" 747 823 aeronlay_lamref=0.6E-6 748 824 call getin_p("aeronlay_lamref",aeronlay_lamref) 749 write(*,*)" aeronlay_lamref = ",aeronlay_lamref 750 751 write(*,*)"Generic n-layer aerosols: Vertical profile choice : & 752 (1) Tau btwn ptop and pbot follows atm. scale height & 753 or (2) Tau above pbot follows its own scale height" 825 if (is_master) write(*,*)trim(rname)//& 826 ": aeronlay_lamref = ",aeronlay_lamref 827 828 if (is_master) then 829 write(*,*)trim(rname)//& 830 ": Generic n-layer aerosols: Vertical profile choice : " 831 write(*,*)trim(rname)//& 832 " (1) Tau btwn ptop and pbot follows atm. scale height" 833 write(*,*)trim(rname)//& 834 " or (2) Tau above pbot follows its own scale height" 835 endif 754 836 aeronlay_choice=1 755 837 call getin_p("aeronlay_choice",aeronlay_choice) 756 write(*,*)" aeronlay_choice = ",aeronlay_choice 757 758 write(*,*)"Generic n-layer aerosols: bottom pressures (Pa)" 838 if (is_master) write(*,*)trim(rname)//& 839 ": aeronlay_choice = ",aeronlay_choice 840 841 if (is_master) write(*,*)trim(rname)//& 842 ": Generic n-layer aerosols: bottom pressures (Pa)" 759 843 aeronlay_pbot=2000.0 760 844 call getin_p("aeronlay_pbot",aeronlay_pbot) 761 write(*,*)" aeronlay_pbot = ",aeronlay_pbot 762 763 write(*,*)"Generic n-layer aerosols: (if choice=1) Top pressures (Pa) " 845 if (is_master) write(*,*)trim(rname)//": aeronlay_pbot = ",aeronlay_pbot 846 847 if (is_master) write(*,*)trim(rname)//& 848 ": Generic n-layer aerosols: (if choice=1) Top pressures (Pa) " 764 849 aeronlay_ptop=300000.0 765 850 call getin_p("aeronlay_ptop",aeronlay_ptop) 766 write(*,*)" aeronlay_ptop = ",aeronlay_ptop 767 768 write(*,*)"Generic n-layer aerosols: (if choice=2) Scale height / atm. scale height" 851 if (is_master) write(*,*)trim(rname)//": aeronlay_ptop = ",aeronlay_ptop 852 853 if (is_master) write(*,*)trim(rname)//& 854 ": Generic n-layer aerosols: (if choice=2) Scale height / atm. scale height" 769 855 aeronlay_sclhght=0.2 770 856 call getin_p("aeronlay_sclhght",aeronlay_sclhght) 771 write(*,*)" aeronlay_sclhght = ",aeronlay_sclhght 772 773 write(*,*)"Generic n-layer aerosols: particles effective radii (m)" 857 if (is_master) write(*,*)trim(rname)//& 858 ": aeronlay_sclhght = ",aeronlay_sclhght 859 860 if (is_master) write(*,*)trim(rname)//& 861 ": Generic n-layer aerosols: particles effective radii (m)" 774 862 aeronlay_size=1.e-6 775 863 call getin_p("aeronlay_size",aeronlay_size) 776 write(*,*)" aeronlay_size = ",aeronlay_size 777 778 write(*,*)"Generic n-layer aerosols: particles radii effective variance" 864 if (is_master) write(*,*)trim(rname)//": aeronlay_size = ",aeronlay_size 865 866 if (is_master) write(*,*)trim(rname)//& 867 ": Generic n-layer aerosols: particles radii effective variance" 779 868 aeronlay_nueff=0.1 780 869 call getin_p("aeronlay_nueff",aeronlay_nueff) 781 write(*,*)" aeronlay_nueff = ",aeronlay_nueff 782 783 write(*,*)"Generic n-layer aerosols: VIS optical properties file" 870 if (is_master) write(*,*)trim(rname)//": aeronlay_nueff = ",aeronlay_nueff 871 872 if (is_master) write(*,*)trim(rname)//& 873 ": Generic n-layer aerosols: VIS optical properties file" 784 874 optprop_aeronlay_vis = 'optprop_saturn_vis_n20.dat' 785 875 call getin_p("optprop_aeronlay_vis",optprop_aeronlay_vis) 786 write(*,*)" optprop_aeronlay_vis = ",optprop_aeronlay_vis 787 788 write(*,*)"Generic n-layer aerosols: IR optical properties file" 876 if (is_master) write(*,*)trim(rname)//& 877 ": optprop_aeronlay_vis = ",optprop_aeronlay_vis 878 879 if (is_master) write(*,*)trim(rname)//& 880 ": Generic n-layer aerosols: IR optical properties file" 789 881 optprop_aeronlay_ir = 'optprop_saturn_ir_n20.dat' 790 882 call getin_p("optprop_aeronlay_ir",optprop_aeronlay_ir) 791 write(*,*)" optprop_aeronlay_ir = ",optprop_aeronlay_ir 883 if (is_master) write(*,*)trim(rname)//& 884 ": optprop_aeronlay_ir = ",optprop_aeronlay_ir 792 885 793 886 794 887 !================================= 795 888 796 write(*,*)"Cloud pressure level (with kastprof only):" 889 if (is_master) write(*,*)trim(rname)//& 890 ": Cloud pressure level (with kastprof only):" 797 891 cloudlvl=0. ! default value 798 892 call getin_p("cloudlvl",cloudlvl) 799 write(*,*)" cloudlvl = ",cloudlvl 800 801 write(*,*)"Is the variable gas species radiatively active?" 893 if (is_master) write(*,*)trim(rname)//": cloudlvl = ",cloudlvl 894 895 if (is_master) write(*,*)trim(rname)//& 896 ": Is the variable gas species radiatively active?" 802 897 Tstrat=167.0 803 898 varactive=.false. 804 899 call getin_p("varactive",varactive) 805 write(*,*)" varactive = ",varactive 806 807 write(*,*)"Is the variable gas species distribution set?" 900 if (is_master) write(*,*)trim(rname)//": varactive = ",varactive 901 902 if (is_master) write(*,*)trim(rname)//& 903 ": Is the variable gas species distribution set?" 808 904 varfixed=.false. 809 905 call getin_p("varfixed",varfixed) 810 write(*,*)" varfixed = ",varfixed 811 812 write(*,*)"What is the saturation % of the variable species?" 906 if (is_master) write(*,*)trim(rname)//": varfixed = ",varfixed 907 908 if (is_master) write(*,*)trim(rname)//& 909 ": What is the saturation % of the variable species?" 813 910 satval=0.8 814 911 call getin_p("satval",satval) 815 write(*,*)"satval = ",satval912 if (is_master) write(*,*)trim(rname)//": satval = ",satval 816 913 817 914 … … 819 916 ! if varactive, then varfixed should be false 820 917 if (varactive.and.varfixed) then 821 print*,'if varactive, varfixed must be OFF!' 822 stop 823 endif 824 825 write(*,*) "Gravitationnal sedimentation ?" 918 call abort_physic(rname,'if varactive, varfixed must be OFF!',1) 919 endif 920 921 if (is_master) write(*,*)trim(rname)//": Gravitationnal sedimentation ?" 826 922 sedimentation=.false. ! default value 827 923 call getin_p("sedimentation",sedimentation) 828 write(*,*) "sedimentation = ",sedimentation829 830 write(*,*) "Compute water cycle ?"924 if (is_master) write(*,*)trim(rname)//": sedimentation = ",sedimentation 925 926 if (is_master) write(*,*)trim(rname)//": Compute water cycle ?" 831 927 water=.false. ! default value 832 928 call getin_p("water",water) 833 write(*,*) "water = ",water929 if (is_master) write(*,*)trim(rname)//": water = ",water 834 930 835 931 ! Test of incompatibility: 836 932 ! if water is true, there should be at least a tracer 837 933 if (water.and.(.not.tracer)) then 838 print*,'if water is ON, tracer must be ON too!' 839 stop 840 endif 841 842 write(*,*) "Include water condensation ?" 934 call abort_physic(rname,'if water is ON, tracer must be ON too!',1) 935 endif 936 937 if (is_master) write(*,*)trim(rname)//": Include water condensation ?" 843 938 watercond=.false. ! default value 844 939 call getin_p("watercond",watercond) 845 write(*,*) "watercond = ",watercond940 if (is_master) write(*,*)trim(rname)//": watercond = ",watercond 846 941 847 942 ! Test of incompatibility: 848 943 ! if watercond is used, then water should be used too 849 944 if (watercond.and.(.not.water)) then 850 print*,'if watercond is used, water should be used too'851 stop852 endif 853 854 write(*,*) "Include water precipitation ?"945 call abort_physic(rname,& 946 'if watercond is used, water should be used too',1) 947 endif 948 949 if (is_master) write(*,*)trim(rname)//": Include water precipitation ?" 855 950 waterrain=.false. ! default value 856 951 call getin_p("waterrain",waterrain) 857 write(*,*) "waterrain = ",waterrain858 859 write(*,*) "Include surface hydrology ?"952 if (is_master) write(*,*)trim(rname)//": waterrain = ",waterrain 953 954 if (is_master) write(*,*)trim(rname)//": Include surface hydrology ?" 860 955 hydrology=.false. ! default value 861 956 call getin_p("hydrology",hydrology) 862 write(*,*) "hydrology = ",hydrology863 864 write(*,*) "Evolve surface water sources ?"957 if (is_master) write(*,*)trim(rname)//": hydrology = ",hydrology 958 959 if (is_master) write(*,*)trim(rname)//": Evolve surface water sources ?" 865 960 sourceevol=.false. ! default value 866 961 call getin_p("sourceevol",sourceevol) 867 write(*,*) "sourceevol = ",sourceevol868 869 write(*,*) "Ice evolution timestep ?"962 if (is_master) write(*,*)trim(rname)//": sourceevol = ",sourceevol 963 964 if (is_master) write(*,*)trim(rname)//": Ice evolution timestep ?" 870 965 icetstep=100.0 ! default value 871 966 call getin_p("icetstep",icetstep) 872 write(*,*) "icetstep = ",icetstep967 if (is_master) write(*,*)trim(rname)//": icetstep = ",icetstep 873 968 874 write(*,*) "Spectral Dependant albedo ?"969 if (is_master) write(*,*)trim(rname)//": Spectral Dependant albedo ?" 875 970 albedo_spectral_mode=.false. ! default value 876 971 call getin_p("albedo_spectral_mode",albedo_spectral_mode) 877 write(*,*) " albedo_spectral_mode = ",albedo_spectral_mode 878 879 write(*,*) "Snow albedo ?" 880 write(*,*) "If albedo_spectral_mode=.true., then this " 881 write(*,*) "corresponds to the 0.5 microns snow albedo." 972 if (is_master) write(*,*)trim(rname)//& 973 ": albedo_spectral_mode = ",albedo_spectral_mode 974 975 if (is_master) then 976 write(*,*)trim(rname)//": Snow albedo ?" 977 write(*,*)trim(rname)//": If albedo_spectral_mode=.true., then this " 978 write(*,*)trim(rname)//": corresponds to the 0.5 microns snow albedo." 979 endif 882 980 albedosnow=0.5 ! default value 883 981 call getin_p("albedosnow",albedosnow) 884 write(*,*) "albedosnow = ",albedosnow982 if (is_master) write(*,*)trim(rname)//": albedosnow = ",albedosnow 885 983 886 write(*,*) "Ocean albedo ?"984 if (is_master) write(*,*)trim(rname)//": Ocean albedo ?" 887 985 alb_ocean=0.07 ! default value 888 986 call getin_p("alb_ocean",alb_ocean) 889 write(*,*) "alb_ocean = ",alb_ocean987 if (is_master) write(*,*)trim(rname)//": alb_ocean = ",alb_ocean 890 988 891 write(*,*) "CO2 ice albedo ?"989 if (is_master) write(*,*)trim(rname)//": CO2 ice albedo ?" 892 990 albedoco2ice=0.5 ! default value 893 991 call getin_p("albedoco2ice",albedoco2ice) 894 write(*,*) "albedoco2ice = ",albedoco2ice895 896 write(*,*) "Maximum ice thickness ?"992 if (is_master) write(*,*)trim(rname)//": albedoco2ice = ",albedoco2ice 993 994 if (is_master) write(*,*)trim(rname)//": Maximum ice thickness ?" 897 995 maxicethick=2.0 ! default value 898 996 call getin_p("maxicethick",maxicethick) 899 write(*,*) "maxicethick = ",maxicethick900 901 write(*,*) "Freezing point of seawater ?"997 if (is_master) write(*,*)trim(rname)//": maxicethick = ",maxicethick 998 999 if (is_master) write(*,*)trim(rname)//": Freezing point of seawater ?" 902 1000 Tsaldiff=-1.8 ! default value 903 1001 call getin_p("Tsaldiff",Tsaldiff) 904 write(*,*) "Tsaldiff = ",Tsaldiff905 906 write(*,*) "Minimum eddy mix coeff in 1D ?"1002 if (is_master) write(*,*)trim(rname)//": Tsaldiff = ",Tsaldiff 1003 1004 if (is_master) write(*,*)trim(rname)//": Minimum eddy mix coeff in 1D ?" 907 1005 kmixmin=1.0e-2 ! default value 908 1006 call getin_p("kmixmin",kmixmin) 909 write(*,*) " kmixmin = ",kmixmin 910 911 write(*,*) "Does user want to force cpp and mugaz?" 1007 if (is_master) write(*,*)trim(rname)//": kmixmin = ",kmixmin 1008 1009 if (is_master) write(*,*)trim(rname)//& 1010 ": Does user want to force cpp and mugaz?" 912 1011 force_cpp=.false. ! default value 913 1012 call getin_p("force_cpp",force_cpp) 914 write(*,*) "force_cpp = ",force_cpp1013 if (is_master) write(*,*)trim(rname)//": force_cpp = ",force_cpp 915 1014 916 1015 if (force_cpp) then 917 1016 mugaz = -99999. 918 PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?' 1017 if (is_master) write(*,*)trim(rname)//& 1018 ": MEAN MOLECULAR MASS in g mol-1 ?" 919 1019 call getin_p("mugaz",mugaz) 920 1020 IF (mugaz.eq.-99999.) THEN 921 PRINT *, "mugaz must be set if force_cpp = T" 922 STOP 1021 call abort_physic(rname,"mugaz must be set if force_cpp = T",1) 923 1022 ELSE 924 write(*,*) "mugaz=",mugaz1023 if (is_master) write(*,*)trim(rname)//": mugaz=",mugaz 925 1024 ENDIF 926 1025 cpp = -99999. 927 PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?' 1026 if (is_master) write(*,*)trim(rname)//& 1027 ": SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?" 928 1028 call getin_p("cpp",cpp) 929 1029 IF (cpp.eq.-99999.) THEN … … 931 1031 STOP 932 1032 ELSE 933 write(*,*) "cpp=",cpp1033 if (is_master) write(*,*)trim(rname)//": cpp=",cpp 934 1034 ENDIF 935 1035 endif ! of if (force_cpp) … … 937 1037 call calc_cpp_mugaz 938 1038 939 PRINT*,'--------------------------------------------' 940 PRINT* 941 PRINT* 1039 if (is_master) then 1040 PRINT*,'--------------------------------------------' 1041 PRINT* 1042 PRINT* 1043 endif 942 1044 ELSE 943 write(*,*) 944 write(*,*) 'Cannot read file callphys.def. Is it here ?' 945 stop 1045 call abort_physic(rname,'Cannot read file callphys.def. Is it here ?',1) 946 1046 ENDIF ! of IF(iscallphys) 947 1047 948 PRINT* 949 PRINT*,'inifis: daysec',daysec 950 PRINT* 951 PRINT*,'inifis: The radiative transfer is computed:' 952 PRINT*,' each ',iradia,' physical time-step' 953 PRINT*,' or each ',iradia*dtphys,' seconds' 954 PRINT* 955 1048 if (is_master) then 1049 PRINT* 1050 PRINT*,'inifis: daysec',daysec 1051 PRINT* 1052 PRINT*,'inifis: The radiative transfer is computed:' 1053 PRINT*,' each ',iradia,' physical time-step' 1054 PRINT*,' or each ',iradia*dtphys,' seconds' 1055 PRINT* 1056 endif 956 1057 957 1058 !-----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.