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

Last change on this file since 3973 was 3973, checked in by tbertrand, 4 weeks ago

PLUTO PCM:
Small fix when adapting N2 ice albedo to match surface pressure in the case of Triton (convention : positive declination in 1989)
TB

File size: 58.5 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, clddprop_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(*,*) "Cloud drop optical properties datafile"
804     clddprop_file="optprop_rannou_r2-200nm_nu003.dat"  ! default file
805     call getin_p("clddprop_file",clddprop_file)
806     if (is_master) write(*,*) trim(rname)//" clddprop_file = ",trim(clddprop_file)
807
808     if (is_master) write(*,*) "Use haze production from CH4 photolysis or production rate?"
809     call_haze_prod_pCH4=.false. ! default value
810     call getin_p("call_haze_prod_pCH4",call_haze_prod_pCH4)
811     if (is_master) write(*,*)" call_haze_prod_pCH4 = ",call_haze_prod_pCH4
812
813     if (is_master) write(*,*) "Pressure level of aerosols production (Pa)?"
814     haze_p_prod=1.0e-2 ! default value
815     call getin_p("haze_p_prod",haze_p_prod)
816     if (is_master) write(*,*)" haze_p_prod = ",haze_p_prod
817
818     if (is_master) write(*,*) "Aerosol production rate (kg.m-2.s-1)?"
819     haze_tx_prod=9.8e-14 ! default value
820     call getin_p("haze_tx_prod",haze_tx_prod)
821     if (is_master) write(*,*)" haze_tx_prod = ",haze_tx_prod
822
823     if (is_master) write(*,*) "Equivalent radius production (m)?"
824     haze_rc_prod=1.0e-9 ! default value
825     call getin_p("haze_rc_prod",haze_rc_prod)
826     if (is_master) write(*,*)" haze_rc_prod = ",haze_rc_prod
827
828     if (is_master) write(*,*) "Monomer radius (m)?"
829     haze_rm=1.0e-8 ! default value
830     call getin_p("haze_rm",haze_rm)
831     if (is_master) write(*,*)" haze_rm = ",haze_rm
832
833     if (is_master) write(*,*) "Aerosol's fractal dimension?"
834     haze_df=2.0 ! default value
835     call getin_p("haze_df",haze_df)
836     if (is_master) write(*,*)" haze_df = ",haze_df
837
838     if (is_master) write(*,*) "Aerosol density (kg.m-3)?"
839     haze_rho=800.0 ! default value
840     call getin_p("haze_rho",haze_rho)
841     if (is_master) write(*,*)" haze_rho = ",haze_rho
842
843     if (is_master) write(*,*) "Radius of air molecule (m)?"
844     air_rad=1.75e-10 ! default value
845     call getin_p("air_rad",air_rad)
846     if (is_master) write(*,*)" air_rad = ",air_rad
847
848     ! Conversion precursor to haze
849     ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
850     if (is_master) write(*,*) "Time constant of conversion in aerosol (s)?"
851     tcon_ch4=1.e7 ! default value
852     call getin_p("tcon_ch4",tcon_ch4)
853     if (is_master) write(*,*)" tcon_ch4 = ",tcon_ch4
854
855     if (is_master) write(*,*) "Ratio CH4 photolysis?"
856     k_ch4=1. ! default value
857     call getin_p("k_ch4",k_ch4)
858     if (is_master) write(*,*)" k_ch4 = ",k_ch4
859
860     if (is_master) write(*,*) "Nitrogen contribution ratio N/C?"
861     ncratio_ch4=0.5 ! default value
862     call getin_p("ncratio_ch4",ncratio_ch4)
863     if (is_master) write(*,*)" ncratio_ch4 = ",ncratio_ch4
864
865     ! Variables for aerosol absorption
866     ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
867     if (is_master) write(*,*) "Sph. aer. absorption correction in VI?"
868     Fabs_aers_VI=1. ! default value
869     call getin_p("Fabs_aers_VI",Fabs_aers_VI)
870     if (is_master) write(*,*)" Fabs_aers_VI = ",Fabs_aers_VI
871
872     if (is_master) write(*,*) "Fra. aer. absorption correction in VI?"
873     Fabs_aerf_VI=1. ! default value
874     call getin_p("Fabs_aerf_VI",Fabs_aerf_VI)
875     if (is_master) write(*,*)" Fabs_aerf_VI = ",Fabs_aerf_VI
876
877     if (is_master) write(*,*) "Sph. aer. absorption correction in IR?"
878     Fabs_aers_IR=1. ! default value
879     call getin_p("Fabs_aers_IR",Fabs_aers_IR)
880     if (is_master) write(*,*)" Fabs_aers_IR = ",Fabs_aers_IR
881
882     if (is_master) write(*,*) "Fra. aer. absorption correction in IR?"
883     Fabs_aerf_IR=1. ! default value
884     call getin_p("Fabs_aerf_IR",Fabs_aerf_IR)
885     if (is_master) write(*,*)" Fabs_aerf_IR = ",Fabs_aerf_IR
886
887     if (is_master) write(*,*) "Cloud drop absorption correction in VI?"
888     Fabs_cldd_VI=1. ! default value
889     call getin_p("Fabs_cldd_VI",Fabs_cldd_VI)
890     if (is_master) write(*,*)" Fabs_cldd_VI = ",Fabs_cldd_VI
891
892     if (is_master) write(*,*) "Cloud drop absorption correction in IR?"
893     Fabs_cldd_IR=1. ! default value
894     call getin_p("Fabs_cldd_IR",Fabs_cldd_IR)
895     if (is_master) write(*,*)" Fabs_cldd_IR = ",Fabs_cldd_IR
896
897     ! Pluto haze model
898     ! ~~~~~~~~~~~~~~~~
899     if (is_master)write(*,*)trim(rname)//&
900     "Production of haze ?"
901     haze=.false. ! default value
902     call getin_p("haze",haze)
903     if (is_master)write(*,*)trim(rname)//&
904     " haze = ",haze
905     hazeconservch4=.false. ! conservch4, by default value ch4 is photolyzed
906     call getin_p("hazeconservch4",hazeconservch4)
907     if (is_master)write(*,*)trim(rname)//&
908     "hazconservch4 = ",hazeconservch4
909     if (is_master) write(*,*) "Production of haze with CH4 fix profile?"
910     haze_ch4proffix=.false. ! default value
911     call getin_p("haze_ch4proffix",haze_ch4proffix)
912     if (is_master) write(*,*)" haze_ch4proffix = ",haze_ch4proffix
913     if (is_master)write(*,*)trim(rname)//&
914     "Production of haze (fast version) ?"
915     fasthaze=.false. ! default value
916     call getin_p("fasthaze",fasthaze)
917     if (is_master)write(*,*)trim(rname)//&
918     "fasthaze = ",fasthaze
919
920     if (is_master)write(*,*)trim(rname)//&
921     "Add source of haze ?"
922     source_haze=.false. ! default value
923     call getin_p("source_haze",source_haze)
924     if (is_master)write(*,*)trim(rname)//&
925     " source_haze = ",source_haze
926     mode_hs=0 ! mode haze source default value
927     call getin_p("mode_hs",mode_hs)
928     if (is_master)write(*,*)trim(rname)//&
929     " mode_hs",mode_hs
930     kfix=1 ! default value
931     call getin_p("kfix",kfix)
932     if (is_master)write(*,*)trim(rname)//&
933     "levels kfix",kfix
934     fracsource=0.1 ! default value
935     call getin_p("fracsource",fracsource)
936     if (is_master)write(*,*)trim(rname)//&
937     " fracsource",fracsource
938     latsource=30. ! default value
939     call getin_p("latsource",latsource)
940     if (is_master)write(*,*)trim(rname)//&
941     " latsource",latsource
942     lonsource=180. ! default value
943     call getin_p("lonsource",lonsource)
944     if (is_master)write(*,*)trim(rname)//&
945     " lonsource",lonsource
946
947     if (is_master)write(*,*)trim(rname)//&
948     "Radiatively active haze ?"
949     optichaze=.false. ! default value
950     call getin_p("optichaze",optichaze)
951     if (is_master)write(*,*)trim(rname)//&
952     "optichaze = ",optichaze
953
954     aerohaze=.false. ! default value
955     call getin_p("aerohaze",aerohaze)
956     if (aerohaze) then
957      if (is_master) write(*,*)trim(rname)//": aerohaze is deprecated.",&
958      "it is now called optichaze=.true."
959      call abort_physic(rname,"aerohaze is deprecated. It is now called optichaze",1)
960     endif
961
962     if (is_master)write(*,*)trim(rname)//&
963     "Haze monomer radius ?"
964     rad_haze=10.e-9 ! default value
965     call getin_p("rad_haze",rad_haze)
966     if (is_master)write(*,*)trim(rname)//&
967     "rad_haze = ",rad_haze
968
969     if (is_master)write(*,*)trim(rname)//&
970     "fractal particle ?"
971     fractal=.false. ! default value
972     call getin_p("fractal",fractal)
973     if (is_master)write(*,*)trim(rname)//&
974     "fractal = ",fractal
975     nb_monomer=10 ! default value
976     call getin_p("nb_monomer",nb_monomer)
977     if (is_master)write(*,*)trim(rname)//&
978     "nb_monomer = ",nb_monomer
979
980     if (is_master)write(*,*)trim(rname)//&
981     "fixed gaseous CH4 mixing ratio for rad transfer ?"
982    ch4fix=.false. ! default value
983    call getin_p("ch4fix",ch4fix)
984    if (is_master)write(*,*)trim(rname)//&
985     " ch4fix = ",ch4fix
986    if (is_master)write(*,*)trim(rname)//&
987     "fixed gaseous CH4 mixing ratio for rad transfer ?"
988    vmrch4fix=0.5 ! default value
989    call getin_p("vmrch4fix",vmrch4fix)
990    if (is_master)write(*,*)trim(rname)//&
991     " vmrch4fix = ",vmrch4fix
992    vmrch4_proffix=.false. ! default value
993    call getin_p("vmrch4_proffix",vmrch4_proffix)
994    if (is_master)write(*,*)trim(rname)//&
995     " vmrch4_proffix = ",vmrch4_proffix
996
997    if (is_master)write(*,*)trim(rname)//&
998     "call specific cooling for radiative transfer ?"
999    cooling=.false.  ! default value
1000    call getin_p("cooling",cooling)
1001    if (is_master)write(*,*)trim(rname)//&
1002     " cooling = ",cooling
1003
1004    if (is_master)write(*,*)trim(rname)//&
1005     "NLTE correction for methane heating rates?"
1006    nlte=.false.  ! default value
1007    call getin_p("nlte",nlte)
1008    if (is_master)write(*,*)trim(rname)//&
1009     " nlte = ",nlte
1010    strobel=.true.  ! default value
1011    call getin_p("strobel",strobel)
1012    if (is_master)write(*,*)trim(rname)//&
1013     " strobel = ",strobel
1014
1015     if (is_master)write(*,*)trim(rname)//&
1016     "fixed radius profile from txt file ?"
1017     haze_radproffix=.false. ! default value
1018     call getin_p("haze_radproffix",haze_radproffix)
1019     if (is_master)write(*,*)trim(rname)//&
1020     "haze_radproffix = ",haze_radproffix
1021     if (is_master)write(*,*)trim(rname)//&
1022     "fixed MMR profile from txt file ?"
1023     haze_proffix=.false. ! default value
1024     call getin_p("haze_proffix",haze_proffix)
1025     if (is_master)write(*,*)trim(rname)//&
1026     "haze_proffix = ",haze_proffix
1027
1028     if (is_master)write(*,*)trim(rname)//&
1029     "Number mix ratio of haze particles for co clouds:"
1030     Nmix_co=100000. ! default value
1031     call getin_p("Nmix_co",Nmix_co)
1032     if (is_master)write(*,*)trim(rname)//&
1033     " Nmix_co = ",Nmix_co
1034
1035     if (is_master)write(*,*)trim(rname)//&
1036     "Number mix ratio of haze particles for ch4 clouds:"
1037     Nmix_ch4=100000. ! default value
1038     call getin_p("Nmix_ch4",Nmix_ch4)
1039     if (is_master)write(*,*)trim(rname)//&
1040     " Nmix_ch4 = ",Nmix_ch4
1041
1042!***************************************************************
1043     !! Surface properties
1044
1045!*********** N2 *********************************
1046
1047     if (is_master)write(*,*)trim(rname)//&
1048      "Mode for changing N2 albedo"
1049     mode_n2=0 ! default value
1050     call getin_p("mode_n2",mode_n2)
1051     if (is_master)write(*,*)trim(rname)//&
1052      " mode_n2 = ",mode_n2
1053     thres_n2ice=1. ! default value
1054     call getin_p("thres_n2ice",thres_n2ice)
1055     if (is_master)write(*,*)trim(rname)//&
1056      " thres_n2ice = ",thres_n2ice
1057
1058     if (is_master)write(*,*)trim(rname)//&
1059      "Diff of N2 albedo with thickness"
1060     deltab=0. ! default value
1061     call getin_p("deltab",deltab)
1062     if (is_master)write(*,*)trim(rname)//&
1063      " deltab = ",deltab
1064
1065     if (is_master)write(*,*)trim(rname)//&
1066      "albedo N2 beta "
1067     alb_n2b=0.7 ! default value
1068     call getin_p("alb_n2b",alb_n2b)
1069     if (is_master)write(*,*)trim(rname)//&
1070      " alb_n2b = ",alb_n2b
1071
1072     if (is_master)write(*,*)trim(rname)//&
1073      "albedo N2 alpha "
1074     alb_n2a=0.7 ! default value
1075     call getin_p("alb_n2a",alb_n2a)
1076     if (is_master)write(*,*)trim(rname)//&
1077      " alb_n2a = ",alb_n2a
1078
1079     if (is_master)write(*,*)trim(rname)//&
1080      "emis N2 beta "
1081     emis_n2b=0.7 ! default value
1082     call getin_p("emis_n2b",emis_n2b)
1083     if (is_master)write(*,*)trim(rname)//&
1084      " emis_n2b = ",emis_n2b
1085
1086     if (is_master)write(*,*)trim(rname)//&
1087      "emis N2 alpha "
1088     emis_n2a=0.7 ! default value
1089     call getin_p("emis_n2a",emis_n2a)
1090     if (is_master)write(*,*)trim(rname)//&
1091      " emis_n2a = ",emis_n2a
1092
1093!*********** CH4 *********************************
1094
1095     if (is_master)write(*,*)trim(rname)//&
1096      "Mode for changing CH4 albedo"
1097     mode_ch4=0 ! default value
1098     call getin_p("mode_ch4",mode_ch4)
1099     if (is_master)write(*,*)trim(rname)//&
1100      " mode_ch4 = ",mode_ch4
1101     feedback_met=0 ! default value
1102     call getin_p("feedback_met",feedback_met)
1103     if (is_master)write(*,*)trim(rname)//&
1104      " feedback_met = ",feedback_met
1105     thres_ch4ice=1. ! default value, mm
1106     call getin_p("thres_ch4ice",thres_ch4ice)
1107     if (is_master)write(*,*)trim(rname)//&
1108      " thres_ch4ice = ",thres_ch4ice
1109     fdch4_latn=-1.
1110     fdch4_lats=0.
1111     fdch4_lone=0.
1112     fdch4_lonw=-1.
1113     fdch4_depalb=0.5
1114     fdch4_finalb=0.95
1115     fdch4_maxalb=0.99
1116     fdch4_ampl=200.
1117     fdch4_maxice=100.
1118     call getin_p("fdch4_latn",fdch4_latn)
1119     call getin_p("fdch4_lats",fdch4_lats)
1120     call getin_p("fdch4_lone",fdch4_lone)
1121     call getin_p("fdch4_lonw",fdch4_lonw)
1122     call getin_p("fdch4_depalb",fdch4_depalb)
1123     call getin_p("fdch4_finalb",fdch4_finalb)
1124     call getin_p("fdch4_maxalb",fdch4_maxalb)
1125     call getin_p("fdch4_maxice",fdch4_maxice)
1126     call getin_p("fdch4_ampl",fdch4_ampl)
1127     if (is_master)write(*,*)trim(rname)//&
1128      " Values for albedo feedback = ",fdch4_latn,&
1129     fdch4_lats,fdch4_lone,fdch4_lonw,fdch4_depalb,&
1130     fdch4_finalb,fdch4_maxalb,fdch4_maxice,fdch4_ampl
1131
1132     if (is_master)write(*,*)trim(rname)//&
1133      "Latitude for diff albedo methane"
1134     metlateq=25. ! default value
1135     call getin_p("metlateq",metlateq)
1136     if (is_master)write(*,*)trim(rname)//&
1137      " metlateq = ",metlateq
1138
1139     if (is_master)write(*,*)trim(rname)//&
1140      "Ls1 and Ls2 for change of ch4 albedo"
1141     metls1=-1. ! default value
1142     metls2=-2. ! default value
1143     call getin_p("metls1",metls1)
1144     call getin_p("metls2",metls2)
1145
1146     if (is_master)write(*,*)trim(rname)//&
1147      "albedo CH4 "
1148     alb_ch4=0.5 ! default value
1149     call getin_p("alb_ch4",alb_ch4)
1150     if (is_master)write(*,*)trim(rname)//&
1151      " alb_ch4 = ",alb_ch4
1152
1153     if (is_master)write(*,*)trim(rname)//&
1154      "albedo equatorial CH4 "
1155     alb_ch4_eq=alb_ch4 ! default value
1156     call getin_p("alb_ch4_eq",alb_ch4_eq)
1157     if (is_master)write(*,*)trim(rname)//&
1158      " alb_ch4_eq = ",alb_ch4_eq
1159
1160     if (is_master)write(*,*)trim(rname)//&
1161      "albedo south hemis CH4 "
1162     alb_ch4_s=alb_ch4 ! default value
1163     call getin_p("alb_ch4_s",alb_ch4_s)
1164     if (is_master)write(*,*)trim(rname)//&
1165      " alb_ch4_s = ",alb_ch4_s
1166
1167     if (is_master)write(*,*)trim(rname)//&
1168      "emis CH4 "
1169     emis_ch4=0.5 ! default value
1170     call getin_p("emis_ch4",emis_ch4)
1171     if (is_master)write(*,*)trim(rname)//&
1172      " emis_ch4 = ",emis_ch4
1173
1174     if (is_master)write(*,*)trim(rname)//&
1175      "CH4 lag for n2 sublimation limitation"
1176     ch4lag=.false. ! default value
1177     latlag=-90. ! default value
1178     vmrlag=1. ! default value
1179     call getin_p("ch4lag",ch4lag)
1180     call getin_p("latlag",latlag)
1181     call getin_p("vmrlag",vmrlag)
1182     if (is_master)write(*,*)trim(rname)//&
1183      " ch4lag = ",ch4lag
1184     if (is_master)write(*,*)trim(rname)//&
1185      " latlag = ",latlag
1186     if (is_master)write(*,*)trim(rname)//&
1187      " vmrlag = ",vmrlag
1188
1189     if (is_master)write(*,*)trim(rname)//&
1190      "Max temperature for surface ?"
1191     tsurfmax=.false. ! default value
1192     albmin_ch4=0.3 ! default value
1193     call getin_p("tsurfmax",tsurfmax)
1194     call getin_p("albmin_ch4",albmin_ch4)
1195     if (is_master)write(*,*)trim(rname)//&
1196      " tsurfmax = ",tsurfmax
1197     if (is_master)write(*,*)trim(rname)//&
1198      " albmin_ch4 = ",albmin_ch4
1199
1200
1201     if (is_master)write(*,*)trim(rname)//&
1202     "fixed gaseous CH4 mixing ratio for rad transfer ?"
1203     ch4fix=.false. ! default value
1204     call getin_p("ch4fix",ch4fix)
1205     if (is_master)write(*,*)trim(rname)//&
1206       " ch4fix = ",ch4fix
1207     if (is_master)write(*,*)trim(rname)//&
1208       "fixed gaseous CH4 mixing ratio for rad transfer ?"
1209     vmrch4fix=0.5 ! default value
1210     call getin_p("vmrch4fix",vmrch4fix)
1211     if (is_master)write(*,*)trim(rname)//&
1212        " vmrch4fix = ",vmrch4fix
1213     vmrch4_proffix=.false. ! default value
1214     call getin_p("vmrch4_proffix",vmrch4_proffix)
1215     if (is_master)write(*,*)trim(rname)//&
1216        " vmrch4_proffix = ",vmrch4_proffix
1217
1218
1219!*********** CO *********************************
1220
1221     if (is_master)write(*,*)trim(rname)//&
1222      "albedo CO "
1223     alb_co=0.4 ! default value
1224     call getin_p("alb_co",alb_co)
1225     if (is_master)write(*,*)trim(rname)//&
1226      " alb_co = ",alb_co
1227     thres_coice=1. ! default value, mm
1228     call getin_p("thres_coice",thres_coice)
1229     if (is_master)write(*,*)trim(rname)//&
1230      " thres_coice = ",thres_coice
1231
1232     if (is_master)write(*,*)trim(rname)//&
1233      "emis CO "
1234     emis_co=0.4 ! default value
1235     call getin_p("emis_co",emis_co)
1236     if (is_master)write(*,*)trim(rname)//&
1237      " emis_co = ",emis_co
1238
1239!*********** THOLINS *********************************
1240     if (is_master)write(*,*)trim(rname)//&
1241      "Mode for changing tholins albedo/emis"
1242     mode_tholins=0 ! default value
1243     call getin_p("mode_tholins",mode_tholins)
1244     if (is_master)write(*,*)trim(rname)//&
1245      " mode_tholins = ",mode_tholins
1246
1247     if (is_master)write(*,*)trim(rname)//&
1248      "albedo tho "
1249     alb_tho=0.1 ! default value
1250     call getin_p("alb_tho",alb_tho)
1251     if (is_master)write(*,*)trim(rname)//&
1252      " alb_tho = ",alb_tho
1253
1254     if (is_master)write(*,*)trim(rname)//&
1255      "albedo tho eq"
1256     alb_tho_eq=0.1 ! default value
1257     call getin_p("alb_tho_eq",alb_tho_eq)
1258     if (is_master)write(*,*)trim(rname)//&
1259      " alb_tho_eq = ",alb_tho_eq
1260
1261     if (is_master)write(*,*)trim(rname)//&
1262      "emis tholins "
1263     emis_tho=1. ! default value
1264     call getin_p("emis_tho",emis_tho)
1265     if (is_master)write(*,*)trim(rname)//&
1266      " emis_tho = ",emis_tho
1267
1268     if (is_master)write(*,*)trim(rname)//&
1269      "emis tholins eq"
1270     emis_tho_eq=1. ! default value
1271     call getin_p("emis_tho_eq",emis_tho_eq)
1272     if (is_master)write(*,*)trim(rname)//&
1273      " emis_tho_eq = ",emis_tho_eq
1274
1275     if (is_master)write(*,*)trim(rname)//&
1276      "Latitude for diff albedo tholins"
1277     tholateq=25. ! default value
1278     call getin_p("tholateq",tholateq)
1279     if (is_master)write(*,*)trim(rname)//&
1280      " tholateq = ",tholateq
1281     tholatn=-1.
1282     tholats=0.
1283     tholone=0.
1284     tholonw=-1.
1285     alb_tho_spe=0.1 ! default value
1286     emis_tho_spe=1. ! default value
1287     call getin_p("  tholatn",tholatn)
1288     call getin_p("  tholats",tholats)
1289     call getin_p("  tholone",tholone)
1290     call getin_p("  tholonw",tholonw)
1291     if (is_master)write(*,*)trim(rname)//&
1292      " Values for special tholins albedo = ",tholatn,&
1293        tholats,tholone,tholonw,alb_tho_spe,emis_tho_spe
1294
1295     if (is_master)write(*,*)trim(rname)//&
1296      "Specific albedo"
1297     spelon1=-180. ! default value
1298     spelon2=180. ! default value
1299     spelat1=-90. ! default value
1300     spelat2=90. ! default value
1301     specalb=.false. ! default value
1302     if (is_master)write(*,*)trim(rname)//&
1303      "albedo/emis spe "
1304     albspe=0.1 ! default value
1305     emispe=1. ! default value
1306     call getin_p("spelon1",spelon1)
1307     call getin_p("spelon2",spelon2)
1308     call getin_p("spelat1",spelat1)
1309     call getin_p("spelat2",spelat2)
1310     call getin_p("specalb",specalb)
1311     call getin_p("albspe",albspe)
1312     call getin_p("emispe",emispe)
1313
1314     if (is_master)write(*,*)trim(rname)//&
1315      " specific = ",specalb
1316
1317!********************** TI *****************************
1318
1319     if (is_master)write(*,*)trim(rname)//&
1320      "Change TI with time"
1321     changeti=.false. ! default value
1322     call getin_p("changeti",changeti)
1323     if (is_master)write(*,*)trim(rname)//&
1324      " changeti = ",changeti
1325     changetid=.false. ! default value for diurnal TI
1326     call getin_p("changetid",changetid)
1327     if (is_master)write(*,*)trim(rname)//&
1328      " changetid = ",changetid
1329
1330     if (is_master)write(*,*)trim(rname)//&
1331      "IT N2 "
1332     ITN2=800. ! default value
1333     call getin_p("ITN2",ITN2)
1334     if (is_master)write(*,*)trim(rname)//&
1335      " ITN2 = ",ITN2
1336     ITN2d=20. ! default value
1337     call getin_p("ITN2d",ITN2d)
1338     if (is_master)write(*,*)trim(rname)//&
1339      " ITN2d = ",ITN2d
1340
1341     if (is_master)write(*,*)trim(rname)//&
1342      "IT CH4"
1343     ITCH4=800. ! default value
1344     call getin_p("ITCH4",ITCH4)
1345     if (is_master)write(*,*)trim(rname)//&
1346      " ITCH4 = ",ITCH4
1347     ITCH4d=20. ! default value
1348     call getin_p("ITCH4d",ITCH4d)
1349     if (is_master)write(*,*)trim(rname)//&
1350      " ITCH4d = ",ITCH4d
1351
1352     if (is_master)write(*,*)trim(rname)//&
1353      "IT H2O"
1354     ITH2O=800. ! default value
1355     call getin_p("ITH2O",ITH2O)
1356     if (is_master)write(*,*)trim(rname)//&
1357      " ITH2O = ",ITH2O
1358     ITH2Od=20. ! default value
1359     call getin_p("ITH2Od",ITH2Od)
1360     if (is_master)write(*,*)trim(rname)//&
1361      " ITH2Od = ",ITH2Od
1362
1363!************** COOLING ***************
1364
1365     alpha_top=5e-11 ! default value
1366     call getin_p("alpha_top",alpha_top)
1367     if (is_master)write(*,*)trim(rname)//&
1368      " alpha_top = ",alpha_top
1369     pref=0.02 ! default value
1370     call getin_p("pref",pref)
1371     if (is_master)write(*,*)trim(rname)//&
1372      " pref = ",pref
1373     deltap=0.1 ! default value
1374     call getin_p("deltap",deltap)
1375     if (is_master)write(*,*)trim(rname)//&
1376      " deltap = ",deltap
1377
1378!=================================
1379
1380     if (is_master) write(*,*)trim(rname)//&
1381       ": Is the variable gas species radiatively active?"
1382     Tstrat=167.0
1383     varactive=.false.
1384     call getin_p("varactive",varactive)
1385     if (is_master) write(*,*)trim(rname)//": varactive = ",varactive
1386
1387     if (is_master) write(*,*)trim(rname)//&
1388       ": Is the variable gas species distribution set?"
1389     varfixed=.false.
1390     call getin_p("varfixed",varfixed)
1391     if (is_master) write(*,*)trim(rname)//": varfixed = ",varfixed
1392
1393! Test of incompatibility:
1394! if varactive, then varfixed should be false
1395     if (varactive.and.varfixed) then
1396       call abort_physic(rname,'if varactive, varfixed must be OFF!',1)
1397     endif
1398
1399     if (is_master) write(*,*)trim(rname)//": Gravitationnal sedimentation ?"
1400     sedimentation=.false. ! default value
1401     call getin_p("sedimentation",sedimentation)
1402     if (is_master) write(*,*)trim(rname)//": sedimentation = ",sedimentation
1403
1404     if (is_master) write(*,*)trim(rname)//": Spectral Dependant albedo ?"
1405     albedo_spectral_mode=.false. ! default value
1406     call getin_p("albedo_spectral_mode",albedo_spectral_mode)
1407     if (is_master) write(*,*)trim(rname)//&
1408     ": albedo_spectral_mode = ",albedo_spectral_mode
1409
1410     if (is_master) then
1411       write(*,*)trim(rname)//": Snow albedo ?"
1412       write(*,*)trim(rname)//": If albedo_spectral_mode=.true., then this "
1413       write(*,*)trim(rname)//": corresponds to the 0.5 microns snow albedo."
1414     endif
1415     albedosnow=0.5         ! default value
1416     call getin_p("albedosnow",albedosnow)
1417     if (is_master) write(*,*)trim(rname)//": albedosnow = ",albedosnow
1418
1419     if (is_master) write(*,*)trim(rname)//": N2 ice albedo ?"
1420     albedon2ice=0.5       ! default value
1421     call getin_p("albedon2ice",albedon2ice)
1422     if (is_master) write(*,*)trim(rname)//": albedon2ice = ",albedon2ice
1423
1424     if (is_master) write(*,*)trim(rname)//": Maximum ice thickness ?"
1425     maxicethick=2.0         ! default value
1426     call getin_p("maxicethick",maxicethick)
1427     if (is_master) write(*,*)trim(rname)//": maxicethick = ",maxicethick
1428
1429     if (is_master) write(*,*)trim(rname)//": Minimum eddy mix coeff in 1D ?"
1430     kmixmin=1.0e-2         ! default value
1431     call getin_p("kmixmin",kmixmin)
1432     if (is_master) write(*,*)trim(rname)//": kmixmin = ",kmixmin
1433     kmix_proffix=.false.  ! default value
1434     call getin_p("kmix_proffix",kmix_proffix)
1435     if (is_master) write(*,*)trim(rname)//": kmix_proffix = ",kmix_proffix
1436
1437     if (is_master) write(*,*)'Predefined Cp from dynamics is ',cpp,'J kg^-1 K^-1'
1438     if (is_master) write(*,*)'Predefined Mg from dynamics is ',mugaz,'amu'
1439
1440     force_cpp=.false. ! default value
1441     call getin_p("force_cpp",force_cpp)
1442     if (force_cpp) then
1443      if (is_master) write(*,*)trim(rname)//": force_cpp = ",force_cpp
1444      if (is_master) write(*,*)trim(rname)//": force_cpp is deprecated.",&
1445      "Set cpp_mugaz_mode=1 in callfis to emulate force_cpp=.true."
1446      call abort_physic(rname,"Anyway, you need to set force_cpp=.false. to continue.",1)
1447     endif
1448
1449     if (is_master) write(*,*)trim(rname)//&
1450     ": where do you want your cpp/mugaz value to come from?",&
1451     "=> 0: dynamics (3d), 1: forced in callfis (1d), 2: computed from gases.def (1d)?"
1452     cpp_mugaz_mode = 0 ! default value
1453     call getin_p("cpp_mugaz_mode",cpp_mugaz_mode)
1454     if (is_master) write(*,*)trim(rname)//": cpp_mugaz_mode = ",cpp_mugaz_mode
1455
1456
1457! Test of incompatibility:
1458
1459     if ((.not.tracer).and.(haze)) then
1460       call abort_physic(rname, 'if haze are on, tracers must be on!', 1)
1461     endif
1462     if ((callmufi).and.(haze)) then
1463       call abort_physic(rname, 'if haze are on, microphysics should be deactivated!', 1)
1464     endif
1465     if ((haze).and.(naerkind.gt.1)) then
1466      call abort_physic(rname, 'if haze are on, naerkind must be equal to 1!', 1)
1467     endif
1468     if ((callmufi).and..not.(naerkind.gt.1)) then
1469      call abort_physic(rname, 'if microphysics is on, naerkind must be > 1!', 1)
1470     endif
1471     if ((callmufi).and..not.(callmuclouds).and.(naerkind.gt.2)) then
1472      call abort_physic(rname, 'Warning: here microphysical clouds are on, naerkind must be = 2!', 1)
1473     endif
1474     if ((callmufi).and.(callmuclouds).and..not.(naerkind.gt.2)) then
1475      call abort_physic(rname, 'if microphysical clouds are on, naerkind must be > 2!', 1)
1476     endif
1477     if (.not.(callmufi.or.haze).and.(optichaze)) then
1478      call abort_physic(rname, 'if microphysics and haze are off, optichaze must be deactivated!', 1)
1479     endif
1480     if ((callmufi.and.call_haze_prod_pCH4).and..not.(methane)) then
1481      call abort_physic(rname, 'if haze production from CH4 photolysis is on, methane must be activated!', 1)
1482     endif
1483     if (haze_proffix.and.sedimentation) then
1484         call abort_physic(rname, 'if haze profile is set, sedimentation must be deactivated', 1)
1485     endif
1486     if (callmolvis.and..not.callconduct) then
1487         call abort_physic(rname, 'if callmolvis is set, callconduct must be true', 1)
1488     endif
1489     if (glaflow.and..not.fast) then
1490         call abort_physic(rname, 'if glaflow is set, fast must be true', 1)
1491     endif
1492     if (paleo.and..not.fast) then
1493         call abort_physic(rname, 'if paleo is set, fast must be true', 1)
1494     endif
1495     if ((haze_proffix.or.haze_radproffix).and..not.optichaze) then
1496         call abort_physic(rname, 'for now, haze/rad_proffix only works with optichaze=T', 1)
1497     endif
1498     ! if (carbox.and.condcosurf.and.no_n2frost) then
1499     !   call abort_physic(rname, "CO surface condensation and no_n2frost are both active which may not be relevant", 1)
1500     ! end if
1501
1502     if ((cpp_mugaz_mode >= 1).and.(is_master).and.(ngrid>1)) then
1503        write(*,*)'    !!!!! Be aware that having different values for cpp and mugaz in the dynamics and physics'
1504        write(*,*)'    in 3D can result in very pathological behavior. You have been warned !!!!!'
1505      else if ((cpp_mugaz_mode == 0).and.(is_master).and.(ngrid==1)) then
1506        ! for this specific 1D error we will remove run.def before aborting JL22
1507        call system("rm -rf run.def")
1508        call abort_physic(rname,"cpp_mugaz_mode must be >= 1 in 1d",1)
1509     endif
1510
1511     if (cpp_mugaz_mode == 1) then
1512       mugaz = -99999.
1513       if (is_master) write(*,*)trim(rname)//&
1514         ": MEAN MOLECULAR MASS in g mol-1 ?"
1515       call getin_p("mugaz",mugaz)
1516       IF (mugaz.eq.-99999.) THEN
1517         call abort_physic(rname,"mugaz must be set if cpp_mugaz_mode = 1",1)
1518       ENDIF
1519       cpp = -99999.
1520       if (is_master) write(*,*)trim(rname)//&
1521         ": SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
1522       call getin_p("cpp",cpp)
1523       IF (cpp.eq.-99999.) THEN
1524           PRINT *, "cpp must be set if cpp_mugaz_mode = 1"
1525           STOP
1526       ENDIF
1527       if (is_master) write(*,*)'New Cp from callfis is ',cpp,'J kg^-1 K^-1'
1528       if (is_master) write(*,*)'New Mg from callfis is ',mugaz,'amu'
1529
1530     endif ! of if (cpp_mugaz_mode == 1)
1531     call su_gases
1532     call calc_cpp_mugaz
1533
1534     if (is_master) then
1535       PRINT*,'--------------------------------------------'
1536       PRINT*
1537       PRINT*
1538     endif
1539  ELSE
1540     call abort_physic(rname,'Cannot read file callphys.def. Is it here ?',1)
1541  ENDIF ! of IF(iscallphys)
1542
1543  if (is_master) then
1544    PRINT*
1545    PRINT*,'inifis: daysec',daysec
1546    PRINT*
1547    PRINT*,'inifis: The radiative transfer is computed:'
1548    PRINT*,'           each ',iradia,' physical time-step'
1549    PRINT*,'        or each ',iradia*dtphys,' seconds'
1550    PRINT*
1551  endif
1552
1553!-----------------------------------------------------------------------
1554!     Some more initialization:
1555!     ------------------------
1556
1557  ! Initializations for comgeomfi_h
1558#ifndef MESOSCALE
1559  totarea=SSUM(ngrid,parea,1)
1560  call planetwide_sumval(parea,totarea_planet)
1561
1562  !! those are defined in comdiurn_h.F90
1563  IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
1564  IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
1565  IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
1566  IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
1567
1568  DO ig=1,ngrid
1569     sinlat(ig)=sin(plat(ig))
1570     coslat(ig)=cos(plat(ig))
1571     sinlon(ig)=sin(plon(ig))
1572     coslon(ig)=cos(plon(ig))
1573  ENDDO
1574#endif
1575  ! initialize variables in radinc_h
1576  call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
1577
1578  ! initialize variables and allocate arrays in radcommon_h
1579  call ini_radcommon_h(naerkind)
1580
1581  ! allocate "comsoil_h" arrays
1582  call ini_comsoil_h(ngrid)
1583
1584  ! allocate arrays in "nonoro_gwd_ran_mod"
1585  call end_nonoro_gwd_ran
1586  call ini_nonoro_gwd_ran(ngrid,nlayer)
1587
1588  ! allocate arrays in "nonoro_gwd_mix_mod"
1589  call end_nonoro_gwd_mix
1590  call ini_nonoro_gwd_mix(ngrid,nlayer,nq)
1591
1592  END SUBROUTINE inifis
1593
1594END MODULE inifis_mod
Note: See TracBrowser for help on using the repository browser.