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

Last change on this file since 3585 was 3585, checked in by debatzbr, 12 days ago

Connecting microphysics to radiative transfer + miscellaneous cleans

File size: 51.1 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
315     if (is_master) write(*,*) trim(rname)//": use turbdiff instead of vdifc ?"
316     UseTurbDiff=.false. ! default value
317     call getin_p("UseTurbDiff",UseTurbDiff)
318     if (is_master) write(*,*) trim(rname)//": UseTurbDiff = ",UseTurbDiff
319
320     if (is_master) write(*,*) trim(rname)//": call convective adjustment ?"
321     calladj=.false. ! default value
322     call getin_p("calladj",calladj)
323     if (is_master) write(*,*) trim(rname)//": calladj = ",calladj
324
325     if (is_master) write(*,*) trim(rname)//&
326     ": Radiative timescale for Newtonian cooling (in Earth days)?"
327     tau_relax=30. ! default value
328     call getin_p("tau_relax",tau_relax)
329     if (is_master) write(*,*) trim(rname)//": tau_relax = ",tau_relax
330     tau_relax=tau_relax*24*3600 ! convert Earth days --> seconds
331
332     if (is_master) write(*,*)trim(rname)//&
333       ": call thermal conduction in the soil ?"
334     callsoil=.true. ! default value
335     call getin_p("callsoil",callsoil)
336     if (is_master) write(*,*) trim(rname)//": callsoil = ",callsoil
337
338     if (is_master) write(*,*)trim(rname)//&
339       ": Rad transfer is computed every iradia", &
340       " physical timestep"
341     iradia=1 ! default value
342     call getin_p("iradia",iradia)
343     if (is_master) write(*,*)trim(rname)//": iradia = ",iradia
344
345     if (is_master) write(*,*)trim(rname)//": Rayleigh scattering ?"
346     rayleigh=.false.
347     call getin_p("rayleigh",rayleigh)
348     if (is_master) write(*,*)trim(rname)//": rayleigh = ",rayleigh
349
350     if (is_master) write(*,*) trim(rname)//&
351       ": Use blackbody for stellar spectrum ?"
352     stelbbody=.false. ! default value
353     call getin_p("stelbbody",stelbbody)
354     if (is_master) write(*,*) trim(rname)//": stelbbody = ",stelbbody
355
356     if (is_master) write(*,*) trim(rname)//": Stellar blackbody temperature ?"
357     stelTbb=5800.0 ! default value
358     call getin_p("stelTbb",stelTbb)
359     if (is_master) write(*,*) trim(rname)//": stelTbb = ",stelTbb
360
361     if (is_master) write(*,*) trim(rname)//": Output mean OLR in 1D?"
362     meanOLR=.false.
363     call getin_p("meanOLR",meanOLR)
364     if (is_master) write(*,*) trim(rname)//": meanOLR = ",meanOLR
365
366     if (is_master) write(*,*)trim(rname)//": Output spectral OLR in 3D?"
367     specOLR=.false.
368     call getin_p("specOLR",specOLR)
369     if (is_master) write(*,*)trim(rname)//": specOLR = ",specOLR
370
371     if (is_master) write(*,*)trim(rname)//&
372       ": Output diagnostic optical thickness attenuation exp(-tau)?"
373     diagdtau=.false.
374     call getin_p("diagdtau",diagdtau)
375     if (is_master) write(*,*)trim(rname)//": diagdtau = ",diagdtau
376
377     if (is_master) write(*,*)trim(rname)//": Operate in kastprof mode?"
378     kastprof=.false.
379     call getin_p("kastprof",kastprof)
380     if (is_master) write(*,*)trim(rname)//": kastprof = ",kastprof
381
382     if (is_master) write(*,*)trim(rname)//&
383       ": Uniform absorption in radiative transfer?"
384     graybody=.false.
385     call getin_p("graybody",graybody)
386     if (is_master) write(*,*)trim(rname)//": graybody = ",graybody
387
388! Soil model
389     if (is_master) write(*,*)trim(rname)//&
390       ": Number of sub-surface layers for soil scheme?"
391     ! default value of nsoilmx set in comsoil_h
392     call getin_p("nsoilmx",nsoilmx)
393     if (is_master) write(*,*)trim(rname)//": nsoilmx=",nsoilmx
394
395     if (is_master) write(*,*)trim(rname)//&
396       ": Thickness of topmost soil layer (m)?"
397     ! default value of lay1_soil set in comsoil_h
398     call getin_p("lay1_soil",lay1_soil)
399     if (is_master) write(*,*)trim(rname)//": lay1_soil = ",lay1_soil
400
401     if (is_master) write(*,*)trim(rname)//&
402       ": Coefficient for soil layer thickness distribution?"
403     ! default value of alpha_soil set in comsoil_h
404     call getin_p("alpha_soil",alpha_soil)
405     if (is_master) write(*,*)trim(rname)//": alpha_soil = ",alpha_soil
406
407! Chemistry in the thermosphere
408     if (is_master) write(*,*) trim(rname)//": Use deposition ?"
409     depos=.false.         ! default value
410     call getin_p("depos",depos)
411     if (is_master) write(*,*) trim(rname)//": depos = ",depos
412
413     if (is_master) write(*,*)trim(rname)//": Production of haze ?"
414     haze=.false. ! default value
415     call getin_p("haze",haze)
416     if (is_master) write(*,*)trim(rname)//": haze = ",haze
417
418
419
420      if (is_master) write(*,*)trim(rname)// "call thermal conduction ?"
421      callconduct=.false. ! default value
422      call getin_p("callconduct",callconduct)
423      if (is_master) write(*,*)trim(rname)// " callconduct = ",callconduct
424
425      if (is_master) write(*,*)trim(rname)// "call phitop ?"
426      phitop=0. ! default value
427      call getin_p("phitop",phitop)
428      if (is_master) write(*,*)trim(rname)// " phitop = ",phitop
429
430      if (is_master) write(*,*)trim(rname)// "call molecular viscosity ?"
431      callmolvis=.false. ! default value
432      call getin_p("callmolvis",callmolvis)
433      if (is_master) write(*,*)trim(rname)// " callmolvis = ",callmolvis
434
435      if (is_master) write(*,*)trim(rname)// "call molecular diffusion ?"
436      callmoldiff=.false. ! default value
437      call getin_p("callmoldiff",callmoldiff)
438      if (is_master) write(*,*)trim(rname)// " callmoldiff = ",callmoldiff
439
440! Global1D mean and solar zenith angle
441
442     if (ngrid.eq.1) then
443      PRINT*, 'Simulate global averaged conditions ?'
444      global1d = .false. ! default value
445      call getin_p("global1d",global1d)
446      write(*,*) "global1d = ",global1d
447
448      ! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
449      if (global1d.and.diurnal) then
450         call abort_physic(rname,'if global1d is true, diurnal must be set to false',1)
451      endif
452
453      if (global1d) then
454         PRINT *,'Solar Zenith angle (deg.) ?'
455         PRINT *,'(assumed for averaged solar flux S/4)'
456         szangle=60.0  ! default value
457         call getin_p("szangle",szangle)
458         write(*,*) "szangle = ",szangle
459      endif
460   endif ! of if (ngrid.eq.1)
461
462! Test of incompatibility:
463! if kastprof used, we must be in 1D
464     if (kastprof.and.(ngrid.gt.1)) then
465       call abort_physic(rname,'kastprof can only be used in 1D!',1)
466     endif
467
468     if (is_master) write(*,*)trim(rname)//&
469       ": Stratospheric temperature for kastprof mode?"
470     Tstrat=167.0
471     call getin_p("Tstrat",Tstrat)
472     if (is_master) write(*,*)trim(rname)//": Tstrat = ",Tstrat
473
474     if (is_master) write(*,*)trim(rname)//": Remove lower boundary?"
475     nosurf=.false.
476     call getin_p("nosurf",nosurf)
477     if (is_master) write(*,*)trim(rname)//": nosurf = ",nosurf
478
479! Tests of incompatibility:
480     if (nosurf.and.callsoil) then
481       if (is_master) then
482         write(*,*)trim(rname)//'nosurf not compatible with soil scheme!'
483         write(*,*)trim(rname)//'... got to make a choice!'
484       endif
485       call abort_physic(rname,"nosurf not compatible with soil scheme!",1)
486     endif
487
488     if (is_master) write(*,*)trim(rname)//": Add an internal heat flux?", &
489                   "... matters only if callsoil=F"
490     intheat=0.
491     call getin_p("intheat",intheat)
492     if (is_master) write(*,*)trim(rname)//": intheat = ",intheat
493
494
495     if (is_master) write(*,*)trim(rname)//": Test physics timescale in 1D?"
496     testradtimes=.false.
497     call getin_p("testradtimes",testradtimes)
498     if (is_master) write(*,*)trim(rname)//": testradtimes = ",testradtimes
499
500! Test of incompatibility:
501! if testradtimes used, we must be in 1D
502     if (testradtimes.and.(ngrid.gt.1)) then
503       call abort_physic(rname,'testradtimes can only be used in 1D!',1)
504     endif
505
506     if (is_master) write(*,*)trim(rname)//": Default planetary temperature?"
507     tplanet=215.0
508     call getin_p("tplanet",tplanet)
509     if (is_master) write(*,*)trim(rname)//": tplanet = ",tplanet
510
511     if (is_master) write(*,*)trim(rname)//": Which star?"
512     startype=1 ! default value = Sol
513     call getin_p("startype",startype)
514     if (is_master) write(*,*)trim(rname)//": startype = ",startype
515
516     if (is_master) write(*,*)trim(rname)//": Value of stellar flux at 1 AU?"
517     Fat1AU=1356.0 ! default value = Sol today
518     call getin_p("Fat1AU",Fat1AU)
519     if (is_master) write(*,*)trim(rname)//": Fat1AU = ",Fat1AU
520
521
522! TRACERS:
523
524     if (is_master) write(*,*)trim(rname)//": fixed radii for Cloud particles?"
525     radfixed=.false. ! default value
526     call getin_p("radfixed",radfixed)
527     if (is_master) write(*,*)trim(rname)//": radfixed = ",radfixed
528
529     if(kastprof)then
530        radfixed=.true.
531     endif
532
533      if (is_master) write(*,*)trim(rname)//&
534       "Number of radiatively active aerosols:"
535     naerkind=0 ! default value
536     call getin_p("naerkind",naerkind)
537     if (is_master) write(*,*)trim(rname)//": naerkind = ",naerkind
538
539
540!***************************************************************
541         !! TRACERS options
542
543     if (is_master)write(*,*)trim(rname)//&
544      "call N2 condensation ?"
545     n2cond=.true. ! default value
546     call getin_p("n2cond",n2cond)
547     if (is_master)write(*,*)trim(rname)//&
548      " n2cond = ",n2cond
549
550     if (is_master)write(*,*)trim(rname)//&
551      "call N2 cloud condensation ?"
552     condensn2=.false. ! default value
553     call getin_p("condensn2",condensn2)
554     if (is_master)write(*,*)trim(rname)//&
555      "condensn2 = ",condensn2
556
557     if (is_master)write(*,*)trim(rname)//&
558      "call no N2 frost formation?"
559     no_n2frost=.false. ! default value
560     call getin_p("no_n2frost",no_n2frost)
561     if (is_master)write(*,*)trim(rname)//&
562      "no_n2frost = ",no_n2frost
563
564     if (is_master)write(*,*)trim(rname)//&
565      "N2 condensation subtimestep?"
566     nbsub=20 ! default value
567     call getin_p("nbsub",nbsub)
568     if (is_master)write(*,*)trim(rname)//&
569      " nbsub = ",nbsub
570
571     if (is_master)write(*,*)trim(rname)//&
572      "Gravitationnal sedimentation ?"
573     sedimentation=.true. ! default value
574     call getin_p("sedimentation",sedimentation)
575     if (is_master)write(*,*)trim(rname)//&
576      " sedimentation = ",sedimentation
577
578     if (is_master)write(*,*)trim(rname)//&
579      "Compute glacial flow ?"
580     glaflow=.false. ! default value
581     call getin_p("glaflow",glaflow)
582     if (is_master)write(*,*)trim(rname)//&
583      " glaflow = ",glaflow
584
585     if (is_master)write(*,*)trim(rname)//&
586      "Compute methane cycle ?"
587     methane=.false. ! default value
588     call getin_p("methane",methane)
589     if (is_master)write(*,*)trim(rname)//&
590      " methane = ",methane
591     condmetsurf=.true. ! default value
592     call getin_p("condmetsurf",condmetsurf)
593     if (is_master)write(*,*)trim(rname)//&
594      " condmetsurf = ",condmetsurf
595
596     if (is_master)write(*,*)trim(rname)//&
597      "Compute methane clouds ?"
598     metcloud=.false. ! default value
599     call getin_p("metcloud",metcloud)
600     if (is_master)write(*,*)trim(rname)//&
601      " metcloud = ",metcloud
602
603     if (is_master)write(*,*)trim(rname)//&
604      "Compute CO cycle ?"
605     carbox=.false. ! default value
606     call getin_p("carbox",carbox)
607     if (is_master)write(*,*)trim(rname)//&
608      " carbox = ",carbox
609     condcosurf=.true. ! default value
610     call getin_p("condcosurf",condcosurf)
611     if (is_master)write(*,*)trim(rname)//&
612      " condcosurf = ",condcosurf
613
614     if (is_master)write(*,*)trim(rname)//&
615      "Compute CO clouds ?"
616     monoxcloud=.false. ! default value
617     call getin_p("monoxcloud",monoxcloud)
618     if (is_master)write(*,*)trim(rname)//&
619      " monoxcloud = ",monoxcloud
620
621     if (is_master)write(*,*)trim(rname)//&
622     "atmospheric redistribution (s):"
623     tau_n2=1. ! default value
624     call getin_p("tau_n2",tau_n2)
625     if (is_master)write(*,*)trim(rname)//&
626     " tau_n2 = ",tau_n2
627     tau_ch4=1.E7 ! default value
628     call getin_p("tau_ch4",tau_ch4)
629     if (is_master)write(*,*)trim(rname)//&
630     " tau_ch4 = ",tau_ch4
631     tau_co=1. ! default value
632     call getin_p("tau_co",tau_co)
633     if (is_master)write(*,*)trim(rname)//&
634     " tau_co = ",tau_co
635     tau_prechaze=1. ! default value
636     call getin_p("tau_prechaze",tau_prechaze)
637     if (is_master)write(*,*)trim(rname)//&
638     " tau_prechaze = ",tau_prechaze
639
640     if (is_master)write(*,*)trim(rname)//&
641      "Day fraction for limited cold trap in SP?"
642     dayfrac=0. ! default value
643     call getin_p("dayfrac",dayfrac)
644     if (is_master)write(*,*)trim(rname)//&
645      " dayfrac = ",dayfrac
646     thresh_non2=0. ! default value
647     call getin_p("thresh_non2",thresh_non2)
648     if (is_master)write(*,*)trim(rname)//&
649      " thresh_non2 = ",thresh_non2
650
651    ! use Pluto.old routines
652     if (is_master) write(*,*) trim(rname)//": use vdifc from old Pluto ?"
653     oldplutovdifc=.false. ! default value
654     call getin_p("oldplutovdifc",oldplutovdifc)
655     if (is_master) write(*,*) trim(rname)//": oldplutovdifc = ",oldplutovdifc
656
657     if (is_master) write(*,*) trim(rname)//&
658       ": call pluto.old correlated-k radiative transfer ?"
659     oldplutocorrk=.false. ! default value
660     call getin_p("oldplutocorrk",oldplutocorrk)
661     if (is_master) write(*,*) trim(rname)//": oldplutocorrk = ",oldplutocorrk
662
663     if (is_master) write(*,*) trim(rname)//&
664       ": call pluto.old sedimentation ?"
665     oldplutosedim=.false. ! default value
666     call getin_p("oldplutosedim",oldplutosedim)
667     if (is_master) write(*,*) trim(rname)//": oldplutosedim = ",oldplutosedim
668
669!***************************************************************
670     !! Haze options
671
672     ! Microphysical moment model
673     ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
674     if (is_master) write(*,*) "Run with or without microphysics?"
675     callmufi=.false. ! default value
676     call getin_p("callmufi",callmufi)
677     if (is_master) write(*,*)" callmufi = ",callmufi
678     
679     ! sanity check
680     if (callmufi.and.(.not.tracer)) then
681       print*,"You are running microphysics without tracer"
682       print*,"Please start again with tracer =.true."
683       stop
684     endif
685 
686     if (is_master) write(*,*) "Path to microphysical config file?"
687     config_mufi='datagcm/microphysics/config.cfg' ! default value
688     call getin_p("config_mufi",config_mufi)
689     if (is_master) write(*,*)" config_mufi = ",config_mufi
690
691     if (is_master) write(*,*) "Spherical aerosol optical properties datafile"
692     aersprop_file="optprop_rannou_r2-200nm_nu003.dat"  ! default file
693     call getin_p("aersprop_file",aersprop_file)
694     if (is_master) write(*,*) trim(rname)//" aersprop_file = ",trim(aersprop_file)
695
696     if (is_master) write(*,*) "Fractal aerosol optical properties datafile"
697     aerfprop_file="optprop_rannou_fractal_r010nm_N1_1e4_d2.dat"  ! default file
698     call getin_p("aerfprop_file",aerfprop_file)
699     if (is_master) write(*,*) trim(rname)//" aerfprop_file = ",trim(aerfprop_file)
700 
701     if (is_master) write(*,*) "Use haze production from CH4 photolysis or production rate?"
702     call_haze_prod_pCH4=.false. ! default value
703     call getin_p("call_haze_prod_pCH4",call_haze_prod_pCH4)
704     if (is_master) write(*,*)" call_haze_prod_pCH4 = ",call_haze_prod_pCH4
705 
706     if (is_master) write(*,*) "Pressure level of aerosols production (Pa)?"
707     haze_p_prod=1.0e-2 ! default value
708     call getin_p("haze_p_prod",haze_p_prod)
709     if (is_master) write(*,*)" haze_p_prod = ",haze_p_prod
710     
711     if (is_master) write(*,*) "Aerosol production rate (kg.m-2.s-1)?"
712     haze_tx_prod=9.8e-14 ! default value
713     call getin_p("haze_tx_prod",haze_tx_prod)
714     if (is_master) write(*,*)" haze_tx_prod = ",haze_tx_prod
715 
716     if (is_master) write(*,*) "Equivalent radius production (m)?"
717     haze_rc_prod=1.0e-9 ! default value
718     call getin_p("haze_rc_prod",haze_rc_prod)
719     if (is_master) write(*,*)" haze_rc_prod = ",haze_rc_prod
720 
721     if (is_master) write(*,*) "Monomer radius (m)?"
722     haze_rm=1.0e-8 ! default value
723     call getin_p("haze_rm",haze_rm)
724     if (is_master) write(*,*)" haze_rm = ",haze_rm
725 
726     if (is_master) write(*,*) "Aerosol's fractal dimension?"
727     haze_df=2.0 ! default value
728     call getin_p("haze_df",haze_df)
729     if (is_master) write(*,*)" haze_df = ",haze_df
730 
731     if (is_master) write(*,*) "Aerosol density (kg.m-3)?"
732     haze_rho=800.0 ! default value
733     call getin_p("haze_rho",haze_rho)
734     if (is_master) write(*,*)" haze_rho = ",haze_rho
735 
736     if (is_master) write(*,*) "Radius of air molecule (m)?"
737     air_rad=1.75e-10 ! default value
738     call getin_p("air_rad",air_rad)
739     if (is_master) write(*,*)" air_rad = ",air_rad
740     
741     ! Pluto haze model
742     ! ~~~~~~~~~~~~~~~~
743     if (is_master)write(*,*)trim(rname)//&
744     "Production of haze ?"
745     haze=.false. ! default value
746     call getin_p("haze",haze)
747     if (is_master)write(*,*)trim(rname)//&
748     " haze = ",haze
749     hazeconservch4=.false. ! conservch4, by default value ch4 is photolyzed
750     call getin_p("hazeconservch4",hazeconservch4)
751     if (is_master)write(*,*)trim(rname)//&
752     "hazconservch4 = ",hazeconservch4
753     if (is_master)write(*,*)trim(rname)//&
754     "Production of haze (fast version) ?"
755     fasthaze=.false. ! default value
756     call getin_p("fasthaze",fasthaze)
757     if (is_master)write(*,*)trim(rname)//&
758     "fasthaze = ",fasthaze
759
760     if (is_master)write(*,*)trim(rname)//&
761     "Add source of haze ?"
762     source_haze=.false. ! default value
763     call getin_p("source_haze",source_haze)
764     if (is_master)write(*,*)trim(rname)//&
765     " source_haze = ",source_haze
766     mode_hs=0 ! mode haze source default value
767     call getin_p("mode_hs",mode_hs)
768     if (is_master)write(*,*)trim(rname)//&
769     " mode_hs",mode_hs
770     kfix=1 ! default value
771     call getin_p("kfix",kfix)
772     if (is_master)write(*,*)trim(rname)//&
773     "levels kfix",kfix
774     fracsource=0.1 ! default value
775     call getin_p("fracsource",fracsource)
776     if (is_master)write(*,*)trim(rname)//&
777     " fracsource",fracsource
778     latsource=30. ! default value
779     call getin_p("latsource",latsource)
780     if (is_master)write(*,*)trim(rname)//&
781     " latsource",latsource
782     lonsource=180. ! default value
783     call getin_p("lonsource",lonsource)
784     if (is_master)write(*,*)trim(rname)//&
785     " lonsource",lonsource
786
787     if (is_master)write(*,*)trim(rname)//&
788     "Radiatively active haze ?"
789     optichaze=.false. ! default value
790     call getin_p("optichaze",optichaze)
791     if (is_master)write(*,*)trim(rname)//&
792     "optichaze = ",optichaze
793
794     if (is_master)write(*,*)trim(rname)//&
795     "Haze monomer radius ?"
796     rad_haze=10.e-9 ! default value
797     call getin_p("rad_haze",rad_haze)
798     if (is_master)write(*,*)trim(rname)//&
799     "rad_haze = ",rad_haze
800
801     if (is_master)write(*,*)trim(rname)//&
802     "fractal particle ?"
803     fractal=.false. ! default value
804     call getin_p("fractal",fractal)
805     if (is_master)write(*,*)trim(rname)//&
806     "fractal = ",fractal
807     nb_monomer=10 ! default value
808     call getin_p("nb_monomer",nb_monomer)
809     if (is_master)write(*,*)trim(rname)//&
810     "nb_monomer = ",nb_monomer
811
812     if (is_master)write(*,*)trim(rname)//&
813     "fixed gaseous CH4 mixing ratio for rad transfer ?"
814    ch4fix=.false. ! default value
815    call getin_p("ch4fix",ch4fix)
816    if (is_master)write(*,*)trim(rname)//&
817     " ch4fix = ",ch4fix
818    if (is_master)write(*,*)trim(rname)//&
819     "fixed gaseous CH4 mixing ratio for rad transfer ?"
820    vmrch4fix=0.5 ! default value
821    call getin_p("vmrch4fix",vmrch4fix)
822    if (is_master)write(*,*)trim(rname)//&
823     " vmrch4fix = ",vmrch4fix
824    vmrch4_proffix=.false. ! default value
825    call getin_p("vmrch4_proffix",vmrch4_proffix)
826    if (is_master)write(*,*)trim(rname)//&
827     " vmrch4_proffix = ",vmrch4_proffix
828
829    if (is_master)write(*,*)trim(rname)//&
830     "call specific cooling for radiative transfer ?"
831    cooling=.false.  ! default value
832    call getin_p("cooling",cooling)
833    if (is_master)write(*,*)trim(rname)//&
834     " cooling = ",cooling
835
836    if (is_master)write(*,*)trim(rname)//&
837     "NLTE correction for methane heating rates?"
838    nlte=.false.  ! default value
839    call getin_p("nlte",nlte)
840    if (is_master)write(*,*)trim(rname)//&
841     " nlte = ",nlte
842    strobel=.false.  ! default value
843    call getin_p("strobel",strobel)
844    if (is_master)write(*,*)trim(rname)//&
845     " strobel = ",strobel
846
847     if (is_master)write(*,*)trim(rname)//&
848     "fixed radius profile from txt file ?"
849     haze_radproffix=.false. ! default value
850     call getin_p("haze_radproffix",haze_radproffix)
851     if (is_master)write(*,*)trim(rname)//&
852     "haze_radproffix = ",haze_radproffix
853     if (is_master)write(*,*)trim(rname)//&
854     "fixed MMR profile from txt file ?"
855     haze_proffix=.false. ! default value
856     call getin_p("haze_proffix",haze_proffix)
857     if (is_master)write(*,*)trim(rname)//&
858     "haze_proffix = ",haze_proffix
859
860     if (is_master)write(*,*)trim(rname)//&
861     "Number mix ratio of haze particles for co clouds:"
862     Nmix_co=100000. ! default value
863     call getin_p("Nmix_co",Nmix_co)
864     if (is_master)write(*,*)trim(rname)//&
865     " Nmix_co = ",Nmix_co
866
867     if (is_master)write(*,*)trim(rname)//&
868     "Number mix ratio of haze particles for ch4 clouds:"
869     Nmix_ch4=100000. ! default value
870     call getin_p("Nmix_ch4",Nmix_ch4)
871     if (is_master)write(*,*)trim(rname)//&
872     " Nmix_ch4 = ",Nmix_ch4
873
874!***************************************************************
875     !! Surface properties
876
877!*********** N2 *********************************
878
879     if (is_master)write(*,*)trim(rname)//&
880      "Mode for changing N2 albedo"
881     mode_n2=0 ! default value
882     call getin_p("mode_n2",mode_n2)
883     if (is_master)write(*,*)trim(rname)//&
884      " mode_n2 = ",mode_n2
885     thres_n2ice=1. ! default value
886     call getin_p("thres_n2ice",thres_n2ice)
887     if (is_master)write(*,*)trim(rname)//&
888      " thres_n2ice = ",thres_n2ice
889
890     if (is_master)write(*,*)trim(rname)//&
891      "Diff of N2 albedo with thickness"
892     deltab=0. ! default value
893     call getin_p("deltab",deltab)
894     if (is_master)write(*,*)trim(rname)//&
895      " deltab = ",deltab
896
897     if (is_master)write(*,*)trim(rname)//&
898      "albedo N2 beta "
899     alb_n2b=0.7 ! default value
900     call getin_p("alb_n2b",alb_n2b)
901     if (is_master)write(*,*)trim(rname)//&
902      " alb_n2b = ",alb_n2b
903
904     if (is_master)write(*,*)trim(rname)//&
905      "albedo N2 alpha "
906     alb_n2a=0.7 ! default value
907     call getin_p("alb_n2a",alb_n2a)
908     if (is_master)write(*,*)trim(rname)//&
909      " alb_n2a = ",alb_n2a
910
911     if (is_master)write(*,*)trim(rname)//&
912      "emis N2 beta "
913     emis_n2b=0.7 ! default value
914     call getin_p("emis_n2b",emis_n2b)
915     if (is_master)write(*,*)trim(rname)//&
916      " emis_n2b = ",emis_n2b
917
918     if (is_master)write(*,*)trim(rname)//&
919      "emis N2 alpha "
920     emis_n2a=0.7 ! default value
921     call getin_p("emis_n2a",emis_n2a)
922     if (is_master)write(*,*)trim(rname)//&
923      " emis_n2a = ",emis_n2a
924
925!*********** CH4 *********************************
926
927     if (is_master)write(*,*)trim(rname)//&
928      "Mode for changing CH4 albedo"
929     mode_ch4=0 ! default value
930     call getin_p("mode_ch4",mode_ch4)
931     if (is_master)write(*,*)trim(rname)//&
932      " mode_ch4 = ",mode_ch4
933     feedback_met=0 ! default value
934     call getin_p("feedback_met",feedback_met)
935     if (is_master)write(*,*)trim(rname)//&
936      " feedback_met = ",feedback_met
937     thres_ch4ice=1. ! default value, mm
938     call getin_p("thres_ch4ice",thres_ch4ice)
939     if (is_master)write(*,*)trim(rname)//&
940      " thres_ch4ice = ",thres_ch4ice
941     fdch4_latn=-1.
942     fdch4_lats=0.
943     fdch4_lone=0.
944     fdch4_lonw=-1.
945     fdch4_depalb=0.5
946     fdch4_finalb=0.95
947     fdch4_maxalb=0.99
948     fdch4_ampl=200.
949     fdch4_maxice=100.
950     call getin_p("fdch4_latn",fdch4_latn)
951     call getin_p("fdch4_lats",fdch4_lats)
952     call getin_p("fdch4_lone",fdch4_lone)
953     call getin_p("fdch4_lonw",fdch4_lonw)
954     call getin_p("fdch4_depalb",fdch4_depalb)
955     call getin_p("fdch4_finalb",fdch4_finalb)
956     call getin_p("fdch4_maxalb",fdch4_maxalb)
957     call getin_p("fdch4_maxice",fdch4_maxice)
958     call getin_p("fdch4_ampl",fdch4_ampl)
959     if (is_master)write(*,*)trim(rname)//&
960      " Values for albedo feedback = ",fdch4_latn,&
961     fdch4_lats,fdch4_lone,fdch4_lonw,fdch4_depalb,&
962     fdch4_finalb,fdch4_maxalb,fdch4_maxice,fdch4_ampl
963
964     if (is_master)write(*,*)trim(rname)//&
965      "Latitude for diff albedo methane"
966     metlateq=25. ! default value
967     call getin_p("metlateq",metlateq)
968     if (is_master)write(*,*)trim(rname)//&
969      " metlateq = ",metlateq
970
971     if (is_master)write(*,*)trim(rname)//&
972      "Ls1 and Ls2 for change of ch4 albedo"
973     metls1=-1. ! default value
974     metls2=-2. ! default value
975     call getin_p("metls1",metls1)
976     call getin_p("metls2",metls2)
977
978     if (is_master)write(*,*)trim(rname)//&
979      "albedo CH4 "
980     alb_ch4=0.5 ! default value
981     call getin_p("alb_ch4",alb_ch4)
982     if (is_master)write(*,*)trim(rname)//&
983      " alb_ch4 = ",alb_ch4
984
985     if (is_master)write(*,*)trim(rname)//&
986      "albedo equatorial CH4 "
987     alb_ch4_eq=alb_ch4 ! default value
988     call getin_p("alb_ch4_eq",alb_ch4_eq)
989     if (is_master)write(*,*)trim(rname)//&
990      " alb_ch4_eq = ",alb_ch4_eq
991
992     if (is_master)write(*,*)trim(rname)//&
993      "albedo south hemis CH4 "
994     alb_ch4_s=alb_ch4 ! default value
995     call getin_p("alb_ch4_s",alb_ch4_s)
996     if (is_master)write(*,*)trim(rname)//&
997      " alb_ch4_s = ",alb_ch4_s
998
999     if (is_master)write(*,*)trim(rname)//&
1000      "emis CH4 "
1001     emis_ch4=0.5 ! default value
1002     call getin_p("emis_ch4",emis_ch4)
1003     if (is_master)write(*,*)trim(rname)//&
1004      " emis_ch4 = ",emis_ch4
1005
1006     if (is_master)write(*,*)trim(rname)//&
1007      "CH4 lag for n2 sublimation limitation"
1008     ch4lag=.false. ! default value
1009     latlag=-90. ! default value
1010     vmrlag=1. ! default value
1011     call getin_p("ch4lag",ch4lag)
1012     call getin_p("latlag",latlag)
1013     call getin_p("vmrlag",vmrlag)
1014     if (is_master)write(*,*)trim(rname)//&
1015      " ch4lag = ",ch4lag
1016     if (is_master)write(*,*)trim(rname)//&
1017      " latlag = ",latlag
1018     if (is_master)write(*,*)trim(rname)//&
1019      " vmrlag = ",vmrlag
1020
1021     if (is_master)write(*,*)trim(rname)//&
1022      "Max temperature for surface ?"
1023     tsurfmax=.false. ! default value
1024     albmin_ch4=0.3 ! default value
1025     call getin_p("tsurfmax",tsurfmax)
1026     call getin_p("albmin_ch4",albmin_ch4)
1027     if (is_master)write(*,*)trim(rname)//&
1028      " tsurfmax = ",tsurfmax
1029     if (is_master)write(*,*)trim(rname)//&
1030      " albmin_ch4 = ",albmin_ch4
1031
1032
1033     if (is_master)write(*,*)trim(rname)//&
1034     "fixed gaseous CH4 mixing ratio for rad transfer ?"
1035     ch4fix=.false. ! default value
1036     call getin_p("ch4fix",ch4fix)
1037     if (is_master)write(*,*)trim(rname)//&
1038       " ch4fix = ",ch4fix
1039     if (is_master)write(*,*)trim(rname)//&
1040       "fixed gaseous CH4 mixing ratio for rad transfer ?"
1041     vmrch4fix=0.5 ! default value
1042     call getin_p("vmrch4fix",vmrch4fix)
1043     if (is_master)write(*,*)trim(rname)//&
1044        " vmrch4fix = ",vmrch4fix
1045     vmrch4_proffix=.false. ! default value
1046     call getin_p("vmrch4_proffix",vmrch4_proffix)
1047     if (is_master)write(*,*)trim(rname)//&
1048        " vmrch4_proffix = ",vmrch4_proffix
1049
1050
1051!*********** CO *********************************
1052
1053     if (is_master)write(*,*)trim(rname)//&
1054      "albedo CO "
1055     alb_co=0.4 ! default value
1056     call getin_p("alb_co",alb_co)
1057     if (is_master)write(*,*)trim(rname)//&
1058      " alb_co = ",alb_co
1059     thres_coice=1. ! default value, mm
1060     call getin_p("thres_coice",thres_coice)
1061     if (is_master)write(*,*)trim(rname)//&
1062      " thres_coice = ",thres_coice
1063
1064     if (is_master)write(*,*)trim(rname)//&
1065      "emis CO "
1066     emis_co=0.4 ! default value
1067     call getin_p("emis_co",emis_co)
1068     if (is_master)write(*,*)trim(rname)//&
1069      " emis_co = ",emis_co
1070
1071!*********** THOLINS *********************************
1072     if (is_master)write(*,*)trim(rname)//&
1073      "Mode for changing tholins albedo/emis"
1074     mode_tholins=0 ! default value
1075     call getin_p("mode_tholins",mode_tholins)
1076     if (is_master)write(*,*)trim(rname)//&
1077      " mode_tholins = ",mode_tholins
1078
1079     if (is_master)write(*,*)trim(rname)//&
1080      "albedo tho "
1081     alb_tho=0.1 ! default value
1082     call getin_p("alb_tho",alb_tho)
1083     if (is_master)write(*,*)trim(rname)//&
1084      " alb_tho = ",alb_tho
1085
1086     if (is_master)write(*,*)trim(rname)//&
1087      "albedo tho eq"
1088     alb_tho_eq=0.1 ! default value
1089     call getin_p("alb_tho_eq",alb_tho_eq)
1090     if (is_master)write(*,*)trim(rname)//&
1091      " alb_tho_eq = ",alb_tho_eq
1092
1093     if (is_master)write(*,*)trim(rname)//&
1094      "emis tholins "
1095     emis_tho=1. ! default value
1096     call getin_p("emis_tho",emis_tho)
1097     if (is_master)write(*,*)trim(rname)//&
1098      " emis_tho = ",emis_tho
1099
1100     if (is_master)write(*,*)trim(rname)//&
1101      "emis tholins eq"
1102     emis_tho_eq=1. ! default value
1103     call getin_p("emis_tho_eq",emis_tho_eq)
1104     if (is_master)write(*,*)trim(rname)//&
1105      " emis_tho_eq = ",emis_tho_eq
1106
1107     if (is_master)write(*,*)trim(rname)//&
1108      "Latitude for diff albedo tholins"
1109     tholateq=25. ! default value
1110     call getin_p("tholateq",tholateq)
1111     if (is_master)write(*,*)trim(rname)//&
1112      " tholateq = ",tholateq
1113     tholatn=-1.
1114     tholats=0.
1115     tholone=0.
1116     tholonw=-1.
1117     alb_tho_spe=0.1 ! default value
1118     emis_tho_spe=1. ! default value
1119     call getin_p("  tholatn",tholatn)
1120     call getin_p("  tholats",tholats)
1121     call getin_p("  tholone",tholone)
1122     call getin_p("  tholonw",tholonw)
1123     if (is_master)write(*,*)trim(rname)//&
1124      " Values for special tholins albedo = ",tholatn,&
1125        tholats,tholone,tholonw,alb_tho_spe,emis_tho_spe
1126
1127     if (is_master)write(*,*)trim(rname)//&
1128      "Specific albedo"
1129     spelon1=-180. ! default value
1130     spelon2=180. ! default value
1131     spelat1=-90. ! default value
1132     spelat2=90. ! default value
1133     specalb=.false. ! default value
1134     if (is_master)write(*,*)trim(rname)//&
1135      "albedo/emis spe "
1136     albspe=0.1 ! default value
1137     emispe=1. ! default value
1138     call getin_p("spelon1",spelon1)
1139     call getin_p("spelon2",spelon2)
1140     call getin_p("spelat1",spelat1)
1141     call getin_p("spelat2",spelat2)
1142     call getin_p("specalb",specalb)
1143     call getin_p("albspe",albspe)
1144     call getin_p("emispe",emispe)
1145
1146     if (is_master)write(*,*)trim(rname)//&
1147      " specific = ",specalb
1148
1149!********************** TI *****************************
1150
1151     if (is_master)write(*,*)trim(rname)//&
1152      "Change TI with time"
1153     changeti=.false. ! default value
1154     call getin_p("changeti",changeti)
1155     if (is_master)write(*,*)trim(rname)//&
1156      " changeti = ",changeti
1157     changetid=.false. ! default value for diurnal TI
1158     call getin_p("changetid",changetid)
1159     if (is_master)write(*,*)trim(rname)//&
1160      " changetid = ",changetid
1161
1162     if (is_master)write(*,*)trim(rname)//&
1163      "IT N2 "
1164     ITN2=800. ! default value
1165     call getin_p("ITN2",ITN2)
1166     if (is_master)write(*,*)trim(rname)//&
1167      " ITN2 = ",ITN2
1168     ITN2d=20. ! default value
1169     call getin_p("ITN2d",ITN2d)
1170     if (is_master)write(*,*)trim(rname)//&
1171      " ITN2d = ",ITN2d
1172
1173     if (is_master)write(*,*)trim(rname)//&
1174      "IT CH4"
1175     ITCH4=800. ! default value
1176     call getin_p("ITCH4",ITCH4)
1177     if (is_master)write(*,*)trim(rname)//&
1178      " ITCH4 = ",ITCH4
1179     ITCH4d=20. ! default value
1180     call getin_p("ITCH4d",ITCH4d)
1181     if (is_master)write(*,*)trim(rname)//&
1182      " ITCH4d = ",ITCH4d
1183
1184     if (is_master)write(*,*)trim(rname)//&
1185      "IT H2O"
1186     ITH2O=800. ! default value
1187     call getin_p("ITH2O",ITH2O)
1188     if (is_master)write(*,*)trim(rname)//&
1189      " ITH2O = ",ITH2O
1190     ITH2Od=20. ! default value
1191     call getin_p("ITH2Od",ITH2Od)
1192     if (is_master)write(*,*)trim(rname)//&
1193      " ITH2Od = ",ITH2Od
1194
1195!************** COOLING ***************
1196
1197     alpha_top=5e-11 ! default value
1198     call getin_p("alpha_top",alpha_top)
1199     if (is_master)write(*,*)trim(rname)//&
1200      " alpha_top = ",alpha_top
1201     pref=0.02 ! default value
1202     call getin_p("pref",pref)
1203     if (is_master)write(*,*)trim(rname)//&
1204      " pref = ",pref
1205     deltap=0.1 ! default value
1206     call getin_p("deltap",deltap)
1207     if (is_master)write(*,*)trim(rname)//&
1208      " deltap = ",deltap
1209
1210!=================================
1211
1212     if (is_master) write(*,*)trim(rname)//&
1213       ": Is the variable gas species radiatively active?"
1214     Tstrat=167.0
1215     varactive=.false.
1216     call getin_p("varactive",varactive)
1217     if (is_master) write(*,*)trim(rname)//": varactive = ",varactive
1218
1219     if (is_master) write(*,*)trim(rname)//&
1220       ": Is the variable gas species distribution set?"
1221     varfixed=.false.
1222     call getin_p("varfixed",varfixed)
1223     if (is_master) write(*,*)trim(rname)//": varfixed = ",varfixed
1224
1225     if (is_master) write(*,*)trim(rname)//&
1226       ": What is the saturation % of the variable species?"
1227     satval=0.8
1228     call getin_p("satval",satval)
1229     if (is_master) write(*,*)trim(rname)//": satval = ",satval
1230
1231
1232! Test of incompatibility:
1233! if varactive, then varfixed should be false
1234     if (varactive.and.varfixed) then
1235       call abort_physic(rname,'if varactive, varfixed must be OFF!',1)
1236     endif
1237
1238     if (is_master) write(*,*)trim(rname)//": Gravitationnal sedimentation ?"
1239     sedimentation=.false. ! default value
1240     call getin_p("sedimentation",sedimentation)
1241     if (is_master) write(*,*)trim(rname)//": sedimentation = ",sedimentation
1242
1243     if (is_master) write(*,*)trim(rname)//": Spectral Dependant albedo ?"
1244     albedo_spectral_mode=.false. ! default value
1245     call getin_p("albedo_spectral_mode",albedo_spectral_mode)
1246     if (is_master) write(*,*)trim(rname)//&
1247     ": albedo_spectral_mode = ",albedo_spectral_mode
1248
1249     if (is_master) then
1250       write(*,*)trim(rname)//": Snow albedo ?"
1251       write(*,*)trim(rname)//": If albedo_spectral_mode=.true., then this "
1252       write(*,*)trim(rname)//": corresponds to the 0.5 microns snow albedo."
1253     endif
1254     albedosnow=0.5         ! default value
1255     call getin_p("albedosnow",albedosnow)
1256     if (is_master) write(*,*)trim(rname)//": albedosnow = ",albedosnow
1257
1258     if (is_master) write(*,*)trim(rname)//": N2 ice albedo ?"
1259     albedon2ice=0.5       ! default value
1260     call getin_p("albedon2ice",albedon2ice)
1261     if (is_master) write(*,*)trim(rname)//": albedon2ice = ",albedon2ice
1262
1263     if (is_master) write(*,*)trim(rname)//": Maximum ice thickness ?"
1264     maxicethick=2.0         ! default value
1265     call getin_p("maxicethick",maxicethick)
1266     if (is_master) write(*,*)trim(rname)//": maxicethick = ",maxicethick
1267
1268     if (is_master) write(*,*)trim(rname)//": Minimum eddy mix coeff in 1D ?"
1269     kmixmin=1.0e-2         ! default value
1270     call getin_p("kmixmin",kmixmin)
1271     if (is_master) write(*,*)trim(rname)//": kmixmin = ",kmixmin
1272
1273     if (is_master) write(*,*)'Predefined Cp from dynamics is ',cpp,'J kg^-1 K^-1'
1274     if (is_master) write(*,*)'Predefined Mg from dynamics is ',mugaz,'amu'
1275
1276     force_cpp=.false. ! default value
1277     call getin_p("force_cpp",force_cpp)
1278     if (force_cpp) then
1279      if (is_master) write(*,*)trim(rname)//": force_cpp = ",force_cpp
1280      if (is_master) write(*,*)trim(rname)//": force_cpp is deprecated.",&
1281      "Set cpp_mugaz_mode=1 in callfis to emulate force_cpp=.true."
1282      call abort_physic(rname,"Anyway, you need to set force_cpp=.false. to continue.",1)
1283     endif
1284
1285     if (is_master) write(*,*)trim(rname)//&
1286     ": where do you want your cpp/mugaz value to come from?",&
1287     "=> 0: dynamics (3d), 1: forced in callfis (1d), 2: computed from gases.def (1d)?"
1288     cpp_mugaz_mode = 0 ! default value
1289     call getin_p("cpp_mugaz_mode",cpp_mugaz_mode)
1290     if (is_master) write(*,*)trim(rname)//": cpp_mugaz_mode = ",cpp_mugaz_mode
1291
1292
1293! Test of incompatibility:
1294
1295     if ((.not.tracer).and.(haze)) then
1296       call abort_physic(rname, 'if haze are on, tracers must be on!', 1)
1297     endif
1298     if ((callmufi).and.(haze)) then
1299       call abort_physic(rname, 'if haze are on, microphysics should be desactivated!', 1)
1300     endif
1301     if ((haze).and.(naerkind.gt.1)) then
1302      call abort_physic(rname, 'if haze are on, naerkind must be equal to 1!', 1)
1303     endif
1304     if ((callmufi).and..not.(naerkind.gt.1)) then
1305      call abort_physic(rname, 'if microphysics is on, naerkind must be > 1!', 1)
1306     endif
1307     if (.not.(callmufi.or.haze).and.(optichaze)) then
1308      call abort_physic(rname, 'if microphysics and haze are off, optichaze must be desactivated!', 1)
1309     endif
1310     if ((callmufi.and.call_haze_prod_pCH4).and..not.(methane)) then
1311      call abort_physic(rname, 'if haze production from CH4 photolysis is on, methane must be activated!', 1)
1312     endif
1313     if (haze_proffix.and.sedimentation) then
1314         call abort_physic(rname, 'if haze profile is set, sedimentation must be deactivated', 1)
1315     endif
1316     if (callmolvis.and..not.callconduct) then
1317         call abort_physic(rname, 'if callmolvis is set, callconduct must be true', 1)
1318     endif
1319     if (glaflow.and..not.fast) then
1320         call abort_physic(rname, 'if glaflow is set, fast must be true', 1)
1321     endif
1322     if (paleo.and..not.fast) then
1323         call abort_physic(rname, 'if paleo is set, fast must be true', 1)
1324     endif
1325     if ((haze_proffix.or.haze_radproffix).and..not.optichaze) then
1326         call abort_physic(rname, 'for now, haze/rad_proffix only works with optichaze=T', 1)
1327     endif
1328      if (carbox.and.condcosurf.and.no_n2frost) then
1329        call abort_physic(rname, "CO surface condensation and no_n2frost are both active which may not be relevant", 1)
1330      end if
1331
1332     if ((cpp_mugaz_mode >= 1).and.(is_master).and.(ngrid>1)) then
1333        write(*,*)'    !!!!! Be aware that having different values for cpp and mugaz in the dynamics and physics'
1334        write(*,*)'    in 3D can result in very pathological behavior. You have been warned !!!!!'
1335      else if ((cpp_mugaz_mode == 0).and.(is_master).and.(ngrid==1)) then
1336        ! for this specific 1D error we will remove run.def before aborting JL22
1337        call system("rm -rf run.def")
1338        call abort_physic(rname,"cpp_mugaz_mode must be >= 1 in 1d",1)
1339     endif
1340
1341     if (cpp_mugaz_mode == 1) then
1342       mugaz = -99999.
1343       if (is_master) write(*,*)trim(rname)//&
1344         ": MEAN MOLECULAR MASS in g mol-1 ?"
1345       call getin_p("mugaz",mugaz)
1346       IF (mugaz.eq.-99999.) THEN
1347         call abort_physic(rname,"mugaz must be set if cpp_mugaz_mode = 1",1)
1348       ENDIF
1349       cpp = -99999.
1350       if (is_master) write(*,*)trim(rname)//&
1351         ": SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
1352       call getin_p("cpp",cpp)
1353       IF (cpp.eq.-99999.) THEN
1354           PRINT *, "cpp must be set if cpp_mugaz_mode = 1"
1355           STOP
1356       ENDIF
1357       if (is_master) write(*,*)'New Cp from callfis is ',cpp,'J kg^-1 K^-1'
1358       if (is_master) write(*,*)'New Mg from callfis is ',mugaz,'amu'
1359
1360     endif ! of if (cpp_mugaz_mode == 1)
1361     call su_gases
1362     call calc_cpp_mugaz
1363
1364     if (is_master) then
1365       PRINT*,'--------------------------------------------'
1366       PRINT*
1367       PRINT*
1368     endif
1369  ELSE
1370     call abort_physic(rname,'Cannot read file callphys.def. Is it here ?',1)
1371  ENDIF ! of IF(iscallphys)
1372
1373  if (is_master) then
1374    PRINT*
1375    PRINT*,'inifis: daysec',daysec
1376    PRINT*
1377    PRINT*,'inifis: The radiative transfer is computed:'
1378    PRINT*,'           each ',iradia,' physical time-step'
1379    PRINT*,'        or each ',iradia*dtphys,' seconds'
1380    PRINT*
1381  endif
1382
1383!-----------------------------------------------------------------------
1384!     Some more initialization:
1385!     ------------------------
1386
1387  ! Initializations for comgeomfi_h
1388#ifndef MESOSCALE
1389  totarea=SSUM(ngrid,parea,1)
1390  call planetwide_sumval(parea,totarea_planet)
1391
1392  !! those are defined in comdiurn_h.F90
1393  IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
1394  IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
1395  IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
1396  IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
1397
1398  DO ig=1,ngrid
1399     sinlat(ig)=sin(plat(ig))
1400     coslat(ig)=cos(plat(ig))
1401     sinlon(ig)=sin(plon(ig))
1402     coslon(ig)=cos(plon(ig))
1403  ENDDO
1404#endif
1405  ! initialize variables in radinc_h
1406  call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
1407
1408  ! initialize variables and allocate arrays in radcommon_h
1409  call ini_radcommon_h(naerkind)
1410
1411  ! allocate "comsoil_h" arrays
1412  call ini_comsoil_h(ngrid)
1413
1414  END SUBROUTINE inifis
1415
1416END MODULE inifis_mod
Note: See TracBrowser for help on using the repository browser.