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

Last change on this file since 3640 was 3627, checked in by tbertrand, 11 months ago

Pluto PCM: a fix for ASR, some cleanups and additional output
TB

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