source: trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90 @ 3949

Last change on this file since 3949 was 3949, checked in by debatzbr, 7 weeks ago

Pluto PCM: Add condensable gas tracers through muphi
BBT

File size: 57.6 KB
Line 
1MODULE inifis_mod
2IMPLICIT NONE
3
4CONTAINS
5
6  SUBROUTINE inifis(ngrid,nlayer,nq, &
7             day_ini,pdaysec,nday,ptimestep, &
8             plat,plon,parea, &
9             prad,pg,pr,pcpp)
10
11  use init_print_control_mod, only: init_print_control
12  use radinc_h, only: ini_radinc_h, naerkind
13  use radcommon_h, only: ini_radcommon_h
14  use radii_mod, only: radfixed
15  use datafile_mod, only: datadir,hazeprop_file,hazerad_file,hazemmr_file,hazedens_file, &
16                          config_mufi, mugasflux_file, aersprop_file, aerfprop_file
17  use comdiurn_h, only: sinlat, coslat, sinlon, coslon
18  use comgeomfi_h, only: totarea, totarea_planet
19  use comsoil_h, only: ini_comsoil_h, nsoilmx, lay1_soil, alpha_soil
20  use time_phylmdz_mod, only: diagfi_output_rate,startfi_output_rate, &
21                              init_time, daysec, dtphys
22  use comcstfi_mod, only: rad, cpp, g, r, rcp, &
23                          mugaz, pi, avocado
24  use planetwide_mod, only: planetwide_sumval
25  use nonoro_gwd_mix_mod, only: calljliu_gwimix
26  use callkeys_mod
27  use nonoro_gwd_ran_mod, only: ini_nonoro_gwd_ran, &
28                                end_nonoro_gwd_ran
29  use nonoro_gwd_mix_mod, only: ini_nonoro_gwd_mix, &
30                                end_nonoro_gwd_mix
31                               
32  use surfdat_h
33  use wstats_mod, only: callstats
34  use writediagsoil_mod, only: diagsoil
35  use ioipsl_getin_p_mod, only : getin_p
36  use mod_phys_lmdz_para, only : is_parallel, is_master, bcast
37
38!=======================================================================
39!
40!   purpose:
41!   -------
42!
43!   Initialisation for the physical parametrisations of the LMD
44!   Generic Model.
45!
46!   author: Frederic Hourdin 15 / 10 /93
47!   -------
48!   modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
49!             Ehouarn Millour (oct. 2008) tracers are now identified
50!              by their names and may not be contiguously
51!              stored in the q(:,:,:,:) array
52!             E.M. (june 2009) use getin routine to load parameters
53!
54!
55!   arguments:
56!   ----------
57!
58!   input:
59!   ------
60!
61!    ngrid                 Size of the horizontal grid.
62!                          All internal loops are performed on that grid.
63!    nlayer                Number of vertical layers.
64!    pdayref               Day of reference for the simulation
65!    pday                  Number of days counted from the North. Spring
66!                          equinoxe.
67!
68!=======================================================================
69!
70!-----------------------------------------------------------------------
71!   declarations:
72!   -------------
73  use datafile_mod, only: datadir
74  use ioipsl_getin_p_mod, only: getin_p
75  IMPLICIT NONE
76
77
78
79  REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep
80  INTEGER,INTENT(IN) :: nday
81  INTEGER,INTENT(IN) :: ngrid,nlayer,nq
82  REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid)
83  integer,intent(in) :: day_ini
84
85  INTEGER ig
86  CHARACTER(len=20) :: rname="inifis" ! routine name, for messages
87
88  EXTERNAL iniorbit,orbite
89  EXTERNAL SSUM
90  REAL SSUM
91
92  ! deprecated parameters
93  REAL :: ecritphy ! to check that this obsolete flag is no longer used...
94  logical aerohaze
95
96  ! Initialize flags lunout, prt_level, debug (in print_control_mod)
97  CALL init_print_control
98
99  ! initialize constants in comcstfi_mod
100  rad=prad
101  cpp=pcpp
102  g=pg
103  r=pr
104  rcp=r/cpp
105  mugaz=8.314*1000./pr ! dummy init
106  pi=2.*asin(1.)
107  avocado = 6.02214179e23   ! added by RW
108
109  ! Initialize some "temporal and calendar" related variables
110#ifndef MESOSCALE
111  CALL init_time(day_ini,pdaysec,nday,ptimestep)
112#endif
113
114
115  ! read in some parameters from "run.def" for physics,
116  ! or shared between dynamics and physics.
117  ecritphy=-666 ! dummy default value
118  call getin_p("ecritphy",ecritphy) ! frequency of outputs in physics,
119                                    ! in dynamical steps
120  if (ecritphy/=-666) then
121    call abort_physic(rname, &
122         "Error: parameter ecritphy is obsolete! Remove it! "//&
123         "And use diagfi_output_rate instead",1)
124  endif
125
126  ! do we read a startphy.nc file? (default: .true.)
127  call getin_p("startphy_file",startphy_file)
128
129! --------------------------------------------------------------
130!  Reading the "callphys.def" file controlling some key options
131! --------------------------------------------------------------
132
133  IF (is_master) THEN
134    ! check if 'callphys.def' file is around
135    inquire(file="callphys.def",exist=iscallphys)
136    write(*,*) trim(rname)//": iscallphys=",iscallphys
137  ENDIF
138  call bcast(iscallphys)
139
140  IF(iscallphys) THEN
141
142     if (is_master) then
143       write(*,*)
144       write(*,*)'--------------------------------------------'
145       write(*,*)trim(rname)//': Parameters for the physics (callphys.def)'
146       write(*,*)'--------------------------------------------'
147     endif
148
149     if (is_master) write(*,*) trim(rname)//&
150       ": Directory where external input files are:"
151     ! default 'datadir' is set in "datadir_mod"
152     call getin_p("datadir",datadir) ! default path
153     if (is_master) write(*,*) trim(rname)//": datadir = ",trim(datadir)
154
155     if (is_master) write(*,*) "Haze optical properties datafile"
156     hazeprop_file="optprop_tholins_fractal050nm.dat"  ! default file
157     call getin_p("hazeprop_file",hazeprop_file)
158     if (is_master) write(*,*) trim(rname)//" hazeprop_file = ",trim(hazeprop_file)
159
160     if (is_master) write(*,*) "Haze radii datafile"
161     hazerad_file="hazerad.txt"  ! default file
162     call getin_p("hazerad_file",hazerad_file)
163     if (is_master) write(*,*) trim(rname)//" hazerad_file = ",trim(hazerad_file)
164
165     if (is_master) write(*,*) "Haze mmr datafile"
166     hazemmr_file="None"  ! default file
167     call getin_p("hazemmr_file",hazemmr_file)
168     if (is_master) write(*,*) trim(rname)//" hazemmr_file = ",trim(hazemmr_file)
169     if (is_master) write(*,*) "Haze density datafile"
170     hazedens_file="None"  ! default file
171     call getin_p("hazedens_file",hazedens_file)
172     if (is_master) write(*,*) trim(rname)//" hazedens_file = ",trim(hazedens_file)
173
174     if (is_master) write(*,*) trim(rname)//&
175       ": Run with or without tracer transport ?"
176     tracer=.false. ! default value
177     call getin_p("tracer",tracer)
178     if (is_master) write(*,*) trim(rname)//": tracer = ",tracer
179
180     if (is_master) write(*,*) trim(rname)//&
181       ": Output rate for diagfi.nc file (in physics steps) ?"
182     diagfi_output_rate=24 !default value
183     call getin_p("diagfi_output_rate",diagfi_output_rate)
184     if (is_master) write(*,*) trim(rname)//": diagfi_output_rate = ",&
185                               diagfi_output_rate
186
187     if (is_master) write(*,*) trim(rname)//&
188       ": Output rate for start/startfi.nc files (in physics steps) ?"
189     startfi_output_rate=0 !default value
190     call getin_p("startfi_output_rate",startfi_output_rate)
191     if (is_master) write(*,*) trim(rname)//": startfi_output_rate = ",&
192                               startfi_output_rate
193
194     if (is_master) write(*,*) trim(rname)//&
195       ": Run with or without atm mass update "//&
196       " due to tracer evaporation/condensation?"
197     mass_redistrib=.false. ! default value
198     call getin_p("mass_redistrib",mass_redistrib)
199     if (is_master) write(*,*) trim(rname)//&
200       ": mass_redistrib = ",mass_redistrib
201
202     if (is_master) then
203       write(*,*) trim(rname)//": Diurnal cycle ?"
204       write(*,*) trim(rname)//&
205       ": (if diurnal=false, diurnal averaged solar heating)"
206     endif
207     diurnal=.true. ! default value
208     call getin_p("diurnal",diurnal)
209     if (is_master) write(*,*) trim(rname)//": diurnal = ",diurnal
210
211     if (is_master) then
212       write(*,*) trim(rname)//": Seasonal cycle ?"
213       write(*,*) trim(rname)//&
214       ": (if season=false, Ls stays constant, to value ", &
215       "set in 'start'"
216     endif
217     season=.true. ! default value
218     call getin_p("season",season)
219     if (is_master) write(*,*) trim(rname)//": season = ",season
220
221     evol1d=.false. ! default value
222     call getin_p("evol1d",evol1d)
223     if (is_master) write(*,*) trim(rname)//": evol1d = ",evol1d
224
225     if (is_master) then
226       write(*,*) trim(rname)//": Fast mode (nogcm) ?"
227     endif
228     fast=.false. ! default value
229     call getin_p("fast",fast)
230     if (is_master) write(*,*) trim(rname)//": fast = ",fast
231
232     if (is_master) write(*,*) trim(rname)//"Run with Triton orbit ?"
233     triton=.false. ! default value
234     call getin_p("triton",triton)
235     if (is_master) write(*,*) trim(rname)//" triton = ",triton
236     convergeps=.false. ! default value
237     call getin_p("convergeps",convergeps)
238     if (is_master) write(*,*) trim(rname)//" convergeps = ",convergeps
239
240     conservn2=.false. ! default value
241     call getin_p("conservn2",conservn2)
242     if (is_master) write(*,*) trim(rname)//" conservn2 = ",conservn2
243     conservch4=.false. ! default value
244     call getin_p("conservch4",conservch4)
245     if (is_master) write(*,*) trim(rname)//" conservch4 = ",conservch4
246
247     if (is_master) write(*,*) trim(rname)//"KBO runs (eris, makemake) ?"
248     kbo=.false. ! default value
249     call getin_p("kbo",kbo)
250     if (is_master) write(*,*) trim(rname)//" kbo = ",kbo
251
252     if (is_master) write(*,*) trim(rname)//"Specific paleo run ?"
253     paleo=.false. ! default value
254     call getin_p("paleo",paleo)
255     if (is_master) write(*,*) trim(rname)//" paleo = ",paleo
256
257     if (is_master) write(*,*) trim(rname)//"paleoclimate step (Earth years) "
258     paleoyears=10000. ! default value
259     call getin_p("paleoyears",paleoyears)
260     if (is_master) write(*,*) trim(rname)//" paleoyears = ",paleoyears
261
262     if (is_master) write(*,*) trim(rname)//&
263       ": No seasonal cycle: initial day to lock the run during restart"
264     noseason_day=0.0 ! default value
265     call getin_p("noseason_day",noseason_day)
266     if (is_master) write(*,*) trim(rname)//": noseason_day=", noseason_day
267
268     if (is_master) write(*,*) trim(rname)//&
269       ": Write some extra output to the screen ?"
270     lwrite=.false. ! default value
271     call getin_p("lwrite",lwrite)
272     if (is_master) write(*,*) trim(rname)//": lwrite = ",lwrite
273
274     if (is_master) write(*,*) trim(rname)//&
275       ": Save statistics in file stats.nc ?"
276     callstats=.true. ! default value
277     call getin_p("callstats",callstats)
278     if (is_master) write(*,*) trim(rname)//": callstats = ",callstats
279
280     if (is_master) write(*,*) trim(rname)//&
281       ": Write sub-surface fields in file diagsoil.nc ?"
282     diagsoil=.false. ! default value
283     call getin_p("diagsoil",diagsoil)
284     if (is_master) write(*,*) trim(rname)//" diagsoil = ",diagsoil
285
286     if (is_master) write(*,*) trim(rname)//&
287       ": Test energy conservation of model physics ?"
288     enertest=.false. ! default value
289     call getin_p("enertest",enertest)
290     if (is_master) write(*,*) trim(rname)//": enertest = ",enertest
291
292     if (is_master) write(*,*) trim(rname)//": call radiative transfer ?"
293     callrad=.true. ! default value
294     call getin_p("callrad",callrad)
295     if (is_master) write(*,*) trim(rname)//": callrad = ",callrad
296
297     if (is_master) write(*,*) trim(rname)//&
298       ": call correlated-k radiative transfer ?"
299     corrk=.true. ! default value
300     call getin_p("corrk",corrk)
301     if (is_master) write(*,*) trim(rname)//": corrk = ",corrk
302
303     if (is_master) write(*,*) trim(rname)//&
304       ": prohibit calculations outside corrk T grid?"
305     strictboundcorrk=.true. ! default value
306     call getin_p("strictboundcorrk",strictboundcorrk)
307     if (is_master) write(*,*) trim(rname)//&
308       ": strictboundcorrk = ",strictboundcorrk
309
310     if (is_master) write(*,*) trim(rname)//&
311       ": prohibit calculations outside CIA T grid?"
312     strictboundcia=.true. ! default value
313     call getin_p("strictboundcia",strictboundcia)
314     if (is_master) write(*,*) trim(rname)//&
315       ": strictboundcia = ",strictboundcia
316
317     if (is_master) write(*,*) trim(rname)//&
318       ": Minimum atmospheric temperature for Planck function integration ?"
319     tplanckmin=30.0 ! default value
320     call getin_p("tplanckmin",tplanckmin)
321     if (is_master) write(*,*) trim(rname)//": tplanckmin = ",tplanckmin
322
323     if (is_master) write(*,*) trim(rname)//&
324       ": Maximum atmospheric temperature for Planck function integration ?"
325     tplanckmax=1500.0 ! default value
326     call getin_p("tplanckmax",tplanckmax)
327     if (is_master) write(*,*) trim(rname)//": tplanckmax = ",tplanckmax
328
329     if (is_master) write(*,*) trim(rname)//&
330       ": Temperature step for Planck function integration ?"
331     dtplanck=0.1 ! default value
332     call getin_p("dtplanck",dtplanck)
333     if (is_master) write(*,*) trim(rname)//": dtplanck = ",dtplanck
334
335     if (is_master) write(*,*) trim(rname)//&
336       ": call gaseous absorption in the visible bands?"//&
337       " (matters only if callrad=T)"
338     callgasvis=.false. ! default value
339     call getin_p("callgasvis",callgasvis)
340     if (is_master) write(*,*) trim(rname)//": callgasvis = ",callgasvis
341
342     if (is_master) write(*,*) trim(rname)//&
343       ": call continuum opacities in radiative transfer ?"//&
344       " (matters only if callrad=T)"
345     continuum=.true. ! default value
346     call getin_p("continuum",continuum)
347     if (is_master) write(*,*) trim(rname)//": continuum = ",continuum
348
349     if (is_master) write(*,*) trim(rname)//&
350       ": call turbulent vertical diffusion ?"
351     calldifv=.false. ! default value
352     call getin_p("calldifv",calldifv)
353     if (is_master) write(*,*) trim(rname)//": calldifv = ",calldifv
354     vertdiff=.true. ! default value
355     call getin_p("vertdiff",vertdiff)
356     if (is_master) write(*,*) trim(rname)//": vertdiff = ",vertdiff
357
358     if (is_master) write(*,*) trim(rname)//": use turbdiff instead of vdifc ?"
359     UseTurbDiff=.false. ! default value
360     call getin_p("UseTurbDiff",UseTurbDiff)
361     if (is_master) write(*,*) trim(rname)//": UseTurbDiff = ",UseTurbDiff
362
363     if (is_master) write(*,*) trim(rname)//": call convective adjustment ?"
364     calladj=.false. ! default value
365     call getin_p("calladj",calladj)
366     if (is_master) write(*,*) trim(rname)//": calladj = ",calladj
367
368     if (is_master) write(*,*) trim(rname)//&
369     ": Radiative timescale for Newtonian cooling (in Earth days)?"
370     tau_relax=30. ! default value
371     call getin_p("tau_relax",tau_relax)
372     if (is_master) write(*,*) trim(rname)//": tau_relax = ",tau_relax
373     tau_relax=tau_relax*24*3600 ! convert Earth days --> seconds
374
375     if (is_master) write(*,*)trim(rname)//&
376       ": call thermal conduction in the soil ?"
377     callsoil=.true. ! default value
378     call getin_p("callsoil",callsoil)
379     if (is_master) write(*,*) trim(rname)//": callsoil = ",callsoil
380
381     if (is_master) write(*,*)trim(rname)//&
382       ": call orographic gravity waves ?"
383     calllott=.false. ! default value
384     call getin_p("calllott",calllott)
385     if (is_master) write(*,*)trim(rname)//": calllott = ",calllott
386
387     if (is_master) write(*,*)trim(rname)//&
388       ": call  non-orographic gravity waves ?"
389     calllott_nonoro=.false. ! default value
390     call getin_p("calllott_nonoro",calllott_nonoro)
391     if (is_master) write(*,*)trim(rname)//&
392       ": calllott_nonoro = ",calllott_nonoro
393
394     if (is_master) write(*,*)trim(rname)//&       
395       ": call jliu's non-orogrphic GW-induced turbulence"
396     calljliu_gwimix=.false. ! default value
397     call getin_p("calljliu_gwimix",calljliu_gwimix)
398     if (is_master) write(*,*)trim(rname)//&   
399       ": calljliu_gwimix = ",calljliu_gwimix
400
401     if (is_master) write(*,*)trim(rname)//&
402       ": Rad transfer is computed every iradia", &
403       " physical timestep"
404     iradia=1 ! default value
405     call getin_p("iradia",iradia)
406     if (is_master) write(*,*)trim(rname)//": iradia = ",iradia
407
408     if (is_master) write(*,*)trim(rname)//": Rayleigh scattering ?"
409     rayleigh=.false.
410     call getin_p("rayleigh",rayleigh)
411     if (is_master) write(*,*)trim(rname)//": rayleigh = ",rayleigh
412
413     if (is_master) write(*,*) trim(rname)//&
414       ": Use blackbody for stellar spectrum ?"
415     stelbbody=.false. ! default value
416     call getin_p("stelbbody",stelbbody)
417     if (is_master) write(*,*) trim(rname)//": stelbbody = ",stelbbody
418
419     if (is_master) write(*,*) trim(rname)//": Stellar blackbody temperature ?"
420     stelTbb=5800.0 ! default value
421     call getin_p("stelTbb",stelTbb)
422     if (is_master) write(*,*) trim(rname)//": stelTbb = ",stelTbb
423
424     if (is_master) write(*,*) trim(rname)//": Output mean OLR in 1D?"
425     meanOLR=.false.
426     call getin_p("meanOLR",meanOLR)
427     if (is_master) write(*,*) trim(rname)//": meanOLR = ",meanOLR
428
429     if (is_master) write(*,*)trim(rname)//": Output spectral OLR in 3D?"
430     specOLR=.false.
431     call getin_p("specOLR",specOLR)
432     if (is_master) write(*,*)trim(rname)//": specOLR = ",specOLR
433
434     if (is_master) write(*,*)trim(rname)//": Operate in kastprof mode?"
435     kastprof=.false.
436     call getin_p("kastprof",kastprof)
437     if (is_master) write(*,*)trim(rname)//": kastprof = ",kastprof
438
439     if (is_master) write(*,*)trim(rname)//&
440       ": Uniform absorption in radiative transfer?"
441     graybody=.false.
442     call getin_p("graybody",graybody)
443     if (is_master) write(*,*)trim(rname)//": graybody = ",graybody
444
445! Soil model
446     if (is_master) write(*,*)trim(rname)//&
447       ": Number of sub-surface layers for soil scheme?"
448     ! default value of nsoilmx set in comsoil_h
449     call getin_p("nsoilmx",nsoilmx)
450     if (is_master) write(*,*)trim(rname)//": nsoilmx=",nsoilmx
451
452     if (is_master) write(*,*)trim(rname)//&
453       ": Thickness of topmost soil layer (m)?"
454     ! default value of lay1_soil set in comsoil_h
455     call getin_p("lay1_soil",lay1_soil)
456     if (is_master) write(*,*)trim(rname)//": lay1_soil = ",lay1_soil
457
458     if (is_master) write(*,*)trim(rname)//&
459       ": Coefficient for soil layer thickness distribution?"
460     ! default value of alpha_soil set in comsoil_h
461     call getin_p("alpha_soil",alpha_soil)
462     if (is_master) write(*,*)trim(rname)//": alpha_soil = ",alpha_soil
463
464     if (is_master) write(*,*)trim(rname)//&
465       "Geothermal flux (W) at the bottom layer"
466     fluxgeo=0. ! default value
467     call getin_p("fluxgeo",fluxgeo)
468     if (is_master) write(*,*)trim(rname)//" fluxgeo = ",fluxgeo
469
470     if (is_master) write(*,*)trim(rname)//&
471       "Assymetry flux (W) at the bottom layer"
472     assymflux=.false. ! default value
473     call getin_p("assymflux",assymflux)
474     if (is_master) write(*,*)trim(rname)//" assymflux = ",assymflux
475
476     if (is_master) write(*,*)trim(rname)//&
477       "Geothermal flux (W) for assymetry"
478     fluxgeo2=fluxgeo ! default value
479     call getin_p("fluxgeo2",fluxgeo2)
480     if (is_master) write(*,*)trim(rname)//" fluxgeo2 = ",fluxgeo2
481
482     if (is_master) write(*,*)trim(rname)//&
483       "Warm patch of flux"
484     patchflux=0 ! default value
485     call getin_p("patchflux",patchflux)
486     if (is_master) write(*,*)trim(rname)//" patchflux = ",patchflux
487
488! Chemistry in the thermosphere
489     if (is_master) write(*,*) trim(rname)//": Use deposition ?"
490     depos=.false.         ! default value
491     call getin_p("depos",depos)
492     if (is_master) write(*,*) trim(rname)//": depos = ",depos
493
494     if (is_master) write(*,*)trim(rname)//": Production of haze ?"
495     haze=.false. ! default value
496     call getin_p("haze",haze)
497     if (is_master) write(*,*)trim(rname)//": haze = ",haze
498
499
500
501      if (is_master) write(*,*)trim(rname)// "call thermal conduction ?"
502      callconduct=.false. ! default value
503      call getin_p("callconduct",callconduct)
504      if (is_master) write(*,*)trim(rname)// " callconduct = ",callconduct
505
506      if (is_master) write(*,*)trim(rname)// "call phitop ?"
507      phitop=0. ! default value
508      call getin_p("phitop",phitop)
509      if (is_master) write(*,*)trim(rname)// " phitop = ",phitop
510
511      if (is_master) write(*,*)trim(rname)// "call molecular viscosity ?"
512      callmolvis=.false. ! default value
513      call getin_p("callmolvis",callmolvis)
514      if (is_master) write(*,*)trim(rname)// " callmolvis = ",callmolvis
515
516      if (is_master) write(*,*)trim(rname)// "call molecular diffusion ?"
517      callmoldiff=.false. ! default value
518      call getin_p("callmoldiff",callmoldiff)
519      if (is_master) write(*,*)trim(rname)// " callmoldiff = ",callmoldiff
520
521! Global1D mean and solar zenith angle
522
523     if (ngrid.eq.1) then
524      PRINT*, 'Simulate global averaged conditions ?'
525      global1d = .false. ! default value
526      call getin_p("global1d",global1d)
527      write(*,*) "global1d = ",global1d
528
529      ! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
530      if (global1d.and.diurnal) then
531         call abort_physic(rname,'if global1d is true, diurnal must be set to false',1)
532      endif
533
534      if (global1d) then
535         PRINT *,'Solar Zenith angle (deg.) ?'
536         PRINT *,'(assumed for averaged solar flux S/4)'
537         szangle=60.0  ! default value
538         call getin_p("szangle",szangle)
539         write(*,*) "szangle = ",szangle
540      endif
541   endif ! of if (ngrid.eq.1)
542
543! Test of incompatibility:
544! if kastprof used, we must be in 1D
545     if (kastprof.and.(ngrid.gt.1)) then
546       call abort_physic(rname,'kastprof can only be used in 1D!',1)
547     endif
548
549     if (is_master) write(*,*)trim(rname)//&
550       ": Stratospheric temperature for kastprof mode?"
551     Tstrat=167.0
552     call getin_p("Tstrat",Tstrat)
553     if (is_master) write(*,*)trim(rname)//": Tstrat = ",Tstrat
554
555     if (is_master) write(*,*)trim(rname)//": Remove lower boundary?"
556     nosurf=.false.
557     call getin_p("nosurf",nosurf)
558     if (is_master) write(*,*)trim(rname)//": nosurf = ",nosurf
559
560! Tests of incompatibility:
561     if (nosurf.and.callsoil) then
562       if (is_master) then
563         write(*,*)trim(rname)//'nosurf not compatible with soil scheme!'
564         write(*,*)trim(rname)//'... got to make a choice!'
565       endif
566       call abort_physic(rname,"nosurf not compatible with soil scheme!",1)
567     endif
568
569     if (is_master) write(*,*)trim(rname)//": Add an internal heat flux?", &
570                   "... matters only if callsoil=F"
571     intheat=0.
572     call getin_p("intheat",intheat)
573     if (is_master) write(*,*)trim(rname)//": intheat = ",intheat
574
575
576     if (is_master) write(*,*)trim(rname)//": Test physics timescale in 1D?"
577     testradtimes=.false.
578     call getin_p("testradtimes",testradtimes)
579     if (is_master) write(*,*)trim(rname)//": testradtimes = ",testradtimes
580
581! Test of incompatibility:
582! if testradtimes used, we must be in 1D
583     if (testradtimes.and.(ngrid.gt.1)) then
584       call abort_physic(rname,'testradtimes can only be used in 1D!',1)
585     endif
586
587     if (is_master) write(*,*)trim(rname)//": Default planetary temperature?"
588     tplanet=215.0
589     call getin_p("tplanet",tplanet)
590     if (is_master) write(*,*)trim(rname)//": tplanet = ",tplanet
591
592     if (is_master) write(*,*)trim(rname)//": Which star?"
593     startype=1 ! default value = Sol
594     call getin_p("startype",startype)
595     if (is_master) write(*,*)trim(rname)//": startype = ",startype
596
597     if (is_master) write(*,*)trim(rname)//": Value of stellar flux at 1 AU?"
598     Fat1AU=1356.0 ! default value = Sol today
599     call getin_p("Fat1AU",Fat1AU)
600     if (is_master) write(*,*)trim(rname)//": Fat1AU = ",Fat1AU
601
602
603! TRACERS:
604
605     if (is_master) write(*,*)trim(rname)//": fixed radii for Cloud particles?"
606     radfixed=.false. ! default value
607     call getin_p("radfixed",radfixed)
608     if (is_master) write(*,*)trim(rname)//": radfixed = ",radfixed
609
610     if(kastprof)then
611        radfixed=.true.
612     endif
613
614      if (is_master) write(*,*)trim(rname)//&
615       "Number of radiatively active aerosols:"
616     naerkind=0 ! default value
617     call getin_p("naerkind",naerkind)
618     if (is_master) write(*,*)trim(rname)//": naerkind = ",naerkind
619
620
621!***************************************************************
622         !! TRACERS options
623
624     if (is_master)write(*,*)trim(rname)//&
625      "call N2 condensation ?"
626     n2cond=.true. ! default value
627     call getin_p("n2cond",n2cond)
628     if (is_master)write(*,*)trim(rname)//&
629      " n2cond = ",n2cond
630
631     if (is_master)write(*,*)trim(rname)//&
632      "call N2 cloud condensation ?"
633     condensn2=.false. ! default value
634     call getin_p("condensn2",condensn2)
635     if (is_master)write(*,*)trim(rname)//&
636      "condensn2 = ",condensn2
637
638     if (is_master)write(*,*)trim(rname)//&
639      "call no N2 frost formation?"
640     no_n2frost=.false. ! default value
641     call getin_p("no_n2frost",no_n2frost)
642     if (is_master)write(*,*)trim(rname)//&
643      "no_n2frost = ",no_n2frost
644
645     if (is_master)write(*,*)trim(rname)//&
646      "N2 condensation subtimestep?"
647     nbsub=20 ! default value
648     call getin_p("nbsub",nbsub)
649     if (is_master)write(*,*)trim(rname)//&
650      " nbsub = ",nbsub
651
652     if (is_master)write(*,*)trim(rname)//&
653      "Gravitationnal sedimentation ?"
654     sedimentation=.true. ! default value
655     call getin_p("sedimentation",sedimentation)
656     if (is_master)write(*,*)trim(rname)//&
657      " sedimentation = ",sedimentation
658
659     if (is_master)write(*,*)trim(rname)//&
660      "Compute glacial flow ?"
661     glaflow=.false. ! default value
662     call getin_p("glaflow",glaflow)
663     if (is_master)write(*,*)trim(rname)//&
664      " glaflow = ",glaflow
665
666     if (is_master)write(*,*)trim(rname)//&
667      "Compute methane cycle ?"
668     methane=.false. ! default value
669     call getin_p("methane",methane)
670     if (is_master)write(*,*)trim(rname)//&
671      " methane = ",methane
672     condmetsurf=.true. ! default value
673     call getin_p("condmetsurf",condmetsurf)
674     if (is_master)write(*,*)trim(rname)//&
675      " condmetsurf = ",condmetsurf
676     if (is_master)write(*,*)trim(rname)//&
677      "call no CH4 frost formation?"
678     no_ch4frost=.false. ! default value
679     call getin_p("no_ch4frost",no_ch4frost)
680     if (is_master)write(*,*)trim(rname)//&
681      "no_ch4frost = ",no_ch4frost
682
683     if (is_master)write(*,*)trim(rname)//&
684      "Compute methane clouds ?"
685     metcloud=.false. ! default value
686     call getin_p("metcloud",metcloud)
687     if (is_master)write(*,*)trim(rname)//&
688      " metcloud = ",metcloud
689
690     if (is_master)write(*,*)trim(rname)//&
691      "Compute CO cycle ?"
692     carbox=.false. ! default value
693     call getin_p("carbox",carbox)
694     if (is_master)write(*,*)trim(rname)//&
695      " carbox = ",carbox
696     condcosurf=.true. ! default value
697     call getin_p("condcosurf",condcosurf)
698     if (is_master)write(*,*)trim(rname)//&
699      " condcosurf = ",condcosurf
700
701     if (is_master)write(*,*)trim(rname)//&
702      "Compute CO clouds ?"
703     monoxcloud=.false. ! default value
704     call getin_p("monoxcloud",monoxcloud)
705     if (is_master)write(*,*)trim(rname)//&
706      " monoxcloud = ",monoxcloud
707
708     if (is_master)write(*,*)trim(rname)//&
709     "atmospheric redistribution (s):"
710     tau_n2=1. ! default value
711     call getin_p("tau_n2",tau_n2)
712     if (is_master)write(*,*)trim(rname)//&
713     " tau_n2 = ",tau_n2
714     tau_ch4=1.E7 ! default value
715     call getin_p("tau_ch4",tau_ch4)
716     if (is_master)write(*,*)trim(rname)//&
717     " tau_ch4 = ",tau_ch4
718     tau_co=1. ! default value
719     call getin_p("tau_co",tau_co)
720     if (is_master)write(*,*)trim(rname)//&
721     " tau_co = ",tau_co
722     tau_prechaze=1. ! default value
723     call getin_p("tau_prechaze",tau_prechaze)
724     if (is_master)write(*,*)trim(rname)//&
725     " tau_prechaze = ",tau_prechaze
726
727     if (is_master)write(*,*)trim(rname)//&
728      "Day fraction for limited cold trap in SP?"
729     dayfrac=0. ! default value
730     call getin_p("dayfrac",dayfrac)
731     if (is_master)write(*,*)trim(rname)//&
732      " dayfrac = ",dayfrac
733     thresh_non2=0. ! default value
734     call getin_p("thresh_non2",thresh_non2)
735     if (is_master)write(*,*)trim(rname)//&
736      " thresh_non2 = ",thresh_non2
737
738    ! use Pluto.old routines
739     if (is_master) write(*,*) trim(rname)//": use vdifc from old Pluto ?"
740     oldplutovdifc=.false. ! default value
741     call getin_p("oldplutovdifc",oldplutovdifc)
742     if (is_master) write(*,*) trim(rname)//": oldplutovdifc = ",oldplutovdifc
743
744     if (is_master) write(*,*) trim(rname)//&
745       ": call pluto.old correlated-k radiative transfer ?"
746     oldplutocorrk=.false. ! default value
747     call getin_p("oldplutocorrk",oldplutocorrk)
748     if (is_master) write(*,*) trim(rname)//": oldplutocorrk = ",oldplutocorrk
749
750     if (is_master) write(*,*) trim(rname)//&
751       ": call pluto.old sedimentation ?"
752     oldplutosedim=.false. ! default value
753     call getin_p("oldplutosedim",oldplutosedim)
754     if (is_master) write(*,*) trim(rname)//": oldplutosedim = ",oldplutosedim
755
756!***************************************************************
757     !! Haze options
758
759     ! Microphysical moment model
760     ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
761     if (is_master) write(*,*) "Run with or without microphysics?"
762     callmufi=.false. ! default value
763     call getin_p("callmufi",callmufi)
764     if (is_master) write(*,*)" callmufi = ",callmufi
765
766     if (is_master) write(*,*) "Run with or without microphysical clouds?"
767     callmuclouds=.false. ! default value
768     call getin_p("callmuclouds",callmuclouds)
769     if (is_master) write(*,*)" callmuclouds = ",callmuclouds
770
771     ! sanity check
772     if (callmufi.and.(.not.tracer)) then
773       print*,"You are running microphysics without tracer"
774       print*,"Please start again with tracer =.true."
775       stop
776     endif
777     if (callmuclouds.and.(.not.callmufi)) then
778       print*,"You are running microphysical clouds without microphysics"
779       print*,"Please start again with callmufi =.true."
780       stop
781     endif
782
783     if (is_master) write(*,*) "Path to microphysical config file?"
784     config_mufi='datagcm/microphysics/config.cfg' ! default value
785     call getin_p("config_mufi",config_mufi)
786     if (is_master) write(*,*) trim(rname)//" config_mufi = ",config_mufi
787
788     if (is_master) write(*,*) "Condensable gas fluxes datafile"
789     mugasflux_file='None' ! default file
790     call getin_p("mugasflux_file",mugasflux_file)
791     if (is_master) write(*,*) trim(rname)//" mugasflux_file = ",trim(mugasflux_file)
792
793     if (is_master) write(*,*) "Spherical aerosol optical properties datafile"
794     aersprop_file="optprop_rannou_r2-200nm_nu003.dat"  ! default file
795     call getin_p("aersprop_file",aersprop_file)
796     if (is_master) write(*,*) trim(rname)//" aersprop_file = ",trim(aersprop_file)
797
798     if (is_master) write(*,*) "Fractal aerosol optical properties datafile"
799     aerfprop_file="optprop_rannou_fractal_r010nm_N1_1e4_d2.dat"  ! default file
800     call getin_p("aerfprop_file",aerfprop_file)
801     if (is_master) write(*,*) trim(rname)//" aerfprop_file = ",trim(aerfprop_file)
802
803     if (is_master) write(*,*) "Use haze production from CH4 photolysis or production rate?"
804     call_haze_prod_pCH4=.false. ! default value
805     call getin_p("call_haze_prod_pCH4",call_haze_prod_pCH4)
806     if (is_master) write(*,*)" call_haze_prod_pCH4 = ",call_haze_prod_pCH4
807
808     if (is_master) write(*,*) "Pressure level of aerosols production (Pa)?"
809     haze_p_prod=1.0e-2 ! default value
810     call getin_p("haze_p_prod",haze_p_prod)
811     if (is_master) write(*,*)" haze_p_prod = ",haze_p_prod
812
813     if (is_master) write(*,*) "Aerosol production rate (kg.m-2.s-1)?"
814     haze_tx_prod=9.8e-14 ! default value
815     call getin_p("haze_tx_prod",haze_tx_prod)
816     if (is_master) write(*,*)" haze_tx_prod = ",haze_tx_prod
817
818     if (is_master) write(*,*) "Equivalent radius production (m)?"
819     haze_rc_prod=1.0e-9 ! default value
820     call getin_p("haze_rc_prod",haze_rc_prod)
821     if (is_master) write(*,*)" haze_rc_prod = ",haze_rc_prod
822
823     if (is_master) write(*,*) "Monomer radius (m)?"
824     haze_rm=1.0e-8 ! default value
825     call getin_p("haze_rm",haze_rm)
826     if (is_master) write(*,*)" haze_rm = ",haze_rm
827
828     if (is_master) write(*,*) "Aerosol's fractal dimension?"
829     haze_df=2.0 ! default value
830     call getin_p("haze_df",haze_df)
831     if (is_master) write(*,*)" haze_df = ",haze_df
832
833     if (is_master) write(*,*) "Aerosol density (kg.m-3)?"
834     haze_rho=800.0 ! default value
835     call getin_p("haze_rho",haze_rho)
836     if (is_master) write(*,*)" haze_rho = ",haze_rho
837
838     if (is_master) write(*,*) "Radius of air molecule (m)?"
839     air_rad=1.75e-10 ! default value
840     call getin_p("air_rad",air_rad)
841     if (is_master) write(*,*)" air_rad = ",air_rad
842
843     ! Conversion precursor to haze
844     ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
845     if (is_master) write(*,*) "Time constant of conversion in aerosol (s)?"
846     tcon_ch4=1.e7 ! default value
847     call getin_p("tcon_ch4",tcon_ch4)
848     if (is_master) write(*,*)" tcon_ch4 = ",tcon_ch4
849
850     if (is_master) write(*,*) "Ratio CH4 photolysis?"
851     k_ch4=1. ! default value
852     call getin_p("k_ch4",k_ch4)
853     if (is_master) write(*,*)" k_ch4 = ",k_ch4
854
855     if (is_master) write(*,*) "Nitrogen contribution ratio N/C?"
856     ncratio_ch4=0.5 ! default value
857     call getin_p("ncratio_ch4",ncratio_ch4)
858     if (is_master) write(*,*)" ncratio_ch4 = ",ncratio_ch4
859
860     ! Variables for aerosol absorption
861     ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
862     if (is_master) write(*,*) "Sph. aer. absorption correction in VI?"
863     Fabs_aers_VI=1. ! default value
864     call getin_p("Fabs_aers_VI",Fabs_aers_VI)
865     if (is_master) write(*,*)" Fabs_aers_VI = ",Fabs_aers_VI
866
867     if (is_master) write(*,*) "Fra. aer. absorption correction in VI?"
868     Fabs_aerf_VI=1. ! default value
869     call getin_p("Fabs_aerf_VI",Fabs_aerf_VI)
870     if (is_master) write(*,*)" Fabs_aerf_VI = ",Fabs_aerf_VI
871
872     if (is_master) write(*,*) "Sph. aer. absorption correction in IR?"
873     Fabs_aers_IR=1. ! default value
874     call getin_p("Fabs_aers_IR",Fabs_aers_IR)
875     if (is_master) write(*,*)" Fabs_aers_IR = ",Fabs_aers_IR
876
877     if (is_master) write(*,*) "Fra. aer. absorption correction in IR?"
878     Fabs_aerf_IR=1. ! default value
879     call getin_p("Fabs_aerf_IR",Fabs_aerf_IR)
880     if (is_master) write(*,*)" Fabs_aerf_IR = ",Fabs_aerf_IR
881
882     ! Pluto haze model
883     ! ~~~~~~~~~~~~~~~~
884     if (is_master)write(*,*)trim(rname)//&
885     "Production of haze ?"
886     haze=.false. ! default value
887     call getin_p("haze",haze)
888     if (is_master)write(*,*)trim(rname)//&
889     " haze = ",haze
890     hazeconservch4=.false. ! conservch4, by default value ch4 is photolyzed
891     call getin_p("hazeconservch4",hazeconservch4)
892     if (is_master)write(*,*)trim(rname)//&
893     "hazconservch4 = ",hazeconservch4
894     if (is_master) write(*,*) "Production of haze with CH4 fix profile?"
895     haze_ch4proffix=.false. ! default value
896     call getin_p("haze_ch4proffix",haze_ch4proffix)
897     if (is_master) write(*,*)" haze_ch4proffix = ",haze_ch4proffix
898     if (is_master)write(*,*)trim(rname)//&
899     "Production of haze (fast version) ?"
900     fasthaze=.false. ! default value
901     call getin_p("fasthaze",fasthaze)
902     if (is_master)write(*,*)trim(rname)//&
903     "fasthaze = ",fasthaze
904
905     if (is_master)write(*,*)trim(rname)//&
906     "Add source of haze ?"
907     source_haze=.false. ! default value
908     call getin_p("source_haze",source_haze)
909     if (is_master)write(*,*)trim(rname)//&
910     " source_haze = ",source_haze
911     mode_hs=0 ! mode haze source default value
912     call getin_p("mode_hs",mode_hs)
913     if (is_master)write(*,*)trim(rname)//&
914     " mode_hs",mode_hs
915     kfix=1 ! default value
916     call getin_p("kfix",kfix)
917     if (is_master)write(*,*)trim(rname)//&
918     "levels kfix",kfix
919     fracsource=0.1 ! default value
920     call getin_p("fracsource",fracsource)
921     if (is_master)write(*,*)trim(rname)//&
922     " fracsource",fracsource
923     latsource=30. ! default value
924     call getin_p("latsource",latsource)
925     if (is_master)write(*,*)trim(rname)//&
926     " latsource",latsource
927     lonsource=180. ! default value
928     call getin_p("lonsource",lonsource)
929     if (is_master)write(*,*)trim(rname)//&
930     " lonsource",lonsource
931
932     if (is_master)write(*,*)trim(rname)//&
933     "Radiatively active haze ?"
934     optichaze=.false. ! default value
935     call getin_p("optichaze",optichaze)
936     if (is_master)write(*,*)trim(rname)//&
937     "optichaze = ",optichaze
938
939     aerohaze=.false. ! default value
940     call getin_p("aerohaze",aerohaze)
941     if (aerohaze) then
942      if (is_master) write(*,*)trim(rname)//": aerohaze is deprecated.",&
943      "it is now called optichaze=.true."
944      call abort_physic(rname,"aerohaze is deprecated. It is now called optichaze",1)
945     endif
946
947     if (is_master)write(*,*)trim(rname)//&
948     "Haze monomer radius ?"
949     rad_haze=10.e-9 ! default value
950     call getin_p("rad_haze",rad_haze)
951     if (is_master)write(*,*)trim(rname)//&
952     "rad_haze = ",rad_haze
953
954     if (is_master)write(*,*)trim(rname)//&
955     "fractal particle ?"
956     fractal=.false. ! default value
957     call getin_p("fractal",fractal)
958     if (is_master)write(*,*)trim(rname)//&
959     "fractal = ",fractal
960     nb_monomer=10 ! default value
961     call getin_p("nb_monomer",nb_monomer)
962     if (is_master)write(*,*)trim(rname)//&
963     "nb_monomer = ",nb_monomer
964
965     if (is_master)write(*,*)trim(rname)//&
966     "fixed gaseous CH4 mixing ratio for rad transfer ?"
967    ch4fix=.false. ! default value
968    call getin_p("ch4fix",ch4fix)
969    if (is_master)write(*,*)trim(rname)//&
970     " ch4fix = ",ch4fix
971    if (is_master)write(*,*)trim(rname)//&
972     "fixed gaseous CH4 mixing ratio for rad transfer ?"
973    vmrch4fix=0.5 ! default value
974    call getin_p("vmrch4fix",vmrch4fix)
975    if (is_master)write(*,*)trim(rname)//&
976     " vmrch4fix = ",vmrch4fix
977    vmrch4_proffix=.false. ! default value
978    call getin_p("vmrch4_proffix",vmrch4_proffix)
979    if (is_master)write(*,*)trim(rname)//&
980     " vmrch4_proffix = ",vmrch4_proffix
981
982    if (is_master)write(*,*)trim(rname)//&
983     "call specific cooling for radiative transfer ?"
984    cooling=.false.  ! default value
985    call getin_p("cooling",cooling)
986    if (is_master)write(*,*)trim(rname)//&
987     " cooling = ",cooling
988
989    if (is_master)write(*,*)trim(rname)//&
990     "NLTE correction for methane heating rates?"
991    nlte=.false.  ! default value
992    call getin_p("nlte",nlte)
993    if (is_master)write(*,*)trim(rname)//&
994     " nlte = ",nlte
995    strobel=.true.  ! default value
996    call getin_p("strobel",strobel)
997    if (is_master)write(*,*)trim(rname)//&
998     " strobel = ",strobel
999
1000     if (is_master)write(*,*)trim(rname)//&
1001     "fixed radius profile from txt file ?"
1002     haze_radproffix=.false. ! default value
1003     call getin_p("haze_radproffix",haze_radproffix)
1004     if (is_master)write(*,*)trim(rname)//&
1005     "haze_radproffix = ",haze_radproffix
1006     if (is_master)write(*,*)trim(rname)//&
1007     "fixed MMR profile from txt file ?"
1008     haze_proffix=.false. ! default value
1009     call getin_p("haze_proffix",haze_proffix)
1010     if (is_master)write(*,*)trim(rname)//&
1011     "haze_proffix = ",haze_proffix
1012
1013     if (is_master)write(*,*)trim(rname)//&
1014     "Number mix ratio of haze particles for co clouds:"
1015     Nmix_co=100000. ! default value
1016     call getin_p("Nmix_co",Nmix_co)
1017     if (is_master)write(*,*)trim(rname)//&
1018     " Nmix_co = ",Nmix_co
1019
1020     if (is_master)write(*,*)trim(rname)//&
1021     "Number mix ratio of haze particles for ch4 clouds:"
1022     Nmix_ch4=100000. ! default value
1023     call getin_p("Nmix_ch4",Nmix_ch4)
1024     if (is_master)write(*,*)trim(rname)//&
1025     " Nmix_ch4 = ",Nmix_ch4
1026
1027!***************************************************************
1028     !! Surface properties
1029
1030!*********** N2 *********************************
1031
1032     if (is_master)write(*,*)trim(rname)//&
1033      "Mode for changing N2 albedo"
1034     mode_n2=0 ! default value
1035     call getin_p("mode_n2",mode_n2)
1036     if (is_master)write(*,*)trim(rname)//&
1037      " mode_n2 = ",mode_n2
1038     thres_n2ice=1. ! default value
1039     call getin_p("thres_n2ice",thres_n2ice)
1040     if (is_master)write(*,*)trim(rname)//&
1041      " thres_n2ice = ",thres_n2ice
1042
1043     if (is_master)write(*,*)trim(rname)//&
1044      "Diff of N2 albedo with thickness"
1045     deltab=0. ! default value
1046     call getin_p("deltab",deltab)
1047     if (is_master)write(*,*)trim(rname)//&
1048      " deltab = ",deltab
1049
1050     if (is_master)write(*,*)trim(rname)//&
1051      "albedo N2 beta "
1052     alb_n2b=0.7 ! default value
1053     call getin_p("alb_n2b",alb_n2b)
1054     if (is_master)write(*,*)trim(rname)//&
1055      " alb_n2b = ",alb_n2b
1056
1057     if (is_master)write(*,*)trim(rname)//&
1058      "albedo N2 alpha "
1059     alb_n2a=0.7 ! default value
1060     call getin_p("alb_n2a",alb_n2a)
1061     if (is_master)write(*,*)trim(rname)//&
1062      " alb_n2a = ",alb_n2a
1063
1064     if (is_master)write(*,*)trim(rname)//&
1065      "emis N2 beta "
1066     emis_n2b=0.7 ! default value
1067     call getin_p("emis_n2b",emis_n2b)
1068     if (is_master)write(*,*)trim(rname)//&
1069      " emis_n2b = ",emis_n2b
1070
1071     if (is_master)write(*,*)trim(rname)//&
1072      "emis N2 alpha "
1073     emis_n2a=0.7 ! default value
1074     call getin_p("emis_n2a",emis_n2a)
1075     if (is_master)write(*,*)trim(rname)//&
1076      " emis_n2a = ",emis_n2a
1077
1078!*********** CH4 *********************************
1079
1080     if (is_master)write(*,*)trim(rname)//&
1081      "Mode for changing CH4 albedo"
1082     mode_ch4=0 ! default value
1083     call getin_p("mode_ch4",mode_ch4)
1084     if (is_master)write(*,*)trim(rname)//&
1085      " mode_ch4 = ",mode_ch4
1086     feedback_met=0 ! default value
1087     call getin_p("feedback_met",feedback_met)
1088     if (is_master)write(*,*)trim(rname)//&
1089      " feedback_met = ",feedback_met
1090     thres_ch4ice=1. ! default value, mm
1091     call getin_p("thres_ch4ice",thres_ch4ice)
1092     if (is_master)write(*,*)trim(rname)//&
1093      " thres_ch4ice = ",thres_ch4ice
1094     fdch4_latn=-1.
1095     fdch4_lats=0.
1096     fdch4_lone=0.
1097     fdch4_lonw=-1.
1098     fdch4_depalb=0.5
1099     fdch4_finalb=0.95
1100     fdch4_maxalb=0.99
1101     fdch4_ampl=200.
1102     fdch4_maxice=100.
1103     call getin_p("fdch4_latn",fdch4_latn)
1104     call getin_p("fdch4_lats",fdch4_lats)
1105     call getin_p("fdch4_lone",fdch4_lone)
1106     call getin_p("fdch4_lonw",fdch4_lonw)
1107     call getin_p("fdch4_depalb",fdch4_depalb)
1108     call getin_p("fdch4_finalb",fdch4_finalb)
1109     call getin_p("fdch4_maxalb",fdch4_maxalb)
1110     call getin_p("fdch4_maxice",fdch4_maxice)
1111     call getin_p("fdch4_ampl",fdch4_ampl)
1112     if (is_master)write(*,*)trim(rname)//&
1113      " Values for albedo feedback = ",fdch4_latn,&
1114     fdch4_lats,fdch4_lone,fdch4_lonw,fdch4_depalb,&
1115     fdch4_finalb,fdch4_maxalb,fdch4_maxice,fdch4_ampl
1116
1117     if (is_master)write(*,*)trim(rname)//&
1118      "Latitude for diff albedo methane"
1119     metlateq=25. ! default value
1120     call getin_p("metlateq",metlateq)
1121     if (is_master)write(*,*)trim(rname)//&
1122      " metlateq = ",metlateq
1123
1124     if (is_master)write(*,*)trim(rname)//&
1125      "Ls1 and Ls2 for change of ch4 albedo"
1126     metls1=-1. ! default value
1127     metls2=-2. ! default value
1128     call getin_p("metls1",metls1)
1129     call getin_p("metls2",metls2)
1130
1131     if (is_master)write(*,*)trim(rname)//&
1132      "albedo CH4 "
1133     alb_ch4=0.5 ! default value
1134     call getin_p("alb_ch4",alb_ch4)
1135     if (is_master)write(*,*)trim(rname)//&
1136      " alb_ch4 = ",alb_ch4
1137
1138     if (is_master)write(*,*)trim(rname)//&
1139      "albedo equatorial CH4 "
1140     alb_ch4_eq=alb_ch4 ! default value
1141     call getin_p("alb_ch4_eq",alb_ch4_eq)
1142     if (is_master)write(*,*)trim(rname)//&
1143      " alb_ch4_eq = ",alb_ch4_eq
1144
1145     if (is_master)write(*,*)trim(rname)//&
1146      "albedo south hemis CH4 "
1147     alb_ch4_s=alb_ch4 ! default value
1148     call getin_p("alb_ch4_s",alb_ch4_s)
1149     if (is_master)write(*,*)trim(rname)//&
1150      " alb_ch4_s = ",alb_ch4_s
1151
1152     if (is_master)write(*,*)trim(rname)//&
1153      "emis CH4 "
1154     emis_ch4=0.5 ! default value
1155     call getin_p("emis_ch4",emis_ch4)
1156     if (is_master)write(*,*)trim(rname)//&
1157      " emis_ch4 = ",emis_ch4
1158
1159     if (is_master)write(*,*)trim(rname)//&
1160      "CH4 lag for n2 sublimation limitation"
1161     ch4lag=.false. ! default value
1162     latlag=-90. ! default value
1163     vmrlag=1. ! default value
1164     call getin_p("ch4lag",ch4lag)
1165     call getin_p("latlag",latlag)
1166     call getin_p("vmrlag",vmrlag)
1167     if (is_master)write(*,*)trim(rname)//&
1168      " ch4lag = ",ch4lag
1169     if (is_master)write(*,*)trim(rname)//&
1170      " latlag = ",latlag
1171     if (is_master)write(*,*)trim(rname)//&
1172      " vmrlag = ",vmrlag
1173
1174     if (is_master)write(*,*)trim(rname)//&
1175      "Max temperature for surface ?"
1176     tsurfmax=.false. ! default value
1177     albmin_ch4=0.3 ! default value
1178     call getin_p("tsurfmax",tsurfmax)
1179     call getin_p("albmin_ch4",albmin_ch4)
1180     if (is_master)write(*,*)trim(rname)//&
1181      " tsurfmax = ",tsurfmax
1182     if (is_master)write(*,*)trim(rname)//&
1183      " albmin_ch4 = ",albmin_ch4
1184
1185
1186     if (is_master)write(*,*)trim(rname)//&
1187     "fixed gaseous CH4 mixing ratio for rad transfer ?"
1188     ch4fix=.false. ! default value
1189     call getin_p("ch4fix",ch4fix)
1190     if (is_master)write(*,*)trim(rname)//&
1191       " ch4fix = ",ch4fix
1192     if (is_master)write(*,*)trim(rname)//&
1193       "fixed gaseous CH4 mixing ratio for rad transfer ?"
1194     vmrch4fix=0.5 ! default value
1195     call getin_p("vmrch4fix",vmrch4fix)
1196     if (is_master)write(*,*)trim(rname)//&
1197        " vmrch4fix = ",vmrch4fix
1198     vmrch4_proffix=.false. ! default value
1199     call getin_p("vmrch4_proffix",vmrch4_proffix)
1200     if (is_master)write(*,*)trim(rname)//&
1201        " vmrch4_proffix = ",vmrch4_proffix
1202
1203
1204!*********** CO *********************************
1205
1206     if (is_master)write(*,*)trim(rname)//&
1207      "albedo CO "
1208     alb_co=0.4 ! default value
1209     call getin_p("alb_co",alb_co)
1210     if (is_master)write(*,*)trim(rname)//&
1211      " alb_co = ",alb_co
1212     thres_coice=1. ! default value, mm
1213     call getin_p("thres_coice",thres_coice)
1214     if (is_master)write(*,*)trim(rname)//&
1215      " thres_coice = ",thres_coice
1216
1217     if (is_master)write(*,*)trim(rname)//&
1218      "emis CO "
1219     emis_co=0.4 ! default value
1220     call getin_p("emis_co",emis_co)
1221     if (is_master)write(*,*)trim(rname)//&
1222      " emis_co = ",emis_co
1223
1224!*********** THOLINS *********************************
1225     if (is_master)write(*,*)trim(rname)//&
1226      "Mode for changing tholins albedo/emis"
1227     mode_tholins=0 ! default value
1228     call getin_p("mode_tholins",mode_tholins)
1229     if (is_master)write(*,*)trim(rname)//&
1230      " mode_tholins = ",mode_tholins
1231
1232     if (is_master)write(*,*)trim(rname)//&
1233      "albedo tho "
1234     alb_tho=0.1 ! default value
1235     call getin_p("alb_tho",alb_tho)
1236     if (is_master)write(*,*)trim(rname)//&
1237      " alb_tho = ",alb_tho
1238
1239     if (is_master)write(*,*)trim(rname)//&
1240      "albedo tho eq"
1241     alb_tho_eq=0.1 ! default value
1242     call getin_p("alb_tho_eq",alb_tho_eq)
1243     if (is_master)write(*,*)trim(rname)//&
1244      " alb_tho_eq = ",alb_tho_eq
1245
1246     if (is_master)write(*,*)trim(rname)//&
1247      "emis tholins "
1248     emis_tho=1. ! default value
1249     call getin_p("emis_tho",emis_tho)
1250     if (is_master)write(*,*)trim(rname)//&
1251      " emis_tho = ",emis_tho
1252
1253     if (is_master)write(*,*)trim(rname)//&
1254      "emis tholins eq"
1255     emis_tho_eq=1. ! default value
1256     call getin_p("emis_tho_eq",emis_tho_eq)
1257     if (is_master)write(*,*)trim(rname)//&
1258      " emis_tho_eq = ",emis_tho_eq
1259
1260     if (is_master)write(*,*)trim(rname)//&
1261      "Latitude for diff albedo tholins"
1262     tholateq=25. ! default value
1263     call getin_p("tholateq",tholateq)
1264     if (is_master)write(*,*)trim(rname)//&
1265      " tholateq = ",tholateq
1266     tholatn=-1.
1267     tholats=0.
1268     tholone=0.
1269     tholonw=-1.
1270     alb_tho_spe=0.1 ! default value
1271     emis_tho_spe=1. ! default value
1272     call getin_p("  tholatn",tholatn)
1273     call getin_p("  tholats",tholats)
1274     call getin_p("  tholone",tholone)
1275     call getin_p("  tholonw",tholonw)
1276     if (is_master)write(*,*)trim(rname)//&
1277      " Values for special tholins albedo = ",tholatn,&
1278        tholats,tholone,tholonw,alb_tho_spe,emis_tho_spe
1279
1280     if (is_master)write(*,*)trim(rname)//&
1281      "Specific albedo"
1282     spelon1=-180. ! default value
1283     spelon2=180. ! default value
1284     spelat1=-90. ! default value
1285     spelat2=90. ! default value
1286     specalb=.false. ! default value
1287     if (is_master)write(*,*)trim(rname)//&
1288      "albedo/emis spe "
1289     albspe=0.1 ! default value
1290     emispe=1. ! default value
1291     call getin_p("spelon1",spelon1)
1292     call getin_p("spelon2",spelon2)
1293     call getin_p("spelat1",spelat1)
1294     call getin_p("spelat2",spelat2)
1295     call getin_p("specalb",specalb)
1296     call getin_p("albspe",albspe)
1297     call getin_p("emispe",emispe)
1298
1299     if (is_master)write(*,*)trim(rname)//&
1300      " specific = ",specalb
1301
1302!********************** TI *****************************
1303
1304     if (is_master)write(*,*)trim(rname)//&
1305      "Change TI with time"
1306     changeti=.false. ! default value
1307     call getin_p("changeti",changeti)
1308     if (is_master)write(*,*)trim(rname)//&
1309      " changeti = ",changeti
1310     changetid=.false. ! default value for diurnal TI
1311     call getin_p("changetid",changetid)
1312     if (is_master)write(*,*)trim(rname)//&
1313      " changetid = ",changetid
1314
1315     if (is_master)write(*,*)trim(rname)//&
1316      "IT N2 "
1317     ITN2=800. ! default value
1318     call getin_p("ITN2",ITN2)
1319     if (is_master)write(*,*)trim(rname)//&
1320      " ITN2 = ",ITN2
1321     ITN2d=20. ! default value
1322     call getin_p("ITN2d",ITN2d)
1323     if (is_master)write(*,*)trim(rname)//&
1324      " ITN2d = ",ITN2d
1325
1326     if (is_master)write(*,*)trim(rname)//&
1327      "IT CH4"
1328     ITCH4=800. ! default value
1329     call getin_p("ITCH4",ITCH4)
1330     if (is_master)write(*,*)trim(rname)//&
1331      " ITCH4 = ",ITCH4
1332     ITCH4d=20. ! default value
1333     call getin_p("ITCH4d",ITCH4d)
1334     if (is_master)write(*,*)trim(rname)//&
1335      " ITCH4d = ",ITCH4d
1336
1337     if (is_master)write(*,*)trim(rname)//&
1338      "IT H2O"
1339     ITH2O=800. ! default value
1340     call getin_p("ITH2O",ITH2O)
1341     if (is_master)write(*,*)trim(rname)//&
1342      " ITH2O = ",ITH2O
1343     ITH2Od=20. ! default value
1344     call getin_p("ITH2Od",ITH2Od)
1345     if (is_master)write(*,*)trim(rname)//&
1346      " ITH2Od = ",ITH2Od
1347
1348!************** COOLING ***************
1349
1350     alpha_top=5e-11 ! default value
1351     call getin_p("alpha_top",alpha_top)
1352     if (is_master)write(*,*)trim(rname)//&
1353      " alpha_top = ",alpha_top
1354     pref=0.02 ! default value
1355     call getin_p("pref",pref)
1356     if (is_master)write(*,*)trim(rname)//&
1357      " pref = ",pref
1358     deltap=0.1 ! default value
1359     call getin_p("deltap",deltap)
1360     if (is_master)write(*,*)trim(rname)//&
1361      " deltap = ",deltap
1362
1363!=================================
1364
1365     if (is_master) write(*,*)trim(rname)//&
1366       ": Is the variable gas species radiatively active?"
1367     Tstrat=167.0
1368     varactive=.false.
1369     call getin_p("varactive",varactive)
1370     if (is_master) write(*,*)trim(rname)//": varactive = ",varactive
1371
1372     if (is_master) write(*,*)trim(rname)//&
1373       ": Is the variable gas species distribution set?"
1374     varfixed=.false.
1375     call getin_p("varfixed",varfixed)
1376     if (is_master) write(*,*)trim(rname)//": varfixed = ",varfixed
1377
1378! Test of incompatibility:
1379! if varactive, then varfixed should be false
1380     if (varactive.and.varfixed) then
1381       call abort_physic(rname,'if varactive, varfixed must be OFF!',1)
1382     endif
1383
1384     if (is_master) write(*,*)trim(rname)//": Gravitationnal sedimentation ?"
1385     sedimentation=.false. ! default value
1386     call getin_p("sedimentation",sedimentation)
1387     if (is_master) write(*,*)trim(rname)//": sedimentation = ",sedimentation
1388
1389     if (is_master) write(*,*)trim(rname)//": Spectral Dependant albedo ?"
1390     albedo_spectral_mode=.false. ! default value
1391     call getin_p("albedo_spectral_mode",albedo_spectral_mode)
1392     if (is_master) write(*,*)trim(rname)//&
1393     ": albedo_spectral_mode = ",albedo_spectral_mode
1394
1395     if (is_master) then
1396       write(*,*)trim(rname)//": Snow albedo ?"
1397       write(*,*)trim(rname)//": If albedo_spectral_mode=.true., then this "
1398       write(*,*)trim(rname)//": corresponds to the 0.5 microns snow albedo."
1399     endif
1400     albedosnow=0.5         ! default value
1401     call getin_p("albedosnow",albedosnow)
1402     if (is_master) write(*,*)trim(rname)//": albedosnow = ",albedosnow
1403
1404     if (is_master) write(*,*)trim(rname)//": N2 ice albedo ?"
1405     albedon2ice=0.5       ! default value
1406     call getin_p("albedon2ice",albedon2ice)
1407     if (is_master) write(*,*)trim(rname)//": albedon2ice = ",albedon2ice
1408
1409     if (is_master) write(*,*)trim(rname)//": Maximum ice thickness ?"
1410     maxicethick=2.0         ! default value
1411     call getin_p("maxicethick",maxicethick)
1412     if (is_master) write(*,*)trim(rname)//": maxicethick = ",maxicethick
1413
1414     if (is_master) write(*,*)trim(rname)//": Minimum eddy mix coeff in 1D ?"
1415     kmixmin=1.0e-2         ! default value
1416     call getin_p("kmixmin",kmixmin)
1417     if (is_master) write(*,*)trim(rname)//": kmixmin = ",kmixmin
1418     kmix_proffix=.false.  ! default value
1419     call getin_p("kmix_proffix",kmix_proffix)
1420     if (is_master) write(*,*)trim(rname)//": kmix_proffix = ",kmix_proffix
1421
1422     if (is_master) write(*,*)'Predefined Cp from dynamics is ',cpp,'J kg^-1 K^-1'
1423     if (is_master) write(*,*)'Predefined Mg from dynamics is ',mugaz,'amu'
1424
1425     force_cpp=.false. ! default value
1426     call getin_p("force_cpp",force_cpp)
1427     if (force_cpp) then
1428      if (is_master) write(*,*)trim(rname)//": force_cpp = ",force_cpp
1429      if (is_master) write(*,*)trim(rname)//": force_cpp is deprecated.",&
1430      "Set cpp_mugaz_mode=1 in callfis to emulate force_cpp=.true."
1431      call abort_physic(rname,"Anyway, you need to set force_cpp=.false. to continue.",1)
1432     endif
1433
1434     if (is_master) write(*,*)trim(rname)//&
1435     ": where do you want your cpp/mugaz value to come from?",&
1436     "=> 0: dynamics (3d), 1: forced in callfis (1d), 2: computed from gases.def (1d)?"
1437     cpp_mugaz_mode = 0 ! default value
1438     call getin_p("cpp_mugaz_mode",cpp_mugaz_mode)
1439     if (is_master) write(*,*)trim(rname)//": cpp_mugaz_mode = ",cpp_mugaz_mode
1440
1441
1442! Test of incompatibility:
1443
1444     if ((.not.tracer).and.(haze)) then
1445       call abort_physic(rname, 'if haze are on, tracers must be on!', 1)
1446     endif
1447     if ((callmufi).and.(haze)) then
1448       call abort_physic(rname, 'if haze are on, microphysics should be deactivated!', 1)
1449     endif
1450     if ((haze).and.(naerkind.gt.1)) then
1451      call abort_physic(rname, 'if haze are on, naerkind must be equal to 1!', 1)
1452     endif
1453     if ((callmufi).and..not.(naerkind.gt.1)) then
1454      call abort_physic(rname, 'if microphysics is on, naerkind must be > 1!', 1)
1455     endif
1456     if ((callmufi).and.(callmuclouds).and..not.(naerkind.gt.2)) then
1457      call abort_physic(rname, 'if microphysical clouds are on, naerkind must be > 2!', 1)
1458     endif
1459     if (.not.(callmufi.or.haze).and.(optichaze)) then
1460      call abort_physic(rname, 'if microphysics and haze are off, optichaze must be deactivated!', 1)
1461     endif
1462     if ((callmufi.and.call_haze_prod_pCH4).and..not.(methane)) then
1463      call abort_physic(rname, 'if haze production from CH4 photolysis is on, methane must be activated!', 1)
1464     endif
1465     if (haze_proffix.and.sedimentation) then
1466         call abort_physic(rname, 'if haze profile is set, sedimentation must be deactivated', 1)
1467     endif
1468     if (callmolvis.and..not.callconduct) then
1469         call abort_physic(rname, 'if callmolvis is set, callconduct must be true', 1)
1470     endif
1471     if (glaflow.and..not.fast) then
1472         call abort_physic(rname, 'if glaflow is set, fast must be true', 1)
1473     endif
1474     if (paleo.and..not.fast) then
1475         call abort_physic(rname, 'if paleo is set, fast must be true', 1)
1476     endif
1477     if ((haze_proffix.or.haze_radproffix).and..not.optichaze) then
1478         call abort_physic(rname, 'for now, haze/rad_proffix only works with optichaze=T', 1)
1479     endif
1480     ! if (carbox.and.condcosurf.and.no_n2frost) then
1481     !   call abort_physic(rname, "CO surface condensation and no_n2frost are both active which may not be relevant", 1)
1482     ! end if
1483
1484     if ((cpp_mugaz_mode >= 1).and.(is_master).and.(ngrid>1)) then
1485        write(*,*)'    !!!!! Be aware that having different values for cpp and mugaz in the dynamics and physics'
1486        write(*,*)'    in 3D can result in very pathological behavior. You have been warned !!!!!'
1487      else if ((cpp_mugaz_mode == 0).and.(is_master).and.(ngrid==1)) then
1488        ! for this specific 1D error we will remove run.def before aborting JL22
1489        call system("rm -rf run.def")
1490        call abort_physic(rname,"cpp_mugaz_mode must be >= 1 in 1d",1)
1491     endif
1492
1493     if (cpp_mugaz_mode == 1) then
1494       mugaz = -99999.
1495       if (is_master) write(*,*)trim(rname)//&
1496         ": MEAN MOLECULAR MASS in g mol-1 ?"
1497       call getin_p("mugaz",mugaz)
1498       IF (mugaz.eq.-99999.) THEN
1499         call abort_physic(rname,"mugaz must be set if cpp_mugaz_mode = 1",1)
1500       ENDIF
1501       cpp = -99999.
1502       if (is_master) write(*,*)trim(rname)//&
1503         ": SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
1504       call getin_p("cpp",cpp)
1505       IF (cpp.eq.-99999.) THEN
1506           PRINT *, "cpp must be set if cpp_mugaz_mode = 1"
1507           STOP
1508       ENDIF
1509       if (is_master) write(*,*)'New Cp from callfis is ',cpp,'J kg^-1 K^-1'
1510       if (is_master) write(*,*)'New Mg from callfis is ',mugaz,'amu'
1511
1512     endif ! of if (cpp_mugaz_mode == 1)
1513     call su_gases
1514     call calc_cpp_mugaz
1515
1516     if (is_master) then
1517       PRINT*,'--------------------------------------------'
1518       PRINT*
1519       PRINT*
1520     endif
1521  ELSE
1522     call abort_physic(rname,'Cannot read file callphys.def. Is it here ?',1)
1523  ENDIF ! of IF(iscallphys)
1524
1525  if (is_master) then
1526    PRINT*
1527    PRINT*,'inifis: daysec',daysec
1528    PRINT*
1529    PRINT*,'inifis: The radiative transfer is computed:'
1530    PRINT*,'           each ',iradia,' physical time-step'
1531    PRINT*,'        or each ',iradia*dtphys,' seconds'
1532    PRINT*
1533  endif
1534
1535!-----------------------------------------------------------------------
1536!     Some more initialization:
1537!     ------------------------
1538
1539  ! Initializations for comgeomfi_h
1540#ifndef MESOSCALE
1541  totarea=SSUM(ngrid,parea,1)
1542  call planetwide_sumval(parea,totarea_planet)
1543
1544  !! those are defined in comdiurn_h.F90
1545  IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
1546  IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
1547  IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
1548  IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
1549
1550  DO ig=1,ngrid
1551     sinlat(ig)=sin(plat(ig))
1552     coslat(ig)=cos(plat(ig))
1553     sinlon(ig)=sin(plon(ig))
1554     coslon(ig)=cos(plon(ig))
1555  ENDDO
1556#endif
1557  ! initialize variables in radinc_h
1558  call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
1559
1560  ! initialize variables and allocate arrays in radcommon_h
1561  call ini_radcommon_h(naerkind)
1562
1563  ! allocate "comsoil_h" arrays
1564  call ini_comsoil_h(ngrid)
1565
1566  ! allocate arrays in "nonoro_gwd_ran_mod"
1567  call end_nonoro_gwd_ran
1568  call ini_nonoro_gwd_ran(ngrid,nlayer)
1569
1570  ! allocate arrays in "nonoro_gwd_mix_mod"
1571  call end_nonoro_gwd_mix
1572  call ini_nonoro_gwd_mix(ngrid,nlayer,nq)
1573
1574  END SUBROUTINE inifis
1575
1576END MODULE inifis_mod
Note: See TracBrowser for help on using the repository browser.