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

Last change on this file since 3380 was 3378, checked in by afalco, 5 months ago

Pluto PCM: changed default params to match aerosol data. AF

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