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

Last change on this file since 3754 was 3754, checked in by afalco, 7 weeks ago

Pluto: imported orographic gravity waves from Mars.
AF

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