source: trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90 @ 2109

Last change on this file since 2109 was 2078, checked in by jvatant, 6 years ago

Cf r2076, fix a missing dummy init of mugaz causing NaNs?
--JVO

File size: 25.7 KB
RevLine 
[1524]1MODULE inifis_mod
2IMPLICIT NONE
[253]3
[1524]4CONTAINS
[253]5
[1524]6  SUBROUTINE inifis(ngrid,nlayer,nq, &
[1525]7             day_ini,pdaysec,nday,ptimestep, &
[1524]8             plat,plon,parea, &
9             prad,pg,pr,pcpp)
10
[1682]11  use init_print_control_mod, only: init_print_control
[1529]12  use radinc_h, only: ini_radinc_h, naerkind
[1524]13  use datafile_mod, only: datadir
14  use comdiurn_h, only: sinlat, coslat, sinlon, coslon
[1542]15  use comgeomfi_h, only: totarea, totarea_planet
[1538]16  use comsoil_h, only: ini_comsoil_h, nsoilmx, lay1_soil, alpha_soil
[1525]17  use time_phylmdz_mod, only: ecritphy,day_step,iphysiq, &
18                              init_time, daysec, dtphys
[1524]19  use comcstfi_mod, only: rad, cpp, g, r, rcp, &
20                          mugaz, pi, avocado
21  use planete_mod, only: nres
22  use planetwide_mod, only: planetwide_sumval
23  use callkeys_mod
[1529]24  use mod_phys_lmdz_para, only : is_parallel
[1524]25
[135]26!=======================================================================
27!
28!   purpose:
29!   -------
30!
31!   Initialisation for the physical parametrisations of the LMD
[305]32!   Generic Model.
[135]33!
34!   author: Frederic Hourdin 15 / 10 /93
35!   -------
36!   modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
37!             Ehouarn Millour (oct. 2008) tracers are now identified
38!              by their names and may not be contiguously
39!              stored in the q(:,:,:,:) array
40!             E.M. (june 2009) use getin routine to load parameters
41!
42!
43!   arguments:
44!   ----------
45!
46!   input:
47!   ------
48!
49!    ngrid                 Size of the horizontal grid.
50!                          All internal loops are performed on that grid.
51!    nlayer                Number of vertical layers.
52!    pdayref               Day of reference for the simulation
53!    pday                  Number of days counted from the North. Spring
54!                          equinoxe.
55!
56!=======================================================================
57!
58!-----------------------------------------------------------------------
59!   declarations:
60!   -------------
[1524]61  use datafile_mod, only: datadir
62  use ioipsl_getin_p_mod, only: getin_p
63  IMPLICIT NONE
[1384]64
[135]65
66
[1524]67  REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep
[1525]68  INTEGER,INTENT(IN) :: nday
[1524]69  INTEGER,INTENT(IN) :: ngrid,nlayer,nq
70  REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid)
71  integer,intent(in) :: day_ini
72  INTEGER ig,ierr
[135]73 
[1524]74  EXTERNAL iniorbit,orbite
75  EXTERNAL SSUM
76  REAL SSUM
[135]77 
[1682]78  ! Initialize flags lunout, prt_level, debug (in print_control_mod)
79  CALL init_print_control
80
[1524]81  ! initialize constants in comcstfi_mod
82  rad=prad
83  cpp=pcpp
84  g=pg
85  r=pr
86  rcp=r/cpp
[2078]87  mugaz=8.314*1000./pr ! dummy init
[1524]88  pi=2.*asin(1.)
89  avocado = 6.02214179e23   ! added by RW
[135]90
[1524]91  ! Initialize some "temporal and calendar" related variables
[1832]92#ifndef MESOSCALE
[1525]93  CALL init_time(day_ini,pdaysec,nday,ptimestep)
[1832]94#endif
[135]95
[1525]96  ! read in some parameters from "run.def" for physics,
97  ! or shared between dynamics and physics.
98  call getin_p("ecritphy",ecritphy) ! frequency of outputs in physics,
99                                    ! in dynamical steps
100  call getin_p("day_step",day_step) ! number of dynamical steps per day
101  call getin_p("iphysiq",iphysiq) ! call physics every iphysiq dyn step
[135]102
[1669]103  ! do we read a startphy.nc file? (default: .true.)
104  call getin_p("startphy_file",startphy_file)
105 
[135]106! --------------------------------------------------------------
107!  Reading the "callphys.def" file controlling some key options
108! --------------------------------------------------------------
109     
[1524]110  ! check that 'callphys.def' file is around
111  OPEN(99,file='callphys.def',status='old',form='formatted',iostat=ierr)
112  CLOSE(99)
113  IF(ierr.EQ.0) iscallphys=.true. !iscallphys initialised as false in callkeys_mod module
[135]114     
[1315]115!!!      IF(ierr.EQ.0) THEN
[1524]116  IF(iscallphys) THEN
117     PRINT*
118     PRINT*
119     PRINT*,'--------------------------------------------'
120     PRINT*,' inifis: Parametres pour la physique (callphys.def)'
121     PRINT*,'--------------------------------------------'
[135]122
[1524]123     write(*,*) "Directory where external input files are:"
124     ! default 'datadir' is set in "datadir_mod"
125     call getin_p("datadir",datadir) ! default path
126     write(*,*) " datadir = ",trim(datadir)
[374]127
[1524]128     write(*,*) "Run with or without tracer transport ?"
129     tracer=.false. ! default value
130     call getin_p("tracer",tracer)
131     write(*,*) " tracer = ",tracer
[135]132
[1524]133     write(*,*) "Run with or without atm mass update ", &
134            " due to tracer evaporation/condensation?"
135     mass_redistrib=.false. ! default value
136     call getin_p("mass_redistrib",mass_redistrib)
137     write(*,*) " mass_redistrib = ",mass_redistrib
[728]138
[1524]139     write(*,*) "Diurnal cycle ?"
140     write(*,*) "(if diurnal=false, diurnal averaged solar heating)"
141     diurnal=.true. ! default value
142     call getin_p("diurnal",diurnal)
143     write(*,*) " diurnal = ",diurnal
[135]144
[1524]145     write(*,*) "Seasonal cycle ?"
146     write(*,*) "(if season=false, Ls stays constant, to value ", &
147         "set in 'start'"
148     season=.true. ! default value
149     call getin_p("season",season)
150     write(*,*) " season = ",season
[135]151
[1524]152     write(*,*) "Tidally resonant rotation ?"
153     tlocked=.false. ! default value
154     call getin_p("tlocked",tlocked)
155     write(*,*) "tlocked = ",tlocked
[135]156
[1524]157     write(*,*) "Saturn ring shadowing ?"
158     rings_shadow = .false.
159     call getin_p("rings_shadow", rings_shadow)
160     write(*,*) "rings_shadow = ", rings_shadow
[1133]161         
[1524]162     write(*,*) "Compute latitude-dependent gravity field?"
163     oblate = .false.
164     call getin_p("oblate", oblate)
165     write(*,*) "oblate = ", oblate
[1194]166
[1524]167     write(*,*) "Flattening of the planet (a-b)/a "
168     flatten = 0.0
169     call getin_p("flatten", flatten)
170     write(*,*) "flatten = ", flatten
[1174]171         
[1194]172
[1524]173     write(*,*) "Needed if oblate=.true.: J2"
174     J2 = 0.0
175     call getin_p("J2", J2)
176     write(*,*) "J2 = ", J2
[1194]177         
[1524]178     write(*,*) "Needed if oblate=.true.: Planet mass (*1e24 kg)"
179     MassPlanet = 0.0
180     call getin_p("MassPlanet", MassPlanet)
181     write(*,*) "MassPlanet = ", MassPlanet         
[1194]182
[1524]183     write(*,*) "Needed if oblate=.true.: Planet mean radius (m)"
184     Rmean = 0.0
185     call getin_p("Rmean", Rmean)
186     write(*,*) "Rmean = ", Rmean
[1194]187         
[135]188! Test of incompatibility:
189! if tlocked, then diurnal should be false
[1524]190     if (tlocked.and.diurnal) then
191       print*,'If diurnal=true, we should turn off tlocked.'
192       stop
193     endif
[135]194
[1524]195     write(*,*) "Tidal resonance ratio ?"
196     nres=0          ! default value
197     call getin_p("nres",nres)
198     write(*,*) "nres = ",nres
[135]199
[1524]200     write(*,*) "Write some extra output to the screen ?"
201     lwrite=.false. ! default value
202     call getin_p("lwrite",lwrite)
203     write(*,*) " lwrite = ",lwrite
[135]204
[1524]205     write(*,*) "Save statistics in file stats.nc ?"
206     callstats=.true. ! default value
207     call getin_p("callstats",callstats)
208     write(*,*) " callstats = ",callstats
[135]209
[1524]210     write(*,*) "Test energy conservation of model physics ?"
211     enertest=.false. ! default value
212     call getin_p("enertest",enertest)
213     write(*,*) " enertest = ",enertest
[135]214
[1524]215     write(*,*) "Check to see if cpp values used match gases.def ?"
216     check_cpp_match=.true. ! default value
217     call getin_p("check_cpp_match",check_cpp_match)
218     write(*,*) " check_cpp_match = ",check_cpp_match
[538]219
[1524]220     write(*,*) "call radiative transfer ?"
221     callrad=.true. ! default value
222     call getin_p("callrad",callrad)
223     write(*,*) " callrad = ",callrad
[135]224
[1524]225     write(*,*) "call correlated-k radiative transfer ?"
226     corrk=.true. ! default value
227     call getin_p("corrk",corrk)
228     write(*,*) " corrk = ",corrk
[135]229
[1524]230     write(*,*) "prohibit calculations outside corrk T grid?"
231     strictboundcorrk=.true. ! default value
232     call getin_p("strictboundcorrk",strictboundcorrk)
233     write(*,*) "strictboundcorrk = ",strictboundcorrk
[1145]234
[1524]235     write(*,*) "call gaseous absorption in the visible bands?", &
236                    "(matters only if callrad=T)"
237     callgasvis=.false. ! default value
238     call getin_p("callgasvis",callgasvis)
239     write(*,*) " callgasvis = ",callgasvis
[538]240       
[1524]241     write(*,*) "call continuum opacities in radiative transfer ?", &
242                    "(matters only if callrad=T)"
243     continuum=.true. ! default value
244     call getin_p("continuum",continuum)
245     write(*,*) " continuum = ",continuum
[538]246
[1524]247     write(*,*) "use analytic function for H2O continuum ?"
248     H2Ocont_simple=.false. ! default value
249     call getin_p("H2Ocont_simple",H2Ocont_simple)
250     write(*,*) " H2Ocont_simple = ",H2Ocont_simple
[538]251 
[1524]252     write(*,*) "call turbulent vertical diffusion ?"
253     calldifv=.true. ! default value
254     call getin_p("calldifv",calldifv)
255     write(*,*) " calldifv = ",calldifv
[135]256
[1524]257     write(*,*) "use turbdiff instead of vdifc ?"
258     UseTurbDiff=.true. ! default value
259     call getin_p("UseTurbDiff",UseTurbDiff)
260     write(*,*) " UseTurbDiff = ",UseTurbDiff
[596]261
[1524]262     write(*,*) "call convective adjustment ?"
263     calladj=.true. ! default value
264     call getin_p("calladj",calladj)
265     write(*,*) " calladj = ",calladj
[135]266
[2060]267     write(*,*) "call thermal plume model ?"
268     calltherm=.false. ! default value
269     call getin_p("calltherm",calltherm)
270     write(*,*) " calltherm = ",calltherm
271
[1524]272     write(*,*) "call CO2 condensation ?"
273     co2cond=.false. ! default value
274     call getin_p("co2cond",co2cond)
275     write(*,*) " co2cond = ",co2cond
[650]276! Test of incompatibility
[1524]277     if (co2cond.and.(.not.tracer)) then
278        print*,'We need a CO2 ice tracer to condense CO2'
279        call abort
280     endif
[716]281 
[1524]282     write(*,*) "CO2 supersaturation level ?"
283     co2supsat=1.0 ! default value
284     call getin_p("co2supsat",co2supsat)
285     write(*,*) " co2supsat = ",co2supsat
[253]286
[1524]287     write(*,*) "Radiative timescale for Newtonian cooling ?"
288     tau_relax=30. ! default value
289     call getin_p("tau_relax",tau_relax)
290     write(*,*) " tau_relax = ",tau_relax
291     tau_relax=tau_relax*24*3600 ! convert Earth days --> seconds
[253]292
[1524]293     write(*,*)"call thermal conduction in the soil ?"
294     callsoil=.true. ! default value
295     call getin_p("callsoil",callsoil)
296     write(*,*) " callsoil = ",callsoil
[135]297         
[1524]298     write(*,*)"Rad transfer is computed every iradia", &
299                   " physical timestep"
300     iradia=1 ! default value
301     call getin_p("iradia",iradia)
302     write(*,*)" iradia = ",iradia
[135]303       
[1524]304     write(*,*)"Rayleigh scattering ?"
305     rayleigh=.false.
306     call getin_p("rayleigh",rayleigh)
307     write(*,*)" rayleigh = ",rayleigh
[135]308
[1524]309     write(*,*) "Use blackbody for stellar spectrum ?"
310     stelbbody=.false. ! default value
311     call getin_p("stelbbody",stelbbody)
312     write(*,*) " stelbbody = ",stelbbody
[135]313
[1524]314     write(*,*) "Stellar blackbody temperature ?"
315     stelTbb=5800.0 ! default value
316     call getin_p("stelTbb",stelTbb)
317     write(*,*) " stelTbb = ",stelTbb
[253]318
[1524]319     write(*,*)"Output mean OLR in 1D?"
320     meanOLR=.false.
321     call getin_p("meanOLR",meanOLR)
322     write(*,*)" meanOLR = ",meanOLR
[135]323
[1524]324     write(*,*)"Output spectral OLR in 3D?"
325     specOLR=.false.
326     call getin_p("specOLR",specOLR)
327     write(*,*)" specOLR = ",specOLR
[135]328
[1524]329     write(*,*)"Operate in kastprof mode?"
330     kastprof=.false.
331     call getin_p("kastprof",kastprof)
332     write(*,*)" kastprof = ",kastprof
[253]333
[1524]334     write(*,*)"Uniform absorption in radiative transfer?"
335     graybody=.false.
336     call getin_p("graybody",graybody)
337     write(*,*)" graybody = ",graybody
[253]338
[1538]339! Soil model
340     write(*,*)"Number of sub-surface layers for soil scheme?"
341     ! default value of nsoilmx set in comsoil_h
342     call getin_p("nsoilmx",nsoilmx)
343     write(*,*)" nsoilmx=",nsoilmx
344     
345     write(*,*)"Thickness of topmost soil layer (m)?"
346     ! default value of lay1_soil set in comsoil_h
347     call getin_p("lay1_soil",lay1_soil)
348     write(*,*)" lay1_soil = ",lay1_soil
349     
350     write(*,*)"Coefficient for soil layer thickness distribution?"
351     ! default value of alpha_soil set in comsoil_h
352     call getin_p("alpha_soil",alpha_soil)
353     write(*,*)" alpha_soil = ",alpha_soil
354
[1297]355! Slab Ocean
[1524]356     write(*,*) "Use slab-ocean ?"
357     ok_slab_ocean=.false.         ! default value
358     call getin_p("ok_slab_ocean",ok_slab_ocean)
359     write(*,*) "ok_slab_ocean = ",ok_slab_ocean
[1529]360     ! Sanity check: for now slab oncean only works in serial mode
361     if (ok_slab_ocean.and.is_parallel) then
362       write(*,*) " Error: slab ocean should only be used in serial mode!"
363       call abort
364     endif
[1297]365
[1524]366     write(*,*) "Use slab-sea-ice ?"
367     ok_slab_sic=.true.         ! default value
368     call getin_p("ok_slab_sic",ok_slab_sic)
369     write(*,*) "ok_slab_sic = ",ok_slab_sic
[1297]370
[1524]371     write(*,*) "Use heat transport for the ocean ?"
372     ok_slab_heat_transp=.true.   ! default value
373     call getin_p("ok_slab_heat_transp",ok_slab_heat_transp)
374     write(*,*) "ok_slab_heat_transp = ",ok_slab_heat_transp
[1297]375
[1801]376! Photochemistry and chemistry in the thermosphere
[1297]377
[1801]378     write(*,*) "Use photochemistry ?"
379     photochem=.false.         ! default value
380     call getin_p("photochem",photochem)
381     write(*,*) "photochem = ",photochem
[1297]382
[1801]383     write(*,*)"Production of haze ?"
384     haze=.false. ! default value
385     call getin_p("haze",haze)
386     write(*,*)" haze = ",haze
387
388
[253]389! Test of incompatibility:
390! if kastprof used, we must be in 1D
[1524]391     if (kastprof.and.(ngrid.gt.1)) then
392       print*,'kastprof can only be used in 1D!'
393       call abort
394     endif
[253]395
[1524]396     write(*,*)"Stratospheric temperature for kastprof mode?"
397     Tstrat=167.0
398     call getin_p("Tstrat",Tstrat)
399     write(*,*)" Tstrat = ",Tstrat
[253]400
[1524]401     write(*,*)"Remove lower boundary?"
402     nosurf=.false.
403     call getin_p("nosurf",nosurf)
404     write(*,*)" nosurf = ",nosurf
[253]405
406! Tests of incompatibility:
[1524]407     if (nosurf.and.callsoil) then
408       print*,'nosurf not compatible with soil scheme!'
409       print*,'... got to make a choice!'
410       call abort
411     endif
[253]412
[1524]413     write(*,*)"Add an internal heat flux?", &
414                   "... matters only if callsoil=F"
415     intheat=0.
416     call getin_p("intheat",intheat)
417     write(*,*)" intheat = ",intheat
[952]418
[1524]419     write(*,*)"Use Newtonian cooling for radiative transfer?"
420     newtonian=.false.
421     call getin_p("newtonian",newtonian)
422     write(*,*)" newtonian = ",newtonian
[253]423
424! Tests of incompatibility:
[1524]425     if (newtonian.and.corrk) then
426       print*,'newtonian not compatible with correlated-k!'
427       call abort
428     endif
429     if (newtonian.and.calladj) then
430       print*,'newtonian not compatible with adjustment!'
431       call abort
432     endif
433     if (newtonian.and.calldifv) then
434       print*,'newtonian not compatible with a boundary layer!'
435       call abort
436     endif
[253]437
[1524]438     write(*,*)"Test physics timescale in 1D?"
439     testradtimes=.false.
440     call getin_p("testradtimes",testradtimes)
441     write(*,*)" testradtimes = ",testradtimes
[253]442
443! Test of incompatibility:
444! if testradtimes used, we must be in 1D
[1524]445     if (testradtimes.and.(ngrid.gt.1)) then
446       print*,'testradtimes can only be used in 1D!'
447       call abort
448     endif
[253]449
[1524]450     write(*,*)"Default planetary temperature?"
451     tplanet=215.0
452     call getin_p("tplanet",tplanet)
453     write(*,*)" tplanet = ",tplanet
[135]454
[1524]455     write(*,*)"Which star?"
456     startype=1 ! default value = Sol
457     call getin_p("startype",startype)
458     write(*,*)" startype = ",startype
[135]459
[1524]460     write(*,*)"Value of stellar flux at 1 AU?"
461     Fat1AU=1356.0 ! default value = Sol today
462     call getin_p("Fat1AU",Fat1AU)
463     write(*,*)" Fat1AU = ",Fat1AU
[135]464
465
466! TRACERS:
467
[1524]468     write(*,*)"Varying H2O cloud fraction?"
469     CLFvarying=.false.     ! default value
470     call getin_p("CLFvarying",CLFvarying)
471     write(*,*)" CLFvarying = ",CLFvarying
[253]472
[1524]473     write(*,*)"Value of fixed H2O cloud fraction?"
474     CLFfixval=1.0                ! default value
475     call getin_p("CLFfixval",CLFfixval)
476     write(*,*)" CLFfixval = ",CLFfixval
[253]477
[1524]478     write(*,*)"fixed radii for Cloud particles?"
479     radfixed=.false. ! default value
480     call getin_p("radfixed",radfixed)
481     write(*,*)" radfixed = ",radfixed
[728]482
[1524]483     if(kastprof)then
484        radfixed=.true.
485     endif 
[728]486
[1524]487     write(*,*)"Number mixing ratio of CO2 ice particles:"
488     Nmix_co2=1.e6 ! default value
489     call getin_p("Nmix_co2",Nmix_co2)
490     write(*,*)" Nmix_co2 = ",Nmix_co2
[135]491
[726]492!         write(*,*)"Number of radiatively active aerosols:"
493!         naerkind=0. ! default value
[1315]494!         call getin_p("naerkind",naerkind)
[726]495!         write(*,*)" naerkind = ",naerkind
496
[1524]497     write(*,*)"Opacity of dust (if used):"
498     dusttau=0. ! default value
499     call getin_p("dusttau",dusttau)
500     write(*,*)" dusttau = ",dusttau
[253]501
[1524]502     write(*,*)"Radiatively active CO2 aerosols?"
503     aeroco2=.false.     ! default value
504     call getin_p("aeroco2",aeroco2)
505     write(*,*)" aeroco2 = ",aeroco2
[253]506
[1524]507     write(*,*)"Fixed CO2 aerosol distribution?"
508     aerofixco2=.false.     ! default value
509     call getin_p("aerofixco2",aerofixco2)
510     write(*,*)" aerofixco2 = ",aerofixco2
[726]511
[1524]512     write(*,*)"Radiatively active water ice?"
513     aeroh2o=.false.     ! default value
514     call getin_p("aeroh2o",aeroh2o)
515     write(*,*)" aeroh2o = ",aeroh2o
[726]516
[1524]517     write(*,*)"Fixed H2O aerosol distribution?"
518     aerofixh2o=.false.     ! default value
519     call getin_p("aerofixh2o",aerofixh2o)
520     write(*,*)" aerofixh2o = ",aerofixh2o
[726]521
[1677]522     write(*,*)"Radiatively active sulfuric acid aerosols?"
[1524]523     aeroh2so4=.false.     ! default value
524     call getin_p("aeroh2so4",aeroh2so4)
525     write(*,*)" aeroh2so4 = ",aeroh2so4
[1026]526         
[1151]527!=================================
528
[1677]529     write(*,*)"Radiatively active two-layer aerosols?"
[1524]530     aeroback2lay=.false.     ! default value
531     call getin_p("aeroback2lay",aeroback2lay)
532     write(*,*)" aeroback2lay = ",aeroback2lay
[726]533
[1677]534     write(*,*)"Radiatively active ammonia cloud?"
535     aeronh3=.false.     ! default value
536     call getin_p("aeronh3",aeronh3)
537     write(*,*)" aeronh3 = ",aeronh3
538
539     write(*,*)"Radiatively active auroral aerosols?"
540     aeroaurora=.false.     ! default value
541     call getin_p("aeroaurora",aeroaurora)
542     write(*,*)" aeroaurora = ",aeroaurora
543
[1524]544     write(*,*)"TWOLAY AEROSOL: total optical depth ", &
545                    "in the tropospheric layer (visible)"
546     obs_tau_col_tropo=8.D0
547     call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
548     write(*,*)" obs_tau_col_tropo = ",obs_tau_col_tropo
[1026]549
[1524]550     write(*,*)"TWOLAY AEROSOL: total optical depth ", &
551                    "in the stratospheric layer (visible)"
552     obs_tau_col_strato=0.08D0
553     call getin_p("obs_tau_col_strato",obs_tau_col_strato)
554     write(*,*)" obs_tau_col_strato = ",obs_tau_col_strato
[1151]555
[2062]556     write(*,*)"TWOLAY AEROSOL: optprop_back2lay_vis?"
557     optprop_back2lay_vis = 'optprop_saturn_vis_n20.dat'
558     call getin_p("optprop_back2lay_vis",optprop_back2lay_vis)
559     write(*,*)" optprop_back2lay_vis = ",optprop_back2lay_vis
560
561     write(*,*)"TWOLAY AEROSOL: optprop_back2lay_ir?"
562     optprop_back2lay_ir = 'optprop_saturn_ir_n20.dat'
563     call getin_p("optprop_back2lay_ir",optprop_back2lay_ir)
564     write(*,*)" optprop_back2lay_ir = ",optprop_back2lay_ir
565     
[1524]566     write(*,*)"TWOLAY AEROSOL: pres_bottom_tropo? in pa"
567     pres_bottom_tropo=66000.0
568     call getin_p("pres_bottom_tropo",pres_bottom_tropo)
569     write(*,*)" pres_bottom_tropo = ",pres_bottom_tropo
[1151]570
[1524]571     write(*,*)"TWOLAY AEROSOL: pres_top_tropo? in pa"
572     pres_top_tropo=18000.0
573     call getin_p("pres_top_tropo",pres_top_tropo)
574     write(*,*)" pres_top_tropo = ",pres_top_tropo
[1151]575
[1524]576     write(*,*)"TWOLAY AEROSOL: pres_bottom_strato? in pa"
577     pres_bottom_strato=2000.0
578     call getin_p("pres_bottom_strato",pres_bottom_strato)
579     write(*,*)" pres_bottom_strato = ",pres_bottom_strato
[1151]580
[1524]581     write(*,*)"TWOLAY AEROSOL: pres_top_strato? in pa"
582     pres_top_strato=100.0
583     call getin_p("pres_top_strato",pres_top_strato)
584     write(*,*)" pres_top_strato = ",pres_top_strato
[1151]585
[1524]586     write(*,*)"TWOLAY AEROSOL: particle size in the ", &
587                    "tropospheric layer, in meters"
588     size_tropo=2.e-6
589     call getin_p("size_tropo",size_tropo)
590     write(*,*)" size_tropo = ",size_tropo
[1151]591
[1524]592     write(*,*)"TWOLAY AEROSOL: particle size in the ", &
593                    "stratospheric layer, in meters"
594     size_strato=1.e-7
595     call getin_p("size_strato",size_strato)
596     write(*,*)" size_strato = ",size_strato
[1151]597
[1677]598     write(*,*)"NH3 (thin) cloud: total optical depth"
599     tau_nh3_cloud=7.D0
600     call getin_p("tau_nh3_cloud",tau_nh3_cloud)
601     write(*,*)" tau_nh3_cloud = ",tau_nh3_cloud
602
603     write(*,*)"NH3 (thin) cloud pressure level"
604     pres_nh3_cloud=7.D0
605     call getin_p("pres_nh3_cloud",pres_nh3_cloud)
606     write(*,*)" pres_nh3_cloud = ",pres_nh3_cloud
607
608     write(*,*)"NH3 (thin) cloud: particle sizes"
609     size_nh3_cloud=3.e-6
610     call getin_p("size_nh3_cloud",size_nh3_cloud)
611     write(*,*)" size_nh3_cloud = ",size_nh3_cloud
612
[1151]613!=================================
614
[1524]615     write(*,*)"Cloud pressure level (with kastprof only):"
616     cloudlvl=0. ! default value
617     call getin_p("cloudlvl",cloudlvl)
618     write(*,*)" cloudlvl = ",cloudlvl
[305]619
[1524]620     write(*,*)"Is the variable gas species radiatively active?"
621     Tstrat=167.0
622     varactive=.false.
623     call getin_p("varactive",varactive)
624     write(*,*)" varactive = ",varactive
[135]625
[1524]626     write(*,*)"Is the variable gas species distribution set?"
627     varfixed=.false.
628     call getin_p("varfixed",varfixed)
629     write(*,*)" varfixed = ",varfixed
[135]630
[1524]631     write(*,*)"What is the saturation % of the variable species?"
632     satval=0.8
633     call getin_p("satval",satval)
634     write(*,*)" satval = ",satval
[135]635
636
637! Test of incompatibility:
638! if varactive, then varfixed should be false
[1524]639     if (varactive.and.varfixed) then
640       print*,'if varactive, varfixed must be OFF!'
641       stop
642     endif
[135]643
[1524]644     write(*,*) "Gravitationnal sedimentation ?"
645     sedimentation=.false. ! default value
646     call getin_p("sedimentation",sedimentation)
647     write(*,*) " sedimentation = ",sedimentation
[135]648
[1524]649     write(*,*) "Compute water cycle ?"
650     water=.false. ! default value
651     call getin_p("water",water)
652     write(*,*) " water = ",water
[135]653         
[622]654! Test of incompatibility:
655! if water is true, there should be at least a tracer
[1524]656     if (water.and.(.not.tracer)) then
657        print*,'if water is ON, tracer must be ON too!'
658        stop
659     endif
[622]660
[1524]661     write(*,*) "Include water condensation ?"
662     watercond=.false. ! default value
663     call getin_p("watercond",watercond)
664     write(*,*) " watercond = ",watercond
[135]665
[622]666! Test of incompatibility:
667! if watercond is used, then water should be used too
[1524]668     if (watercond.and.(.not.water)) then
669        print*,'if watercond is used, water should be used too'
670        stop
671     endif
[622]672
[1524]673     write(*,*) "Include water precipitation ?"
674     waterrain=.false. ! default value
675     call getin_p("waterrain",waterrain)
676     write(*,*) " waterrain = ",waterrain
[135]677
[1524]678     write(*,*) "Include surface hydrology ?"
679     hydrology=.false. ! default value
680     call getin_p("hydrology",hydrology)
681     write(*,*) " hydrology = ",hydrology
[135]682
[1524]683     write(*,*) "Evolve surface water sources ?"
684     sourceevol=.false. ! default value
685     call getin_p("sourceevol",sourceevol)
686     write(*,*) " sourceevol = ",sourceevol
[135]687
[1524]688     write(*,*) "Ice evolution timestep ?"
689     icetstep=100.0 ! default value
690     call getin_p("icetstep",icetstep)
691     write(*,*) " icetstep = ",icetstep
[1498]692         
[1524]693     write(*,*) "Spectral Dependant albedo ?"
694     albedo_spectral_mode=.false. ! default value
695     call getin_p("albedo_spectral_mode",albedo_spectral_mode)
696     write(*,*) " albedo_spectral_mode = ",albedo_spectral_mode
[486]697
[1524]698     write(*,*) "Snow albedo ?"
699     write(*,*) "If albedo_spectral_mode=.true., then this "
700     write(*,*) "corresponds to the 0.5 microns snow albedo."
701     albedosnow=0.5         ! default value
702     call getin_p("albedosnow",albedosnow)
703     write(*,*) " albedosnow = ",albedosnow
[1482]704         
[1524]705     write(*,*) "CO2 ice albedo ?"
706     albedoco2ice=0.5       ! default value
707     call getin_p("albedoco2ice",albedoco2ice)
708     write(*,*) " albedoco2ice = ",albedoco2ice
[135]709
[1524]710     write(*,*) "Maximum ice thickness ?"
711     maxicethick=2.0         ! default value
712     call getin_p("maxicethick",maxicethick)
713     write(*,*) " maxicethick = ",maxicethick
[135]714
[1524]715     write(*,*) "Freezing point of seawater ?"
716     Tsaldiff=-1.8          ! default value
717     call getin_p("Tsaldiff",Tsaldiff)
718     write(*,*) " Tsaldiff = ",Tsaldiff
[135]719
[1524]720     write(*,*) "Does user want to force cpp and mugaz?"
721     force_cpp=.false. ! default value
722     call getin_p("force_cpp",force_cpp)
723     write(*,*) " force_cpp = ",force_cpp
[589]724
[1524]725     if (force_cpp) then
726       mugaz = -99999.
727       PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
728       call getin_p("mugaz",mugaz)
729       IF (mugaz.eq.-99999.) THEN
730           PRINT *, "mugaz must be set if force_cpp = T"
731           STOP
732       ELSE
733           write(*,*) "mugaz=",mugaz
734       ENDIF
735       cpp = -99999.
736       PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
737       call getin_p("cpp",cpp)
738       IF (cpp.eq.-99999.) THEN
739           PRINT *, "cpp must be set if force_cpp = T"
740           STOP
741       ELSE
742           write(*,*) "cpp=",cpp
743       ENDIF
744     endif ! of if (force_cpp)
745     call su_gases
746     call calc_cpp_mugaz
[135]747
[1524]748     PRINT*,'--------------------------------------------'
749     PRINT*
750     PRINT*
751  ELSE
752     write(*,*)
753     write(*,*) 'Cannot read file callphys.def. Is it here ?'
754     stop
755  ENDIF ! of IF(iscallphys)
[135]756
[1524]757  PRINT*
758  PRINT*,'inifis: daysec',daysec
759  PRINT*
760  PRINT*,'inifis: The radiative transfer is computed:'
761  PRINT*,'           each ',iradia,' physical time-step'
762  PRINT*,'        or each ',iradia*dtphys,' seconds'
763  PRINT*
[135]764
765
766!-----------------------------------------------------------------------
767!     Some more initialization:
768!     ------------------------
769
[1542]770  ! Initializations for comgeomfi_h
[1832]771#ifndef MESOSCALE
[1542]772  totarea=SSUM(ngrid,parea,1)
773  call planetwide_sumval(parea,totarea_planet)
[787]774
[1524]775  !! those are defined in comdiurn_h.F90
776  IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
777  IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
778  IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
779  IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
[787]780
[1524]781  DO ig=1,ngrid
782     sinlat(ig)=sin(plat(ig))
783     coslat(ig)=cos(plat(ig))
784     sinlon(ig)=sin(plon(ig))
785     coslon(ig)=cos(plon(ig))
786  ENDDO
[1832]787#endif
[1529]788  ! initialize variables in radinc_h
789  call ini_radinc_h(nlayer)
790 
[1524]791  ! allocate "comsoil_h" arrays
792  call ini_comsoil_h(ngrid)
[1801]793   
[1524]794  END SUBROUTINE inifis
[135]795
[1524]796END MODULE inifis_mod
Note: See TracBrowser for help on using the repository browser.