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

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

Pluto: ecritphy changed into diagfi_output_rate (followup of generic model), which defines the output rate of the diagfi in physical timesteps rather than dynamical ones.
AF

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