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

Last change on this file since 3477 was 3477, checked in by afalco, 4 weeks ago

Pluto PCM: fixed some variables imports for 'libphy' compilation (to compile with DYNAMICO).
AF

File size: 57.4 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, Nmix_n2
15  use datafile_mod, only: datadir,hazeprop_file,hazerad_file,hazemmr_file,hazedens_file
16  use comdiurn_h, only: sinlat, coslat, sinlon, coslon
17  use comgeomfi_h, only: totarea, totarea_planet
18  use comsoil_h, only: ini_comsoil_h, nsoilmx, lay1_soil, alpha_soil
19  use time_phylmdz_mod, only: ecritphy,day_step,iphysiq, &
20                              init_time, daysec, dtphys
21  use comcstfi_mod, only: rad, cpp, g, r, rcp, &
22                          mugaz, pi, avocado
23  use planete_mod, only: nres
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     conservn2=.false. ! default value
206     call getin_p("conservn2",conservn2)
207     if (is_master) write(*,*) trim(rname)//" conservn2 = ",conservn2
208
209     if (is_master) write(*,*) trim(rname)//"KBO runs (eris, makemake) ?"
210     kbo=.false. ! default value
211     call getin_p("kbo",kbo)
212     if (is_master) write(*,*) trim(rname)//" kbo = ",kbo
213
214     if (is_master) write(*,*) trim(rname)//"Specific paleo run ?"
215     paleo=.false. ! default value
216     call getin_p("paleo",paleo)
217     if (is_master) write(*,*) trim(rname)//" paleo = ",paleo
218
219     if (is_master) write(*,*) trim(rname)//"paleoclimate step (Earth years) "
220     paleoyears=10000. ! default value
221     call getin_p("paleoyears",paleoyears)
222     if (is_master) write(*,*) trim(rname)//" paleoyears = ",paleoyears
223
224     if (is_master) write(*,*) trim(rname)//&
225       ": No seasonal cycle: initial day to lock the run during restart"
226     noseason_day=0.0 ! default value
227     call getin_p("noseason_day",noseason_day)
228     if (is_master) write(*,*) trim(rname)//": noseason_day=", noseason_day
229
230     if (is_master) write(*,*) trim(rname)//": Tidal resonance ratio ?"
231     nres=0          ! default value
232     call getin_p("nres",nres)
233     if (is_master) write(*,*) trim(rname)//": nres = ",nres
234
235     if (is_master) write(*,*) trim(rname)//&
236       ": Write some extra output to the screen ?"
237     lwrite=.false. ! default value
238     call getin_p("lwrite",lwrite)
239     if (is_master) write(*,*) trim(rname)//": lwrite = ",lwrite
240
241     if (is_master) write(*,*) trim(rname)//&
242       ": Save statistics in file stats.nc ?"
243     callstats=.true. ! default value
244     call getin_p("callstats",callstats)
245     if (is_master) write(*,*) trim(rname)//": callstats = ",callstats
246
247     if (is_master) write(*,*) trim(rname)//&
248       ": Test energy conservation of model physics ?"
249     enertest=.false. ! default value
250     call getin_p("enertest",enertest)
251     if (is_master) write(*,*) trim(rname)//": enertest = ",enertest
252
253     if (is_master) write(*,*) trim(rname)//": call radiative transfer ?"
254     callrad=.true. ! default value
255     call getin_p("callrad",callrad)
256     if (is_master) write(*,*) trim(rname)//": callrad = ",callrad
257
258     if (is_master) write(*,*) trim(rname)//&
259       ": call correlated-k radiative transfer ?"
260     corrk=.true. ! default value
261     call getin_p("corrk",corrk)
262     if (is_master) write(*,*) trim(rname)//": corrk = ",corrk
263
264     if (is_master) write(*,*) trim(rname)//&
265       ": prohibit calculations outside corrk T grid?"
266     strictboundcorrk=.true. ! default value
267     call getin_p("strictboundcorrk",strictboundcorrk)
268     if (is_master) write(*,*) trim(rname)//&
269       ": strictboundcorrk = ",strictboundcorrk
270
271     if (is_master) write(*,*) trim(rname)//&
272       ": prohibit calculations outside CIA T grid?"
273     strictboundcia=.true. ! default value
274     call getin_p("strictboundcia",strictboundcia)
275     if (is_master) write(*,*) trim(rname)//&
276       ": strictboundcia = ",strictboundcia
277
278     if (is_master) write(*,*) trim(rname)//&
279       ": Minimum atmospheric temperature for Planck function integration ?"
280     tplanckmin=30.0 ! default value
281     call getin_p("tplanckmin",tplanckmin)
282     if (is_master) write(*,*) trim(rname)//": tplanckmin = ",tplanckmin
283
284     if (is_master) write(*,*) trim(rname)//&
285       ": Maximum atmospheric temperature for Planck function integration ?"
286     tplanckmax=1500.0 ! default value
287     call getin_p("tplanckmax",tplanckmax)
288     if (is_master) write(*,*) trim(rname)//": tplanckmax = ",tplanckmax
289
290     if (is_master) write(*,*) trim(rname)//&
291       ": Temperature step for Planck function integration ?"
292     dtplanck=0.1 ! default value
293     call getin_p("dtplanck",dtplanck)
294     if (is_master) write(*,*) trim(rname)//": dtplanck = ",dtplanck
295
296     if (is_master) write(*,*) trim(rname)//&
297       ": call gaseous absorption in the visible bands?"//&
298       " (matters only if callrad=T)"
299     callgasvis=.false. ! default value
300     call getin_p("callgasvis",callgasvis)
301     if (is_master) write(*,*) trim(rname)//": callgasvis = ",callgasvis
302
303     if (is_master) write(*,*) trim(rname)//&
304       ": call continuum opacities in radiative transfer ?"//&
305       " (matters only if callrad=T)"
306     continuum=.true. ! default value
307     call getin_p("continuum",continuum)
308     if (is_master) write(*,*) trim(rname)//": continuum = ",continuum
309
310     if (is_master) write(*,*) trim(rname)//&
311       ": call turbulent vertical diffusion ?"
312     calldifv=.true. ! default value
313     call getin_p("calldifv",calldifv)
314     if (is_master) write(*,*) trim(rname)//": calldifv = ",calldifv
315
316     if (is_master) write(*,*) trim(rname)//": use turbdiff instead of vdifc ?"
317     UseTurbDiff=.true. ! default value
318     call getin_p("UseTurbDiff",UseTurbDiff)
319     if (is_master) write(*,*) trim(rname)//": UseTurbDiff = ",UseTurbDiff
320
321     if (is_master) write(*,*) trim(rname)//": call convective adjustment ?"
322     calladj=.true. ! default value
323     call getin_p("calladj",calladj)
324     if (is_master) write(*,*) trim(rname)//": calladj = ",calladj
325
326     if (is_master) write(*,*) trim(rname)//&
327     ": Radiative timescale for Newtonian cooling (in Earth days)?"
328     tau_relax=30. ! default value
329     call getin_p("tau_relax",tau_relax)
330     if (is_master) write(*,*) trim(rname)//": tau_relax = ",tau_relax
331     tau_relax=tau_relax*24*3600 ! convert Earth days --> seconds
332
333     if (is_master) write(*,*)trim(rname)//&
334       ": call thermal conduction in the soil ?"
335     callsoil=.true. ! default value
336     call getin_p("callsoil",callsoil)
337     if (is_master) write(*,*) trim(rname)//": callsoil = ",callsoil
338
339     if (is_master) write(*,*)trim(rname)//&
340       ": Rad transfer is computed every iradia", &
341       " physical timestep"
342     iradia=1 ! default value
343     call getin_p("iradia",iradia)
344     if (is_master) write(*,*)trim(rname)//": iradia = ",iradia
345
346     if (is_master) write(*,*)trim(rname)//": Rayleigh scattering ?"
347     rayleigh=.false.
348     call getin_p("rayleigh",rayleigh)
349     if (is_master) write(*,*)trim(rname)//": rayleigh = ",rayleigh
350
351     if (is_master) write(*,*) trim(rname)//&
352       ": Use blackbody for stellar spectrum ?"
353     stelbbody=.false. ! default value
354     call getin_p("stelbbody",stelbbody)
355     if (is_master) write(*,*) trim(rname)//": stelbbody = ",stelbbody
356
357     if (is_master) write(*,*) trim(rname)//": Stellar blackbody temperature ?"
358     stelTbb=5800.0 ! default value
359     call getin_p("stelTbb",stelTbb)
360     if (is_master) write(*,*) trim(rname)//": stelTbb = ",stelTbb
361
362     if (is_master) write(*,*) trim(rname)//": Output mean OLR in 1D?"
363     meanOLR=.false.
364     call getin_p("meanOLR",meanOLR)
365     if (is_master) write(*,*) trim(rname)//": meanOLR = ",meanOLR
366
367     if (is_master) write(*,*)trim(rname)//": Output spectral OLR in 3D?"
368     specOLR=.false.
369     call getin_p("specOLR",specOLR)
370     if (is_master) write(*,*)trim(rname)//": specOLR = ",specOLR
371
372     if (is_master) write(*,*)trim(rname)//&
373       ": Output diagnostic optical thickness attenuation exp(-tau)?"
374     diagdtau=.false.
375     call getin_p("diagdtau",diagdtau)
376     if (is_master) write(*,*)trim(rname)//": diagdtau = ",diagdtau
377
378     if (is_master) write(*,*)trim(rname)//": Operate in kastprof mode?"
379     kastprof=.false.
380     call getin_p("kastprof",kastprof)
381     if (is_master) write(*,*)trim(rname)//": kastprof = ",kastprof
382
383     if (is_master) write(*,*)trim(rname)//&
384       ": Uniform absorption in radiative transfer?"
385     graybody=.false.
386     call getin_p("graybody",graybody)
387     if (is_master) write(*,*)trim(rname)//": graybody = ",graybody
388
389! Soil model
390     if (is_master) write(*,*)trim(rname)//&
391       ": Number of sub-surface layers for soil scheme?"
392     ! default value of nsoilmx set in comsoil_h
393     call getin_p("nsoilmx",nsoilmx)
394     if (is_master) write(*,*)trim(rname)//": nsoilmx=",nsoilmx
395
396     if (is_master) write(*,*)trim(rname)//&
397       ": Thickness of topmost soil layer (m)?"
398     ! default value of lay1_soil set in comsoil_h
399     call getin_p("lay1_soil",lay1_soil)
400     if (is_master) write(*,*)trim(rname)//": lay1_soil = ",lay1_soil
401
402     if (is_master) write(*,*)trim(rname)//&
403       ": Coefficient for soil layer thickness distribution?"
404     ! default value of alpha_soil set in comsoil_h
405     call getin_p("alpha_soil",alpha_soil)
406     if (is_master) write(*,*)trim(rname)//": alpha_soil = ",alpha_soil
407
408! Chemistry in the thermosphere
409     if (is_master) write(*,*) trim(rname)//": Use deposition ?"
410     depos=.false.         ! default value
411     call getin_p("depos",depos)
412     if (is_master) write(*,*) trim(rname)//": depos = ",depos
413
414     if (is_master) write(*,*)trim(rname)//": Production of haze ?"
415     haze=.false. ! default value
416     call getin_p("haze",haze)
417     if (is_master) write(*,*)trim(rname)//": haze = ",haze
418
419
420
421      if (is_master) write(*,*)trim(rname)// "call thermal conduction ?"
422      callconduct=.false. ! default value
423      call getin_p("callconduct",callconduct)
424      if (is_master) write(*,*)trim(rname)// " callconduct = ",callconduct
425
426      if (is_master) write(*,*)trim(rname)// "call phitop ?"
427      phitop=0. ! default value
428      call getin_p("phitop",phitop)
429      if (is_master) write(*,*)trim(rname)// " phitop = ",phitop
430
431      if (is_master) write(*,*)trim(rname)// "call molecular viscosity ?"
432      callmolvis=.false. ! default value
433      call getin_p("callmolvis",callmolvis)
434      if (is_master) write(*,*)trim(rname)// " callmolvis = ",callmolvis
435
436      if (is_master) write(*,*)trim(rname)// "call molecular diffusion ?"
437      callmoldiff=.false. ! default value
438      call getin_p("callmoldiff",callmoldiff)
439      if (is_master) write(*,*)trim(rname)// " callmoldiff = ",callmoldiff
440
441! Global1D mean and solar zenith angle
442
443     if (ngrid.eq.1) then
444      PRINT*, 'Simulate global averaged conditions ?'
445      global1d = .false. ! default value
446      call getin_p("global1d",global1d)
447      write(*,*) "global1d = ",global1d
448
449      ! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
450      if (global1d.and.diurnal) then
451         call abort_physic(rname,'if global1d is true, diurnal must be set to false',1)
452      endif
453
454      if (global1d) then
455         PRINT *,'Solar Zenith angle (deg.) ?'
456         PRINT *,'(assumed for averaged solar flux S/4)'
457         szangle=60.0  ! default value
458         call getin_p("szangle",szangle)
459         write(*,*) "szangle = ",szangle
460      endif
461   endif ! of if (ngrid.eq.1)
462
463! Test of incompatibility:
464! if kastprof used, we must be in 1D
465     if (kastprof.and.(ngrid.gt.1)) then
466       call abort_physic(rname,'kastprof can only be used in 1D!',1)
467     endif
468
469     if (is_master) write(*,*)trim(rname)//&
470       ": Stratospheric temperature for kastprof mode?"
471     Tstrat=167.0
472     call getin_p("Tstrat",Tstrat)
473     if (is_master) write(*,*)trim(rname)//": Tstrat = ",Tstrat
474
475     if (is_master) write(*,*)trim(rname)//": Remove lower boundary?"
476     nosurf=.false.
477     call getin_p("nosurf",nosurf)
478     if (is_master) write(*,*)trim(rname)//": nosurf = ",nosurf
479
480! Tests of incompatibility:
481     if (nosurf.and.callsoil) then
482       if (is_master) then
483         write(*,*)trim(rname)//'nosurf not compatible with soil scheme!'
484         write(*,*)trim(rname)//'... got to make a choice!'
485       endif
486       call abort_physic(rname,"nosurf not compatible with soil scheme!",1)
487     endif
488
489     if (is_master) write(*,*)trim(rname)//": Add an internal heat flux?", &
490                   "... matters only if callsoil=F"
491     intheat=0.
492     call getin_p("intheat",intheat)
493     if (is_master) write(*,*)trim(rname)//": intheat = ",intheat
494
495
496     if (is_master) write(*,*)trim(rname)//": Test physics timescale in 1D?"
497     testradtimes=.false.
498     call getin_p("testradtimes",testradtimes)
499     if (is_master) write(*,*)trim(rname)//": testradtimes = ",testradtimes
500
501! Test of incompatibility:
502! if testradtimes used, we must be in 1D
503     if (testradtimes.and.(ngrid.gt.1)) then
504       call abort_physic(rname,'testradtimes can only be used in 1D!',1)
505     endif
506
507     if (is_master) write(*,*)trim(rname)//": Default planetary temperature?"
508     tplanet=215.0
509     call getin_p("tplanet",tplanet)
510     if (is_master) write(*,*)trim(rname)//": tplanet = ",tplanet
511
512     if (is_master) write(*,*)trim(rname)//": Which star?"
513     startype=1 ! default value = Sol
514     call getin_p("startype",startype)
515     if (is_master) write(*,*)trim(rname)//": startype = ",startype
516
517     if (is_master) write(*,*)trim(rname)//": Value of stellar flux at 1 AU?"
518     Fat1AU=1356.0 ! default value = Sol today
519     call getin_p("Fat1AU",Fat1AU)
520     if (is_master) write(*,*)trim(rname)//": Fat1AU = ",Fat1AU
521
522
523! TRACERS:
524
525     if (is_master) write(*,*)trim(rname)//": fixed radii for Cloud particles?"
526     radfixed=.false. ! default value
527     call getin_p("radfixed",radfixed)
528     if (is_master) write(*,*)trim(rname)//": radfixed = ",radfixed
529
530     if(kastprof)then
531        radfixed=.true.
532     endif
533
534      if (is_master) write(*,*)trim(rname)//&
535       "Number of radiatively active aerosols:"
536     naerkind=0 ! default value
537     call getin_p("naerkind",naerkind)
538     if (is_master) write(*,*)trim(rname)//": naerkind = ",naerkind
539
540
541!***************************************************************
542         !! TRACERS options
543
544     if (is_master)write(*,*)trim(rname)//&
545      "call N2 condensation ?"
546     n2cond=.true. ! default value
547     call getin_p("n2cond",n2cond)
548     if (is_master)write(*,*)trim(rname)//&
549      " n2cond = ",n2cond
550
551     if (is_master)write(*,*)trim(rname)//&
552      "call N2 cloud condensation ?"
553     condensn2=.false. ! default value
554     call getin_p("condensn2",condensn2)
555     if (is_master)write(*,*)trim(rname)//&
556      "condensn2 = ",condensn2
557
558     if (is_master)write(*,*)trim(rname)//&
559      "call no N2 frost formation?"
560     no_n2frost=.false. ! default value
561     call getin_p("no_n2frost",no_n2frost)
562     if (is_master)write(*,*)trim(rname)//&
563      "no_n2frost = ",no_n2frost
564
565     if (is_master)write(*,*)trim(rname)//&
566      "N2 condensation subtimestep?"
567     nbsub=20 ! default value
568     call getin_p("nbsub",nbsub)
569     if (is_master)write(*,*)trim(rname)//&
570      " nbsub = ",nbsub
571
572     if (is_master)write(*,*)trim(rname)//&
573      "Gravitationnal sedimentation ?"
574     sedimentation=.true. ! default value
575     call getin_p("sedimentation",sedimentation)
576     if (is_master)write(*,*)trim(rname)//&
577      " sedimentation = ",sedimentation
578
579     if (is_master)write(*,*)trim(rname)//&
580      "Compute glacial flow ?"
581     glaflow=.false. ! default value
582     call getin_p("glaflow",glaflow)
583     if (is_master)write(*,*)trim(rname)//&
584      " glaflow = ",glaflow
585
586     if (is_master)write(*,*)trim(rname)//&
587      "Compute methane cycle ?"
588     methane=.false. ! default value
589     call getin_p("methane",methane)
590     if (is_master)write(*,*)trim(rname)//&
591      " methane = ",methane
592     condmetsurf=.true. ! default value
593     call getin_p("condmetsurf",condmetsurf)
594     if (is_master)write(*,*)trim(rname)//&
595      " condmetsurf = ",condmetsurf
596
597     if (is_master)write(*,*)trim(rname)//&
598      "Compute methane clouds ?"
599     metcloud=.false. ! default value
600     call getin_p("metcloud",metcloud)
601     if (is_master)write(*,*)trim(rname)//&
602      " metcloud = ",metcloud
603
604     if (is_master)write(*,*)trim(rname)//&
605      "Compute CO cycle ?"
606     carbox=.false. ! default value
607     call getin_p("carbox",carbox)
608     if (is_master)write(*,*)trim(rname)//&
609      " carbox = ",carbox
610     condcosurf=.true. ! default value
611     call getin_p("condcosurf",condcosurf)
612     if (is_master)write(*,*)trim(rname)//&
613      " condcosurf = ",condcosurf
614
615     if (is_master)write(*,*)trim(rname)//&
616      "Compute CO clouds ?"
617     monoxcloud=.false. ! default value
618     call getin_p("monoxcloud",monoxcloud)
619     if (is_master)write(*,*)trim(rname)//&
620      " monoxcloud = ",monoxcloud
621
622     if (is_master)write(*,*)trim(rname)//&
623     "atmospheric redistribution (s):"
624     tau_n2=1. ! default value
625     call getin_p("tau_n2",tau_n2)
626     if (is_master)write(*,*)trim(rname)//&
627     " tau_n2 = ",tau_n2
628     tau_ch4=1.E7 ! default value
629     call getin_p("tau_ch4",tau_ch4)
630     if (is_master)write(*,*)trim(rname)//&
631     " tau_ch4 = ",tau_ch4
632     tau_co=1. ! default value
633     call getin_p("tau_co",tau_co)
634     if (is_master)write(*,*)trim(rname)//&
635     " tau_co = ",tau_co
636     tau_prechaze=1. ! default value
637     call getin_p("tau_prechaze",tau_prechaze)
638     if (is_master)write(*,*)trim(rname)//&
639     " tau_prechaze = ",tau_prechaze
640
641     if (is_master)write(*,*)trim(rname)//&
642      "Day fraction for limited cold trap in SP?"
643     dayfrac=0. ! default value
644     call getin_p("dayfrac",dayfrac)
645     if (is_master)write(*,*)trim(rname)//&
646      " dayfrac = ",dayfrac
647     thresh_non2=0. ! default value
648     call getin_p("thresh_non2",thresh_non2)
649     if (is_master)write(*,*)trim(rname)//&
650      " thresh_non2 = ",thresh_non2
651
652    ! use Pluto.old routines
653     if (is_master) write(*,*) trim(rname)//": use vdifc from old Pluto ?"
654     oldplutovdifc=.false. ! default value
655     call getin_p("oldplutovdifc",oldplutovdifc)
656     if (is_master) write(*,*) trim(rname)//": oldplutovdifc = ",oldplutovdifc
657
658     if (is_master) write(*,*) trim(rname)//&
659       ": call pluto.old correlated-k radiative transfer ?"
660     oldplutocorrk=.false. ! default value
661     call getin_p("oldplutocorrk",oldplutocorrk)
662     if (is_master) write(*,*) trim(rname)//": oldplutocorrk = ",oldplutocorrk
663
664     if (is_master) write(*,*) trim(rname)//&
665       ": call pluto.old sedimentation ?"
666     oldplutosedim=.false. ! default value
667     call getin_p("oldplutosedim",oldplutosedim)
668     if (is_master) write(*,*) trim(rname)//": oldplutosedim = ",oldplutosedim
669
670!***************************************************************
671     !! Haze options
672
673     if (is_master)write(*,*)trim(rname)//&
674     "Production of haze ?"
675     haze=.false. ! default value
676     call getin_p("haze",haze)
677     if (is_master)write(*,*)trim(rname)//&
678     " haze = ",haze
679     hazeconservch4=.false. ! conservch4, by default value ch4 is photolyzed
680     call getin_p("hazeconservch4",hazeconservch4)
681     if (is_master)write(*,*)trim(rname)//&
682     "hazconservch4 = ",hazeconservch4
683     if (is_master)write(*,*)trim(rname)//&
684     "Production of haze (fast version) ?"
685     fasthaze=.false. ! default value
686     call getin_p("fasthaze",fasthaze)
687     if (is_master)write(*,*)trim(rname)//&
688     "fasthaze = ",fasthaze
689
690     if (is_master)write(*,*)trim(rname)//&
691     "Add source of haze ?"
692     source_haze=.false. ! default value
693     call getin_p("source_haze",source_haze)
694     if (is_master)write(*,*)trim(rname)//&
695     " source_haze = ",source_haze
696     mode_hs=0 ! mode haze source default value
697     call getin_p("mode_hs",mode_hs)
698     if (is_master)write(*,*)trim(rname)//&
699     " mode_hs",mode_hs
700     kfix=1 ! default value
701     call getin_p("kfix",kfix)
702     if (is_master)write(*,*)trim(rname)//&
703     "levels kfix",kfix
704     fracsource=0.1 ! default value
705     call getin_p("fracsource",fracsource)
706     if (is_master)write(*,*)trim(rname)//&
707     " fracsource",fracsource
708     latsource=30. ! default value
709     call getin_p("latsource",latsource)
710     if (is_master)write(*,*)trim(rname)//&
711     " latsource",latsource
712     lonsource=180. ! default value
713     call getin_p("lonsource",lonsource)
714     if (is_master)write(*,*)trim(rname)//&
715     " lonsource",lonsource
716
717     if (is_master)write(*,*)trim(rname)//&
718     "Radiatively active haze ?"
719     aerohaze=.false. ! default value
720     call getin_p("aerohaze",aerohaze)
721     if (is_master)write(*,*)trim(rname)//&
722     "aerohaze = ",aerohaze
723
724     if (is_master)write(*,*)trim(rname)//&
725     "Haze monomer radius ?"
726     rad_haze=10.e-9 ! default value
727     call getin_p("rad_haze",rad_haze)
728     if (is_master)write(*,*)trim(rname)//&
729     "rad_haze = ",rad_haze
730
731     if (is_master)write(*,*)trim(rname)//&
732     "fractal particle ?"
733     fractal=.false. ! default value
734     call getin_p("fractal",fractal)
735     if (is_master)write(*,*)trim(rname)//&
736     "fractal = ",fractal
737     nb_monomer=10 ! default value
738     call getin_p("nb_monomer",nb_monomer)
739     if (is_master)write(*,*)trim(rname)//&
740     "nb_monomer = ",nb_monomer
741
742     if (is_master)write(*,*)trim(rname)//&
743     "fixed gaseous CH4 mixing ratio for rad transfer ?"
744    ch4fix=.false. ! default value
745    call getin_p("ch4fix",ch4fix)
746    if (is_master)write(*,*)trim(rname)//&
747     " ch4fix = ",ch4fix
748    if (is_master)write(*,*)trim(rname)//&
749     "fixed gaseous CH4 mixing ratio for rad transfer ?"
750    vmrch4fix=0.5 ! default value
751    call getin_p("vmrch4fix",vmrch4fix)
752    if (is_master)write(*,*)trim(rname)//&
753     " vmrch4fix = ",vmrch4fix
754    vmrch4_proffix=.false. ! default value
755    call getin_p("vmrch4_proffix",vmrch4_proffix)
756    if (is_master)write(*,*)trim(rname)//&
757     " vmrch4_proffix = ",vmrch4_proffix
758
759    if (is_master)write(*,*)trim(rname)//&
760     "call specific cooling for radiative transfer ?"
761    cooling=.false.  ! default value
762    call getin_p("cooling",cooling)
763    if (is_master)write(*,*)trim(rname)//&
764     " cooling = ",cooling
765
766    if (is_master)write(*,*)trim(rname)//&
767     "NLTE correction for methane heating rates?"
768    nlte=.false.  ! default value
769    call getin_p("nlte",nlte)
770    if (is_master)write(*,*)trim(rname)//&
771     " nlte = ",nlte
772    strobel=.false.  ! default value
773    call getin_p("strobel",strobel)
774    if (is_master)write(*,*)trim(rname)//&
775     " strobel = ",strobel
776
777     if (is_master)write(*,*)trim(rname)//&
778     "fixed radius profile from txt file ?"
779     haze_radproffix=.false. ! default value
780     call getin_p("haze_radproffix",haze_radproffix)
781     if (is_master)write(*,*)trim(rname)//&
782     "haze_radproffix = ",haze_radproffix
783     if (is_master)write(*,*)trim(rname)//&
784     "fixed MMR profile from txt file ?"
785     haze_proffix=.false. ! default value
786     call getin_p("haze_proffix",haze_proffix)
787     if (is_master)write(*,*)trim(rname)//&
788     "haze_proffix = ",haze_proffix
789
790     if (is_master)write(*,*)trim(rname)//&
791     "Number mix ratio of haze particles for co clouds:"
792     Nmix_co=100000. ! default value
793     call getin_p("Nmix_co",Nmix_co)
794     if (is_master)write(*,*)trim(rname)//&
795     " Nmix_co = ",Nmix_co
796
797     if (is_master)write(*,*)trim(rname)//&
798     "Number mix ratio of haze particles for ch4 clouds:"
799     Nmix_ch4=100000. ! default value
800     call getin_p("Nmix_ch4",Nmix_ch4)
801     if (is_master)write(*,*)trim(rname)//&
802     " Nmix_ch4 = ",Nmix_ch4
803
804!***************************************************************
805     !! Surface properties
806
807!*********** N2 *********************************
808
809     if (is_master)write(*,*)trim(rname)//&
810      "Mode for changing N2 albedo"
811     mode_n2=0 ! default value
812     call getin_p("mode_n2",mode_n2)
813     if (is_master)write(*,*)trim(rname)//&
814      " mode_n2 = ",mode_n2
815     thres_n2ice=1. ! default value
816     call getin_p("thres_n2ice",thres_n2ice)
817     if (is_master)write(*,*)trim(rname)//&
818      " thres_n2ice = ",thres_n2ice
819
820     if (is_master)write(*,*)trim(rname)//&
821      "Diff of N2 albedo with thickness"
822     deltab=0. ! default value
823     call getin_p("deltab",deltab)
824     if (is_master)write(*,*)trim(rname)//&
825      " deltab = ",deltab
826
827     if (is_master)write(*,*)trim(rname)//&
828      "albedo N2 beta "
829     alb_n2b=0.7 ! default value
830     call getin_p("alb_n2b",alb_n2b)
831     if (is_master)write(*,*)trim(rname)//&
832      " alb_n2b = ",alb_n2b
833
834     if (is_master)write(*,*)trim(rname)//&
835      "albedo N2 alpha "
836     alb_n2a=0.7 ! default value
837     call getin_p("alb_n2a",alb_n2a)
838     if (is_master)write(*,*)trim(rname)//&
839      " alb_n2a = ",alb_n2a
840
841     if (is_master)write(*,*)trim(rname)//&
842      "emis N2 beta "
843     emis_n2b=0.7 ! default value
844     call getin_p("emis_n2b",emis_n2b)
845     if (is_master)write(*,*)trim(rname)//&
846      " emis_n2b = ",emis_n2b
847
848     if (is_master)write(*,*)trim(rname)//&
849      "emis N2 alpha "
850     emis_n2a=0.7 ! default value
851     call getin_p("emis_n2a",emis_n2a)
852     if (is_master)write(*,*)trim(rname)//&
853      " emis_n2a = ",emis_n2a
854
855!*********** CH4 *********************************
856
857     if (is_master)write(*,*)trim(rname)//&
858      "Mode for changing CH4 albedo"
859     mode_ch4=0 ! default value
860     call getin_p("mode_ch4",mode_ch4)
861     if (is_master)write(*,*)trim(rname)//&
862      " mode_ch4 = ",mode_ch4
863     feedback_met=0 ! default value
864     call getin_p("feedback_met",feedback_met)
865     if (is_master)write(*,*)trim(rname)//&
866      " feedback_met = ",feedback_met
867     thres_ch4ice=1. ! default value, mm
868     call getin_p("thres_ch4ice",thres_ch4ice)
869     if (is_master)write(*,*)trim(rname)//&
870      " thres_ch4ice = ",thres_ch4ice
871     fdch4_latn=-1.
872     fdch4_lats=0.
873     fdch4_lone=0.
874     fdch4_lonw=-1.
875     fdch4_depalb=0.5
876     fdch4_finalb=0.95
877     fdch4_maxalb=0.99
878     fdch4_ampl=200.
879     fdch4_maxice=100.
880     call getin_p("fdch4_latn",fdch4_latn)
881     call getin_p("fdch4_lats",fdch4_lats)
882     call getin_p("fdch4_lone",fdch4_lone)
883     call getin_p("fdch4_lonw",fdch4_lonw)
884     call getin_p("fdch4_depalb",fdch4_depalb)
885     call getin_p("fdch4_finalb",fdch4_finalb)
886     call getin_p("fdch4_maxalb",fdch4_maxalb)
887     call getin_p("fdch4_maxice",fdch4_maxice)
888     call getin_p("fdch4_ampl",fdch4_ampl)
889     if (is_master)write(*,*)trim(rname)//&
890      " Values for albedo feedback = ",fdch4_latn,&
891     fdch4_lats,fdch4_lone,fdch4_lonw,fdch4_depalb,&
892     fdch4_finalb,fdch4_maxalb,fdch4_maxice,fdch4_ampl
893
894     if (is_master)write(*,*)trim(rname)//&
895      "Latitude for diff albedo methane"
896     metlateq=25. ! default value
897     call getin_p("metlateq",metlateq)
898     if (is_master)write(*,*)trim(rname)//&
899      " metlateq = ",metlateq
900
901     if (is_master)write(*,*)trim(rname)//&
902      "Ls1 and Ls2 for change of ch4 albedo"
903     metls1=-1. ! default value
904     metls2=-2. ! default value
905     call getin_p("metls1",metls1)
906     call getin_p("metls2",metls2)
907
908     if (is_master)write(*,*)trim(rname)//&
909      "albedo CH4 "
910     alb_ch4=0.5 ! default value
911     call getin_p("alb_ch4",alb_ch4)
912     if (is_master)write(*,*)trim(rname)//&
913      " alb_ch4 = ",alb_ch4
914
915     if (is_master)write(*,*)trim(rname)//&
916      "albedo equatorial CH4 "
917     alb_ch4_eq=alb_ch4 ! default value
918     call getin_p("alb_ch4_eq",alb_ch4_eq)
919     if (is_master)write(*,*)trim(rname)//&
920      " alb_ch4_eq = ",alb_ch4_eq
921
922     if (is_master)write(*,*)trim(rname)//&
923      "albedo south hemis CH4 "
924     alb_ch4_s=alb_ch4 ! default value
925     call getin_p("alb_ch4_s",alb_ch4_s)
926     if (is_master)write(*,*)trim(rname)//&
927      " alb_ch4_s = ",alb_ch4_s
928
929     if (is_master)write(*,*)trim(rname)//&
930      "emis CH4 "
931     emis_ch4=0.5 ! default value
932     call getin_p("emis_ch4",emis_ch4)
933     if (is_master)write(*,*)trim(rname)//&
934      " emis_ch4 = ",emis_ch4
935
936     if (is_master)write(*,*)trim(rname)//&
937      "CH4 lag for n2 sublimation limitation"
938     ch4lag=.false. ! default value
939     latlag=-90. ! default value
940     vmrlag=1. ! default value
941     call getin_p("ch4lag",ch4lag)
942     call getin_p("latlag",latlag)
943     call getin_p("vmrlag",vmrlag)
944     if (is_master)write(*,*)trim(rname)//&
945      " ch4lag = ",ch4lag
946     if (is_master)write(*,*)trim(rname)//&
947      " latlag = ",latlag
948     if (is_master)write(*,*)trim(rname)//&
949      " vmrlag = ",vmrlag
950
951     if (is_master)write(*,*)trim(rname)//&
952      "Max temperature for surface ?"
953     tsurfmax=.false. ! default value
954     albmin_ch4=0.3 ! default value
955     call getin_p("tsurfmax",tsurfmax)
956     call getin_p("albmin_ch4",albmin_ch4)
957     if (is_master)write(*,*)trim(rname)//&
958      " tsurfmax = ",tsurfmax
959     if (is_master)write(*,*)trim(rname)//&
960      " albmin_ch4 = ",albmin_ch4
961
962
963     if (is_master)write(*,*)trim(rname)//&
964     "fixed gaseous CH4 mixing ratio for rad transfer ?"
965     ch4fix=.false. ! default value
966     call getin_p("ch4fix",ch4fix)
967     if (is_master)write(*,*)trim(rname)//&
968       " ch4fix = ",ch4fix
969     if (is_master)write(*,*)trim(rname)//&
970       "fixed gaseous CH4 mixing ratio for rad transfer ?"
971     vmrch4fix=0.5 ! default value
972     call getin_p("vmrch4fix",vmrch4fix)
973     if (is_master)write(*,*)trim(rname)//&
974        " vmrch4fix = ",vmrch4fix
975     vmrch4_proffix=.false. ! default value
976     call getin_p("vmrch4_proffix",vmrch4_proffix)
977     if (is_master)write(*,*)trim(rname)//&
978        " vmrch4_proffix = ",vmrch4_proffix
979
980
981!*********** CO *********************************
982
983     if (is_master)write(*,*)trim(rname)//&
984      "albedo CO "
985     alb_co=0.4 ! default value
986     call getin_p("alb_co",alb_co)
987     if (is_master)write(*,*)trim(rname)//&
988      " alb_co = ",alb_co
989     thres_coice=1. ! default value, mm
990     call getin_p("thres_coice",thres_coice)
991     if (is_master)write(*,*)trim(rname)//&
992      " thres_coice = ",thres_coice
993
994     if (is_master)write(*,*)trim(rname)//&
995      "emis CO "
996     emis_co=0.4 ! default value
997     call getin_p("emis_co",emis_co)
998     if (is_master)write(*,*)trim(rname)//&
999      " emis_co = ",emis_co
1000
1001!*********** THOLINS *********************************
1002     if (is_master)write(*,*)trim(rname)//&
1003      "Mode for changing tholins albedo/emis"
1004     mode_tholins=0 ! default value
1005     call getin_p("mode_tholins",mode_tholins)
1006     if (is_master)write(*,*)trim(rname)//&
1007      " mode_tholins = ",mode_tholins
1008
1009     if (is_master)write(*,*)trim(rname)//&
1010      "albedo tho "
1011     alb_tho=0.1 ! default value
1012     call getin_p("alb_tho",alb_tho)
1013     if (is_master)write(*,*)trim(rname)//&
1014      " alb_tho = ",alb_tho
1015
1016     if (is_master)write(*,*)trim(rname)//&
1017      "albedo tho eq"
1018     alb_tho_eq=0.1 ! default value
1019     call getin_p("alb_tho_eq",alb_tho_eq)
1020     if (is_master)write(*,*)trim(rname)//&
1021      " alb_tho_eq = ",alb_tho_eq
1022
1023     if (is_master)write(*,*)trim(rname)//&
1024      "emis tholins "
1025     emis_tho=1. ! default value
1026     call getin_p("emis_tho",emis_tho)
1027     if (is_master)write(*,*)trim(rname)//&
1028      " emis_tho = ",emis_tho
1029
1030     if (is_master)write(*,*)trim(rname)//&
1031      "emis tholins eq"
1032     emis_tho_eq=1. ! default value
1033     call getin_p("emis_tho_eq",emis_tho_eq)
1034     if (is_master)write(*,*)trim(rname)//&
1035      " emis_tho_eq = ",emis_tho_eq
1036
1037     if (is_master)write(*,*)trim(rname)//&
1038      "Latitude for diff albedo tholins"
1039     tholateq=25. ! default value
1040     call getin_p("tholateq",tholateq)
1041     if (is_master)write(*,*)trim(rname)//&
1042      " tholateq = ",tholateq
1043     tholatn=-1.
1044     tholats=0.
1045     tholone=0.
1046     tholonw=-1.
1047     alb_tho_spe=0.1 ! default value
1048     emis_tho_spe=1. ! default value
1049     call getin_p("  tholatn",tholatn)
1050     call getin_p("  tholats",tholats)
1051     call getin_p("  tholone",tholone)
1052     call getin_p("  tholonw",tholonw)
1053     if (is_master)write(*,*)trim(rname)//&
1054      " Values for special tholins albedo = ",tholatn,&
1055        tholats,tholone,tholonw,alb_tho_spe,emis_tho_spe
1056
1057     if (is_master)write(*,*)trim(rname)//&
1058      "Specific albedo"
1059     spelon1=-180. ! default value
1060     spelon2=180. ! default value
1061     spelat1=-90. ! default value
1062     spelat2=90. ! default value
1063     specalb=.false. ! default value
1064     if (is_master)write(*,*)trim(rname)//&
1065      "albedo/emis spe "
1066     albspe=0.1 ! default value
1067     emispe=1. ! default value
1068     call getin_p("spelon1",spelon1)
1069     call getin_p("spelon2",spelon2)
1070     call getin_p("spelat1",spelat1)
1071     call getin_p("spelat2",spelat2)
1072     call getin_p("specalb",specalb)
1073     call getin_p("albspe",albspe)
1074     call getin_p("emispe",emispe)
1075
1076     if (is_master)write(*,*)trim(rname)//&
1077      " specific = ",specalb
1078
1079!********************** TI *****************************
1080
1081     if (is_master)write(*,*)trim(rname)//&
1082      "Change TI with time"
1083     changeti=.false. ! default value
1084     call getin_p("changeti",changeti)
1085     if (is_master)write(*,*)trim(rname)//&
1086      " changeti = ",changeti
1087     changetid=.false. ! default value for diurnal TI
1088     call getin_p("changetid",changetid)
1089     if (is_master)write(*,*)trim(rname)//&
1090      " changetid = ",changetid
1091
1092     if (is_master)write(*,*)trim(rname)//&
1093      "IT N2 "
1094     ITN2=800. ! default value
1095     call getin_p("ITN2",ITN2)
1096     if (is_master)write(*,*)trim(rname)//&
1097      " ITN2 = ",ITN2
1098     ITN2d=20. ! default value
1099     call getin_p("ITN2d",ITN2d)
1100     if (is_master)write(*,*)trim(rname)//&
1101      " ITN2d = ",ITN2d
1102
1103     if (is_master)write(*,*)trim(rname)//&
1104      "IT CH4"
1105     ITCH4=800. ! default value
1106     call getin_p("ITCH4",ITCH4)
1107     if (is_master)write(*,*)trim(rname)//&
1108      " ITCH4 = ",ITCH4
1109     ITCH4d=20. ! default value
1110     call getin_p("ITCH4d",ITCH4d)
1111     if (is_master)write(*,*)trim(rname)//&
1112      " ITCH4d = ",ITCH4d
1113
1114     if (is_master)write(*,*)trim(rname)//&
1115      "IT H2O"
1116     ITH2O=800. ! default value
1117     call getin_p("ITH2O",ITH2O)
1118     if (is_master)write(*,*)trim(rname)//&
1119      " ITH2O = ",ITH2O
1120     ITH2Od=20. ! default value
1121     call getin_p("ITH2Od",ITH2Od)
1122     if (is_master)write(*,*)trim(rname)//&
1123      " ITH2Od = ",ITH2Od
1124
1125!************** COOLING ***************
1126
1127     alpha_top=5e-11 ! default value
1128     call getin_p("alpha_top",alpha_top)
1129     if (is_master)write(*,*)trim(rname)//&
1130      " alpha_top = ",alpha_top
1131     pref=0.02 ! default value
1132     call getin_p("pref",pref)
1133     if (is_master)write(*,*)trim(rname)//&
1134      " pref = ",pref
1135     deltap=0.1 ! default value
1136     call getin_p("deltap",deltap)
1137     if (is_master)write(*,*)trim(rname)//&
1138      " deltap = ",deltap
1139
1140!=================================
1141! TWOLAY scheme and NH3 cloudare left for retrocompatibility only,
1142! You should now use N-LAYER scheme (see below).
1143
1144     if (is_master) write(*,*)trim(rname)//&
1145       ": Radiatively active two-layer aerosols?"
1146     aeroback2lay=.false.     ! default value
1147     call getin_p("aeroback2lay",aeroback2lay)
1148     if (is_master) write(*,*)trim(rname)//": aeroback2lay = ",aeroback2lay
1149
1150     if (aeroback2lay.and.is_master) then
1151       print*,'Warning : The TWOLAY AEROSOL scheme is deprecated and buggy...'
1152       print*,'You should use the generic n-layer scheme (see aeronlay).'
1153     endif
1154
1155     if (is_master) write(*,*)trim(rname)//": Radiatively active ammonia cloud?"
1156     aeronh3=.false.     ! default value
1157     call getin_p("aeronh3",aeronh3)
1158     if (is_master) write(*,*)trim(rname)//": aeronh3 = ",aeronh3
1159
1160     if (aeronh3.and.is_master) then
1161       print*,'Warning : You are using specific NH3 cloud scheme ...'
1162       print*,'You should use the generic n-layer scheme (see aeronlay).'
1163     endif
1164
1165     if (is_master) write(*,*)trim(rname)//&
1166       ": TWOLAY AEROSOL: total optical depth "//&
1167       "in the tropospheric layer (visible)"
1168     obs_tau_col_tropo=8.D0
1169     call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
1170     if (is_master) write(*,*)trim(rname)//&
1171       ": obs_tau_col_tropo = ",obs_tau_col_tropo
1172
1173     if (is_master) write(*,*)trim(rname)//&
1174       ": TWOLAY AEROSOL: total optical depth "//&
1175       "in the stratospheric layer (visible)"
1176     obs_tau_col_strato=0.08D0
1177     call getin_p("obs_tau_col_strato",obs_tau_col_strato)
1178     if (is_master) write(*,*)trim(rname)//&
1179       ": obs_tau_col_strato = ",obs_tau_col_strato
1180
1181     if (is_master) write(*,*)trim(rname)//&
1182       ": TWOLAY AEROSOL: optprop_back2lay_vis?"
1183     optprop_back2lay_vis = 'optprop_saturn_vis_n20.dat'
1184     call getin_p("optprop_back2lay_vis",optprop_back2lay_vis)
1185     if (is_master) write(*,*)trim(rname)//&
1186       ": optprop_back2lay_vis = ",trim(optprop_back2lay_vis)
1187
1188     if (is_master) write(*,*)trim(rname)//&
1189       ": TWOLAY AEROSOL: optprop_back2lay_ir?"
1190     optprop_back2lay_ir = 'optprop_saturn_ir_n20.dat'
1191     call getin_p("optprop_back2lay_ir",optprop_back2lay_ir)
1192     if (is_master) write(*,*)trim(rname)//&
1193       ": optprop_back2lay_ir = ",trim(optprop_back2lay_ir)
1194
1195     if (is_master) write(*,*)trim(rname)//&
1196       ": TWOLAY AEROSOL: pres_bottom_tropo? in pa"
1197     pres_bottom_tropo=66000.0
1198     call getin_p("pres_bottom_tropo",pres_bottom_tropo)
1199     if (is_master) write(*,*)trim(rname)//&
1200       ": pres_bottom_tropo = ",pres_bottom_tropo
1201
1202     if (is_master) write(*,*)trim(rname)//&
1203       ": TWOLAY AEROSOL: pres_top_tropo? in pa"
1204     pres_top_tropo=18000.0
1205     call getin_p("pres_top_tropo",pres_top_tropo)
1206     if (is_master) write(*,*)trim(rname)//&
1207       ": pres_top_tropo = ",pres_top_tropo
1208
1209     if (is_master) write(*,*)trim(rname)//&
1210       ": TWOLAY AEROSOL: pres_bottom_strato? in pa"
1211     pres_bottom_strato=2000.0
1212     call getin_p("pres_bottom_strato",pres_bottom_strato)
1213     if (is_master) write(*,*)trim(rname)//&
1214       ": pres_bottom_strato = ",pres_bottom_strato
1215
1216     ! Sanity check
1217     if (pres_bottom_strato .gt. pres_top_tropo) then
1218       if(is_master) then
1219       print*,'Error : TWOLAY AEROSOL, Please ensure that in callphys.def'
1220       print*,'you have pres_top_tropo > pres_bottom_strato !'
1221       endif
1222       call abort_physic(rname," pres_top_tropo > pres_bottom_strato!",1)
1223     endif
1224
1225     if (is_master) write(*,*)trim(rname)//&
1226       ": TWOLAY AEROSOL: pres_top_strato? in pa"
1227     pres_top_strato=100.0
1228     call getin_p("pres_top_strato",pres_top_strato)
1229     if (is_master) write(*,*)trim(rname)//&
1230       ": pres_top_strato = ",pres_top_strato
1231
1232     if (is_master) write(*,*)trim(rname)//&
1233       ": TWOLAY AEROSOL: particle size in the ", &
1234       "tropospheric layer, in meters"
1235     size_tropo=2.e-6
1236     call getin_p("size_tropo",size_tropo)
1237     if (is_master) write(*,*)trim(rname)//": size_tropo = ",size_tropo
1238
1239     if (is_master) write(*,*)trim(rname)//&
1240       ": TWOLAY AEROSOL: particle size in the ", &
1241       "stratospheric layer, in meters"
1242     size_strato=1.e-7
1243     call getin_p("size_strato",size_strato)
1244     if (is_master) write(*,*)trim(rname)//": size_strato = ",size_strato
1245
1246     if (is_master) write(*,*)trim(rname)//&
1247       ": NH3 (thin) cloud: total optical depth"
1248     tau_nh3_cloud=7.D0
1249     call getin_p("tau_nh3_cloud",tau_nh3_cloud)
1250     if (is_master) write(*,*)trim(rname)//": tau_nh3_cloud = ",tau_nh3_cloud
1251
1252     if (is_master) write(*,*)trim(rname)//": NH3 (thin) cloud pressure level"
1253     pres_nh3_cloud=7.D0
1254     call getin_p("pres_nh3_cloud",pres_nh3_cloud)
1255     if (is_master) write(*,*)trim(rname)//": pres_nh3_cloud = ",pres_nh3_cloud
1256
1257     if (is_master) write(*,*)trim(rname)//": NH3 (thin) cloud: particle sizes"
1258     size_nh3_cloud=3.e-6
1259     call getin_p("size_nh3_cloud",size_nh3_cloud)
1260     if (is_master) write(*,*)trim(rname)//": size_nh3_cloud = ",size_nh3_cloud
1261
1262!=================================
1263! Generic N-LAYER aerosol scheme
1264
1265     if (is_master) write(*,*)trim(rname)//&
1266       ": Radiatively active generic n-layer aerosols?"
1267     aeronlay=.false.     ! default value
1268     call getin_p("aeronlay",aeronlay)
1269     if (is_master) write(*,*)trim(rname)//": aeronlay = ",aeronlay
1270
1271     if (is_master) write(*,*)trim(rname)//&
1272       ": Number of generic aerosols layers?"
1273     nlayaero=1     ! default value
1274     call getin_p("nlayaero",nlayaero)
1275     ! Avoid to allocate arrays of size 0
1276     if (aeronlay .and. nlayaero.lt.1) then
1277       if (is_master) then
1278       print*, " You are trying to set no generic aerosols..."
1279       print*, " Set aeronlay=.false. instead ! I abort."
1280       endif
1281       call abort_physic(rname,"no generic aerosols...",1)
1282     endif
1283     if (.not. aeronlay) nlayaero=1
1284     if (is_master) write(*,*)trim(rname)//": nlayaero = ",nlayaero
1285
1286     ! This is necessary, we just set the number of aerosol layers
1287     IF(.NOT.ALLOCATED(aeronlay_tauref))      ALLOCATE(aeronlay_tauref(nlayaero))
1288     IF(.NOT.ALLOCATED(aeronlay_lamref))      ALLOCATE(aeronlay_lamref(nlayaero))
1289     IF(.NOT.ALLOCATED(aeronlay_choice))      ALLOCATE(aeronlay_choice(nlayaero))
1290     IF(.NOT.ALLOCATED(aeronlay_pbot))        ALLOCATE(aeronlay_pbot(nlayaero))
1291     IF(.NOT.ALLOCATED(aeronlay_ptop))        ALLOCATE(aeronlay_ptop(nlayaero))
1292     IF(.NOT.ALLOCATED(aeronlay_sclhght))     ALLOCATE(aeronlay_sclhght(nlayaero))
1293     IF(.NOT.ALLOCATED(aeronlay_size))        ALLOCATE(aeronlay_size(nlayaero))
1294     IF(.NOT.ALLOCATED(aeronlay_nueff))       ALLOCATE(aeronlay_nueff(nlayaero))
1295     IF(.NOT.ALLOCATED(optprop_aeronlay_ir))  ALLOCATE(optprop_aeronlay_ir(nlayaero))
1296     IF(.NOT.ALLOCATED(optprop_aeronlay_vis)) ALLOCATE(optprop_aeronlay_vis(nlayaero))
1297
1298     if (is_master) write(*,*)trim(rname)//&
1299       ": Generic n-layer aerosols: Optical depth at reference wavelenght"
1300     aeronlay_tauref=1.0E-1
1301     call getin_p("aeronlay_tauref",aeronlay_tauref)
1302     if (is_master) write(*,*)trim(rname)//&
1303       ": aeronlay_tauref = ",aeronlay_tauref
1304
1305     if (is_master) write(*,*)trim(rname)//&
1306       ": Generic n-layer aerosols: Reference wavelenght for optical depths (m)"
1307     aeronlay_lamref=0.6E-6
1308     call getin_p("aeronlay_lamref",aeronlay_lamref)
1309     if (is_master) write(*,*)trim(rname)//&
1310       ": aeronlay_lamref = ",aeronlay_lamref
1311
1312     if (is_master) then
1313       write(*,*)trim(rname)//&
1314       ": Generic n-layer aerosols: Vertical profile choice : "
1315       write(*,*)trim(rname)//&
1316       "             (1) Tau btwn ptop and pbot follows atm. scale height"
1317       write(*,*)trim(rname)//&
1318       "         or  (2) Tau above pbot follows its own scale height"
1319     endif
1320     aeronlay_choice=1
1321     call getin_p("aeronlay_choice",aeronlay_choice)
1322     if (is_master) write(*,*)trim(rname)//&
1323       ": aeronlay_choice = ",aeronlay_choice
1324
1325     if (is_master) write(*,*)trim(rname)//&
1326       ": Generic n-layer aerosols: bottom pressures (Pa)"
1327     aeronlay_pbot=2000.0
1328     call getin_p("aeronlay_pbot",aeronlay_pbot)
1329     if (is_master) write(*,*)trim(rname)//": aeronlay_pbot = ",aeronlay_pbot
1330
1331     if (is_master) write(*,*)trim(rname)//&
1332       ": Generic n-layer aerosols: (if choice=1) Top pressures (Pa) "
1333     aeronlay_ptop=300000.0
1334     call getin_p("aeronlay_ptop",aeronlay_ptop)
1335     if (is_master) write(*,*)trim(rname)//": aeronlay_ptop = ",aeronlay_ptop
1336
1337     if (is_master) write(*,*)trim(rname)//&
1338       ": Generic n-layer aerosols: (if choice=2) Scale height / atm. scale height"
1339     aeronlay_sclhght=0.2
1340     call getin_p("aeronlay_sclhght",aeronlay_sclhght)
1341     if (is_master) write(*,*)trim(rname)//&
1342       ": aeronlay_sclhght = ",aeronlay_sclhght
1343
1344     if (is_master) write(*,*)trim(rname)//&
1345       ": Generic n-layer aerosols: particles effective radii (m)"
1346     aeronlay_size=1.e-6
1347     call getin_p("aeronlay_size",aeronlay_size)
1348     if (is_master) write(*,*)trim(rname)//": aeronlay_size = ",aeronlay_size
1349
1350     if (is_master) write(*,*)trim(rname)//&
1351       ": Generic n-layer aerosols: particles radii effective variance"
1352     aeronlay_nueff=0.1
1353     call getin_p("aeronlay_nueff",aeronlay_nueff)
1354     if (is_master) write(*,*)trim(rname)//": aeronlay_nueff = ",aeronlay_nueff
1355
1356     if (is_master) write(*,*)trim(rname)//&
1357       ": Generic n-layer aerosols: VIS optical properties file"
1358     optprop_aeronlay_vis = 'optprop_saturn_vis_n20.dat'
1359     call getin_p("optprop_aeronlay_vis",optprop_aeronlay_vis)
1360     if (is_master) write(*,*)trim(rname)//&
1361       ": optprop_aeronlay_vis = ",optprop_aeronlay_vis
1362
1363     if (is_master) write(*,*)trim(rname)//&
1364     ": Generic n-layer aerosols: IR optical properties file"
1365     optprop_aeronlay_ir = 'optprop_saturn_ir_n20.dat'
1366     call getin_p("optprop_aeronlay_ir",optprop_aeronlay_ir)
1367     if (is_master) write(*,*)trim(rname)//&
1368       ": optprop_aeronlay_ir = ",optprop_aeronlay_ir
1369
1370
1371!=================================
1372
1373     if (is_master) write(*,*)trim(rname)//&
1374       ": Is the variable gas species radiatively active?"
1375     Tstrat=167.0
1376     varactive=.false.
1377     call getin_p("varactive",varactive)
1378     if (is_master) write(*,*)trim(rname)//": varactive = ",varactive
1379
1380     if (is_master) write(*,*)trim(rname)//&
1381       ": Is the variable gas species distribution set?"
1382     varfixed=.false.
1383     call getin_p("varfixed",varfixed)
1384     if (is_master) write(*,*)trim(rname)//": varfixed = ",varfixed
1385
1386     if (is_master) write(*,*)trim(rname)//&
1387       ": What is the saturation % of the variable species?"
1388     satval=0.8
1389     call getin_p("satval",satval)
1390     if (is_master) write(*,*)trim(rname)//": satval = ",satval
1391
1392
1393! Test of incompatibility:
1394! if varactive, then varfixed should be false
1395     if (varactive.and.varfixed) then
1396       call abort_physic(rname,'if varactive, varfixed must be OFF!',1)
1397     endif
1398
1399     if (is_master) write(*,*)trim(rname)//": Gravitationnal sedimentation ?"
1400     sedimentation=.false. ! default value
1401     call getin_p("sedimentation",sedimentation)
1402     if (is_master) write(*,*)trim(rname)//": sedimentation = ",sedimentation
1403
1404     if (is_master) write(*,*)trim(rname)//": Generic Condensation of tracers ?"
1405     generic_condensation=.false. !default value
1406     call getin_p("generic_condensation",generic_condensation)
1407     if (is_master) write(*,*)trim(rname)//": generic_condensation = ",generic_condensation
1408
1409     if (is_master) write(*,*)trim(rname)//": Spectral Dependant albedo ?"
1410     albedo_spectral_mode=.false. ! default value
1411     call getin_p("albedo_spectral_mode",albedo_spectral_mode)
1412     if (is_master) write(*,*)trim(rname)//&
1413     ": albedo_spectral_mode = ",albedo_spectral_mode
1414
1415     if (is_master) then
1416       write(*,*)trim(rname)//": Snow albedo ?"
1417       write(*,*)trim(rname)//": If albedo_spectral_mode=.true., then this "
1418       write(*,*)trim(rname)//": corresponds to the 0.5 microns snow albedo."
1419     endif
1420     albedosnow=0.5         ! default value
1421     call getin_p("albedosnow",albedosnow)
1422     if (is_master) write(*,*)trim(rname)//": albedosnow = ",albedosnow
1423
1424     if (is_master) write(*,*)trim(rname)//": N2 ice albedo ?"
1425     albedon2ice=0.5       ! default value
1426     call getin_p("albedon2ice",albedon2ice)
1427     if (is_master) write(*,*)trim(rname)//": albedon2ice = ",albedon2ice
1428
1429     if (is_master) write(*,*)trim(rname)//": Maximum ice thickness ?"
1430     maxicethick=2.0         ! default value
1431     call getin_p("maxicethick",maxicethick)
1432     if (is_master) write(*,*)trim(rname)//": maxicethick = ",maxicethick
1433
1434     if (is_master) write(*,*)trim(rname)//": Minimum eddy mix coeff in 1D ?"
1435     kmixmin=1.0e-2         ! default value
1436     call getin_p("kmixmin",kmixmin)
1437     if (is_master) write(*,*)trim(rname)//": kmixmin = ",kmixmin
1438
1439     if (is_master) write(*,*)'Predefined Cp from dynamics is ',cpp,'J kg^-1 K^-1'
1440     if (is_master) write(*,*)'Predefined Mg from dynamics is ',mugaz,'amu'
1441
1442     force_cpp=.false. ! default value
1443     call getin_p("force_cpp",force_cpp)
1444     if (force_cpp) then
1445      if (is_master) write(*,*)trim(rname)//": force_cpp = ",force_cpp
1446      if (is_master) write(*,*)trim(rname)//": force_cpp is deprecated.",&
1447      "Set cpp_mugaz_mode=1 in callfis to emulate force_cpp=.true."
1448      call abort_physic(rname,"Anyway, you need to set force_cpp=.false. to continue.",1)
1449     endif
1450
1451     if (is_master) write(*,*)trim(rname)//&
1452     ": where do you want your cpp/mugaz value to come from?",&
1453     "=> 0: dynamics (3d), 1: forced in callfis (1d), 2: computed from gases.def (1d)?"
1454     cpp_mugaz_mode = 0 ! default value
1455     call getin_p("cpp_mugaz_mode",cpp_mugaz_mode)
1456     if (is_master) write(*,*)trim(rname)//": cpp_mugaz_mode = ",cpp_mugaz_mode
1457
1458
1459! Test of incompatibility:
1460
1461     if ((.not.tracer).and.(haze)) then
1462       call abort_physic(rname, 'if haze are on, tracers must be on!', 1)
1463     endif
1464     if (haze_proffix.and.sedimentation) then
1465         call abort_physic(rname, 'if haze profile is set, sedimentation must be deactivated', 1)
1466     endif
1467     if (callmolvis.and..not.callconduct) then
1468         call abort_physic(rname, 'if callmolvis is set, callconduct must be true', 1)
1469     endif
1470     if (glaflow.and..not.fast) then
1471         call abort_physic(rname, 'if glaflow is set, fast must be true', 1)
1472     endif
1473     if (paleo.and..not.fast) then
1474         call abort_physic(rname, 'if paleo is set, fast must be true', 1)
1475     endif
1476     if ((haze_proffix.or.haze_radproffix).and..not.aerohaze) then
1477         call abort_physic(rname, 'for now, haze/rad proffix only works w aerohaze=T', 1)
1478     endif
1479      if (condcosurf.and.no_n2frost) then
1480        call abort_physic(rname, "CO surface condensation and no_n2frost are both active which may not be relevant", 1)
1481      end if
1482
1483     if ((cpp_mugaz_mode >= 1).and.(is_master).and.(ngrid>1)) then
1484        write(*,*)'    !!!!! Be aware that having different values for cpp and mugaz in the dynamics and physics'
1485        write(*,*)'    in 3D can result in very pathological behavior. You have been warned !!!!!'
1486      else if ((cpp_mugaz_mode == 0).and.(is_master).and.(ngrid==1)) then
1487        ! for this specific 1D error we will remove run.def before aborting JL22
1488        call system("rm -rf run.def")
1489        call abort_physic(rname,"cpp_mugaz_mode must be >= 1 in 1d",1)
1490     endif
1491
1492     if (cpp_mugaz_mode == 1) then
1493       mugaz = -99999.
1494       if (is_master) write(*,*)trim(rname)//&
1495         ": MEAN MOLECULAR MASS in g mol-1 ?"
1496       call getin_p("mugaz",mugaz)
1497       IF (mugaz.eq.-99999.) THEN
1498         call abort_physic(rname,"mugaz must be set if cpp_mugaz_mode = 1",1)
1499       ENDIF
1500       cpp = -99999.
1501       if (is_master) write(*,*)trim(rname)//&
1502         ": SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
1503       call getin_p("cpp",cpp)
1504       IF (cpp.eq.-99999.) THEN
1505           PRINT *, "cpp must be set if cpp_mugaz_mode = 1"
1506           STOP
1507       ENDIF
1508       if (is_master) write(*,*)'New Cp from callfis is ',cpp,'J kg^-1 K^-1'
1509       if (is_master) write(*,*)'New Mg from callfis is ',mugaz,'amu'
1510
1511     endif ! of if (cpp_mugaz_mode == 1)
1512     call su_gases
1513     call calc_cpp_mugaz
1514
1515     if (is_master) then
1516       PRINT*,'--------------------------------------------'
1517       PRINT*
1518       PRINT*
1519     endif
1520  ELSE
1521     call abort_physic(rname,'Cannot read file callphys.def. Is it here ?',1)
1522  ENDIF ! of IF(iscallphys)
1523
1524  if (is_master) then
1525    PRINT*
1526    PRINT*,'inifis: daysec',daysec
1527    PRINT*
1528    PRINT*,'inifis: The radiative transfer is computed:'
1529    PRINT*,'           each ',iradia,' physical time-step'
1530    PRINT*,'        or each ',iradia*dtphys,' seconds'
1531    PRINT*
1532  endif
1533
1534!-----------------------------------------------------------------------
1535!     Some more initialization:
1536!     ------------------------
1537
1538  ! Initializations for comgeomfi_h
1539#ifndef MESOSCALE
1540  totarea=SSUM(ngrid,parea,1)
1541  call planetwide_sumval(parea,totarea_planet)
1542
1543  !! those are defined in comdiurn_h.F90
1544  IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
1545  IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
1546  IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
1547  IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
1548
1549  DO ig=1,ngrid
1550     sinlat(ig)=sin(plat(ig))
1551     coslat(ig)=cos(plat(ig))
1552     sinlon(ig)=sin(plon(ig))
1553     coslon(ig)=cos(plon(ig))
1554  ENDDO
1555#endif
1556  ! initialize variables in radinc_h
1557  call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
1558
1559  ! initialize variables and allocate arrays in radcommon_h
1560  call ini_radcommon_h(naerkind)
1561
1562  ! allocate "comsoil_h" arrays
1563  call ini_comsoil_h(ngrid)
1564
1565  END SUBROUTINE inifis
1566
1567END MODULE inifis_mod
Note: See TracBrowser for help on using the repository browser.