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

Last change on this file since 3684 was 3684, checked in by debatzbr, 3 months ago

Cleaning for new diagnostics of optical thickness

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