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

Last change on this file since 3945 was 3940, checked in by tbertrand, 8 weeks ago

PLUTO PCM :
Preparing for Triton simulations, implementation of internal heat flux and new North-South convention for the subsolar point
TB

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