source: trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90 @ 1682

Last change on this file since 1682 was 1682, checked in by emillour, 8 years ago

All GCMs: set things up to enable pluging physics with dynamico

  • dyn3d
  • gcm.F90 : move I/O initialization (dates) to be done before physics

initialization

  • dyn3dpar
  • gcm.F : move I/O initialization (dates) to be done before physics

initialization

  • dynphy_lonlat:
  • inigeomphy_mod.F90 : add ind_cell_glo computation and transfer

to init_geometry

  • phy_common:
  • geometry_mod.F90 : add ind_cell_glo module variable to store global

column index

  • print_control_mod.F90 : make initialization occur via init_print_control_mod

to avoid circular module dependencies

  • init_print_control_mod.F90 : added to initialize print_control_mod module

variables

  • mod_phys_lmdz_mpi_data.F90 : use print_control_mod (rather than iniprint.h)
  • mod_phys_lmdz_para.F90 : use print_control_mod (rather than iniprint.h)
  • mod_phys_lmdz_omp_data.F90 : add is_omp_master (alias of is_omp_root) module

variable and use print_control_mod (rather than
iniprint.h)

  • physics_distribution_mod.F90 : add call to init_dimphy in

init_physics_distribution

  • xios_writefield.F90 : generic routine to output field with XIOS (for debug)
  • misc:
  • handle_err_m.F90 : call abort_physic, rather than abort_gcm
  • wxios.F90 : updates to enable unstructured grids

set module variable g_ctx_name to "LMDZ"
wxios_init(): remove call to wxios_context_init
wxios_context_init(): call xios_context_initialize with COMM_LMDZ_PHY
add routine wxios_set_context() to get handle and set context to XIOS
wxios_domain_param(): change arguments and generate the domain in-place
add wxios_domain_param_unstructured(): generate domain for unstructured case

NB: access is via "domain group" (whereas it is via "domain" in

wxios_domain_param)

  • dynphy_lonlat/phy[std|mars|venus|titan]:
  • iniphysiq_mod.F90 : Remove call to init_dimphy (which is now done in

phy_common/physics_distribution_mod.F90)

EM
File size: 24.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 datafile_mod, only: datadir
15  use comdiurn_h, only: sinlat, coslat, sinlon, coslon
16  use comgeomfi_h, only: totarea, totarea_planet
17  use comsoil_h, only: ini_comsoil_h, nsoilmx, lay1_soil, alpha_soil
18  use time_phylmdz_mod, only: ecritphy,day_step,iphysiq, &
19                              init_time, daysec, dtphys
20  use comcstfi_mod, only: rad, cpp, g, r, rcp, &
21                          mugaz, pi, avocado
22  use planete_mod, only: nres
23  use planetwide_mod, only: planetwide_sumval
24  use callkeys_mod
25  use mod_phys_lmdz_para, only : is_parallel
26
27!=======================================================================
28!
29!   purpose:
30!   -------
31!
32!   Initialisation for the physical parametrisations of the LMD
33!   Generic Model.
34!
35!   author: Frederic Hourdin 15 / 10 /93
36!   -------
37!   modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
38!             Ehouarn Millour (oct. 2008) tracers are now identified
39!              by their names and may not be contiguously
40!              stored in the q(:,:,:,:) array
41!             E.M. (june 2009) use getin routine to load parameters
42!
43!
44!   arguments:
45!   ----------
46!
47!   input:
48!   ------
49!
50!    ngrid                 Size of the horizontal grid.
51!                          All internal loops are performed on that grid.
52!    nlayer                Number of vertical layers.
53!    pdayref               Day of reference for the simulation
54!    pday                  Number of days counted from the North. Spring
55!                          equinoxe.
56!
57!=======================================================================
58!
59!-----------------------------------------------------------------------
60!   declarations:
61!   -------------
62  use datafile_mod, only: datadir
63  use ioipsl_getin_p_mod, only: getin_p
64  IMPLICIT NONE
65
66
67
68  REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep
69  INTEGER,INTENT(IN) :: nday
70  INTEGER,INTENT(IN) :: ngrid,nlayer,nq
71  REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid)
72  integer,intent(in) :: day_ini
73  INTEGER ig,ierr
74 
75  EXTERNAL iniorbit,orbite
76  EXTERNAL SSUM
77  REAL SSUM
78 
79  ! Initialize flags lunout, prt_level, debug (in print_control_mod)
80  CALL init_print_control
81
82  ! initialize constants in comcstfi_mod
83  rad=prad
84  cpp=pcpp
85  g=pg
86  r=pr
87  rcp=r/cpp
88  pi=2.*asin(1.)
89  avocado = 6.02214179e23   ! added by RW
90
91  ! Initialize some "temporal and calendar" related variables
92  CALL init_time(day_ini,pdaysec,nday,ptimestep)
93
94  ! read in some parameters from "run.def" for physics,
95  ! or shared between dynamics and physics.
96  call getin_p("ecritphy",ecritphy) ! frequency of outputs in physics,
97                                    ! in dynamical steps
98  call getin_p("day_step",day_step) ! number of dynamical steps per day
99  call getin_p("iphysiq",iphysiq) ! call physics every iphysiq dyn step
100
101  ! do we read a startphy.nc file? (default: .true.)
102  call getin_p("startphy_file",startphy_file)
103 
104! --------------------------------------------------------------
105!  Reading the "callphys.def" file controlling some key options
106! --------------------------------------------------------------
107     
108  ! check that 'callphys.def' file is around
109  OPEN(99,file='callphys.def',status='old',form='formatted',iostat=ierr)
110  CLOSE(99)
111  IF(ierr.EQ.0) iscallphys=.true. !iscallphys initialised as false in callkeys_mod module
112     
113!!!      IF(ierr.EQ.0) THEN
114  IF(iscallphys) THEN
115     PRINT*
116     PRINT*
117     PRINT*,'--------------------------------------------'
118     PRINT*,' inifis: Parametres pour la physique (callphys.def)'
119     PRINT*,'--------------------------------------------'
120
121     write(*,*) "Directory where external input files are:"
122     ! default 'datadir' is set in "datadir_mod"
123     call getin_p("datadir",datadir) ! default path
124     write(*,*) " datadir = ",trim(datadir)
125
126     write(*,*) "Run with or without tracer transport ?"
127     tracer=.false. ! default value
128     call getin_p("tracer",tracer)
129     write(*,*) " tracer = ",tracer
130
131     write(*,*) "Run with or without atm mass update ", &
132            " due to tracer evaporation/condensation?"
133     mass_redistrib=.false. ! default value
134     call getin_p("mass_redistrib",mass_redistrib)
135     write(*,*) " mass_redistrib = ",mass_redistrib
136
137     write(*,*) "Diurnal cycle ?"
138     write(*,*) "(if diurnal=false, diurnal averaged solar heating)"
139     diurnal=.true. ! default value
140     call getin_p("diurnal",diurnal)
141     write(*,*) " diurnal = ",diurnal
142
143     write(*,*) "Seasonal cycle ?"
144     write(*,*) "(if season=false, Ls stays constant, to value ", &
145         "set in 'start'"
146     season=.true. ! default value
147     call getin_p("season",season)
148     write(*,*) " season = ",season
149
150     write(*,*) "Tidally resonant rotation ?"
151     tlocked=.false. ! default value
152     call getin_p("tlocked",tlocked)
153     write(*,*) "tlocked = ",tlocked
154
155     write(*,*) "Saturn ring shadowing ?"
156     rings_shadow = .false.
157     call getin_p("rings_shadow", rings_shadow)
158     write(*,*) "rings_shadow = ", rings_shadow
159         
160     write(*,*) "Compute latitude-dependent gravity field?"
161     oblate = .false.
162     call getin_p("oblate", oblate)
163     write(*,*) "oblate = ", oblate
164
165     write(*,*) "Flattening of the planet (a-b)/a "
166     flatten = 0.0
167     call getin_p("flatten", flatten)
168     write(*,*) "flatten = ", flatten
169         
170
171     write(*,*) "Needed if oblate=.true.: J2"
172     J2 = 0.0
173     call getin_p("J2", J2)
174     write(*,*) "J2 = ", J2
175         
176     write(*,*) "Needed if oblate=.true.: Planet mass (*1e24 kg)"
177     MassPlanet = 0.0
178     call getin_p("MassPlanet", MassPlanet)
179     write(*,*) "MassPlanet = ", MassPlanet         
180
181     write(*,*) "Needed if oblate=.true.: Planet mean radius (m)"
182     Rmean = 0.0
183     call getin_p("Rmean", Rmean)
184     write(*,*) "Rmean = ", Rmean
185         
186! Test of incompatibility:
187! if tlocked, then diurnal should be false
188     if (tlocked.and.diurnal) then
189       print*,'If diurnal=true, we should turn off tlocked.'
190       stop
191     endif
192
193     write(*,*) "Tidal resonance ratio ?"
194     nres=0          ! default value
195     call getin_p("nres",nres)
196     write(*,*) "nres = ",nres
197
198     write(*,*) "Write some extra output to the screen ?"
199     lwrite=.false. ! default value
200     call getin_p("lwrite",lwrite)
201     write(*,*) " lwrite = ",lwrite
202
203     write(*,*) "Save statistics in file stats.nc ?"
204     callstats=.true. ! default value
205     call getin_p("callstats",callstats)
206     write(*,*) " callstats = ",callstats
207
208     write(*,*) "Test energy conservation of model physics ?"
209     enertest=.false. ! default value
210     call getin_p("enertest",enertest)
211     write(*,*) " enertest = ",enertest
212
213     write(*,*) "Check to see if cpp values used match gases.def ?"
214     check_cpp_match=.true. ! default value
215     call getin_p("check_cpp_match",check_cpp_match)
216     write(*,*) " check_cpp_match = ",check_cpp_match
217
218     write(*,*) "call radiative transfer ?"
219     callrad=.true. ! default value
220     call getin_p("callrad",callrad)
221     write(*,*) " callrad = ",callrad
222
223     write(*,*) "call correlated-k radiative transfer ?"
224     corrk=.true. ! default value
225     call getin_p("corrk",corrk)
226     write(*,*) " corrk = ",corrk
227
228     write(*,*) "prohibit calculations outside corrk T grid?"
229     strictboundcorrk=.true. ! default value
230     call getin_p("strictboundcorrk",strictboundcorrk)
231     write(*,*) "strictboundcorrk = ",strictboundcorrk
232
233     write(*,*) "call gaseous absorption in the visible bands?", &
234                    "(matters only if callrad=T)"
235     callgasvis=.false. ! default value
236     call getin_p("callgasvis",callgasvis)
237     write(*,*) " callgasvis = ",callgasvis
238       
239     write(*,*) "call continuum opacities in radiative transfer ?", &
240                    "(matters only if callrad=T)"
241     continuum=.true. ! default value
242     call getin_p("continuum",continuum)
243     write(*,*) " continuum = ",continuum
244
245     write(*,*) "use analytic function for H2O continuum ?"
246     H2Ocont_simple=.false. ! default value
247     call getin_p("H2Ocont_simple",H2Ocont_simple)
248     write(*,*) " H2Ocont_simple = ",H2Ocont_simple
249 
250     write(*,*) "call turbulent vertical diffusion ?"
251     calldifv=.true. ! default value
252     call getin_p("calldifv",calldifv)
253     write(*,*) " calldifv = ",calldifv
254
255     write(*,*) "use turbdiff instead of vdifc ?"
256     UseTurbDiff=.true. ! default value
257     call getin_p("UseTurbDiff",UseTurbDiff)
258     write(*,*) " UseTurbDiff = ",UseTurbDiff
259
260     write(*,*) "call convective adjustment ?"
261     calladj=.true. ! default value
262     call getin_p("calladj",calladj)
263     write(*,*) " calladj = ",calladj
264
265     write(*,*) "call CO2 condensation ?"
266     co2cond=.false. ! default value
267     call getin_p("co2cond",co2cond)
268     write(*,*) " co2cond = ",co2cond
269! Test of incompatibility
270     if (co2cond.and.(.not.tracer)) then
271        print*,'We need a CO2 ice tracer to condense CO2'
272        call abort
273     endif
274 
275     write(*,*) "CO2 supersaturation level ?"
276     co2supsat=1.0 ! default value
277     call getin_p("co2supsat",co2supsat)
278     write(*,*) " co2supsat = ",co2supsat
279
280     write(*,*) "Radiative timescale for Newtonian cooling ?"
281     tau_relax=30. ! default value
282     call getin_p("tau_relax",tau_relax)
283     write(*,*) " tau_relax = ",tau_relax
284     tau_relax=tau_relax*24*3600 ! convert Earth days --> seconds
285
286     write(*,*)"call thermal conduction in the soil ?"
287     callsoil=.true. ! default value
288     call getin_p("callsoil",callsoil)
289     write(*,*) " callsoil = ",callsoil
290         
291     write(*,*)"Rad transfer is computed every iradia", &
292                   " physical timestep"
293     iradia=1 ! default value
294     call getin_p("iradia",iradia)
295     write(*,*)" iradia = ",iradia
296       
297     write(*,*)"Rayleigh scattering ?"
298     rayleigh=.false.
299     call getin_p("rayleigh",rayleigh)
300     write(*,*)" rayleigh = ",rayleigh
301
302     write(*,*) "Use blackbody for stellar spectrum ?"
303     stelbbody=.false. ! default value
304     call getin_p("stelbbody",stelbbody)
305     write(*,*) " stelbbody = ",stelbbody
306
307     write(*,*) "Stellar blackbody temperature ?"
308     stelTbb=5800.0 ! default value
309     call getin_p("stelTbb",stelTbb)
310     write(*,*) " stelTbb = ",stelTbb
311
312     write(*,*)"Output mean OLR in 1D?"
313     meanOLR=.false.
314     call getin_p("meanOLR",meanOLR)
315     write(*,*)" meanOLR = ",meanOLR
316
317     write(*,*)"Output spectral OLR in 3D?"
318     specOLR=.false.
319     call getin_p("specOLR",specOLR)
320     write(*,*)" specOLR = ",specOLR
321
322     write(*,*)"Operate in kastprof mode?"
323     kastprof=.false.
324     call getin_p("kastprof",kastprof)
325     write(*,*)" kastprof = ",kastprof
326
327     write(*,*)"Uniform absorption in radiative transfer?"
328     graybody=.false.
329     call getin_p("graybody",graybody)
330     write(*,*)" graybody = ",graybody
331
332! Soil model
333     write(*,*)"Number of sub-surface layers for soil scheme?"
334     ! default value of nsoilmx set in comsoil_h
335     call getin_p("nsoilmx",nsoilmx)
336     write(*,*)" nsoilmx=",nsoilmx
337     
338     write(*,*)"Thickness of topmost soil layer (m)?"
339     ! default value of lay1_soil set in comsoil_h
340     call getin_p("lay1_soil",lay1_soil)
341     write(*,*)" lay1_soil = ",lay1_soil
342     
343     write(*,*)"Coefficient for soil layer thickness distribution?"
344     ! default value of alpha_soil set in comsoil_h
345     call getin_p("alpha_soil",alpha_soil)
346     write(*,*)" alpha_soil = ",alpha_soil
347
348! Slab Ocean
349     write(*,*) "Use slab-ocean ?"
350     ok_slab_ocean=.false.         ! default value
351     call getin_p("ok_slab_ocean",ok_slab_ocean)
352     write(*,*) "ok_slab_ocean = ",ok_slab_ocean
353     ! Sanity check: for now slab oncean only works in serial mode
354     if (ok_slab_ocean.and.is_parallel) then
355       write(*,*) " Error: slab ocean should only be used in serial mode!"
356       call abort
357     endif
358
359     write(*,*) "Use slab-sea-ice ?"
360     ok_slab_sic=.true.         ! default value
361     call getin_p("ok_slab_sic",ok_slab_sic)
362     write(*,*) "ok_slab_sic = ",ok_slab_sic
363
364     write(*,*) "Use heat transport for the ocean ?"
365     ok_slab_heat_transp=.true.   ! default value
366     call getin_p("ok_slab_heat_transp",ok_slab_heat_transp)
367     write(*,*) "ok_slab_heat_transp = ",ok_slab_heat_transp
368
369
370
371! Test of incompatibility:
372! if kastprof used, we must be in 1D
373     if (kastprof.and.(ngrid.gt.1)) then
374       print*,'kastprof can only be used in 1D!'
375       call abort
376     endif
377
378     write(*,*)"Stratospheric temperature for kastprof mode?"
379     Tstrat=167.0
380     call getin_p("Tstrat",Tstrat)
381     write(*,*)" Tstrat = ",Tstrat
382
383     write(*,*)"Remove lower boundary?"
384     nosurf=.false.
385     call getin_p("nosurf",nosurf)
386     write(*,*)" nosurf = ",nosurf
387
388! Tests of incompatibility:
389     if (nosurf.and.callsoil) then
390       print*,'nosurf not compatible with soil scheme!'
391       print*,'... got to make a choice!'
392       call abort
393     endif
394
395     write(*,*)"Add an internal heat flux?", &
396                   "... matters only if callsoil=F"
397     intheat=0.
398     call getin_p("intheat",intheat)
399     write(*,*)" intheat = ",intheat
400
401     write(*,*)"Use Newtonian cooling for radiative transfer?"
402     newtonian=.false.
403     call getin_p("newtonian",newtonian)
404     write(*,*)" newtonian = ",newtonian
405
406! Tests of incompatibility:
407     if (newtonian.and.corrk) then
408       print*,'newtonian not compatible with correlated-k!'
409       call abort
410     endif
411     if (newtonian.and.calladj) then
412       print*,'newtonian not compatible with adjustment!'
413       call abort
414     endif
415     if (newtonian.and.calldifv) then
416       print*,'newtonian not compatible with a boundary layer!'
417       call abort
418     endif
419
420     write(*,*)"Test physics timescale in 1D?"
421     testradtimes=.false.
422     call getin_p("testradtimes",testradtimes)
423     write(*,*)" testradtimes = ",testradtimes
424
425! Test of incompatibility:
426! if testradtimes used, we must be in 1D
427     if (testradtimes.and.(ngrid.gt.1)) then
428       print*,'testradtimes can only be used in 1D!'
429       call abort
430     endif
431
432     write(*,*)"Default planetary temperature?"
433     tplanet=215.0
434     call getin_p("tplanet",tplanet)
435     write(*,*)" tplanet = ",tplanet
436
437     write(*,*)"Which star?"
438     startype=1 ! default value = Sol
439     call getin_p("startype",startype)
440     write(*,*)" startype = ",startype
441
442     write(*,*)"Value of stellar flux at 1 AU?"
443     Fat1AU=1356.0 ! default value = Sol today
444     call getin_p("Fat1AU",Fat1AU)
445     write(*,*)" Fat1AU = ",Fat1AU
446
447
448! TRACERS:
449
450     write(*,*)"Varying H2O cloud fraction?"
451     CLFvarying=.false.     ! default value
452     call getin_p("CLFvarying",CLFvarying)
453     write(*,*)" CLFvarying = ",CLFvarying
454
455     write(*,*)"Value of fixed H2O cloud fraction?"
456     CLFfixval=1.0                ! default value
457     call getin_p("CLFfixval",CLFfixval)
458     write(*,*)" CLFfixval = ",CLFfixval
459
460     write(*,*)"fixed radii for Cloud particles?"
461     radfixed=.false. ! default value
462     call getin_p("radfixed",radfixed)
463     write(*,*)" radfixed = ",radfixed
464
465     if(kastprof)then
466        radfixed=.true.
467     endif 
468
469     write(*,*)"Number mixing ratio of CO2 ice particles:"
470     Nmix_co2=1.e6 ! default value
471     call getin_p("Nmix_co2",Nmix_co2)
472     write(*,*)" Nmix_co2 = ",Nmix_co2
473
474!         write(*,*)"Number of radiatively active aerosols:"
475!         naerkind=0. ! default value
476!         call getin_p("naerkind",naerkind)
477!         write(*,*)" naerkind = ",naerkind
478
479     write(*,*)"Opacity of dust (if used):"
480     dusttau=0. ! default value
481     call getin_p("dusttau",dusttau)
482     write(*,*)" dusttau = ",dusttau
483
484     write(*,*)"Radiatively active CO2 aerosols?"
485     aeroco2=.false.     ! default value
486     call getin_p("aeroco2",aeroco2)
487     write(*,*)" aeroco2 = ",aeroco2
488
489     write(*,*)"Fixed CO2 aerosol distribution?"
490     aerofixco2=.false.     ! default value
491     call getin_p("aerofixco2",aerofixco2)
492     write(*,*)" aerofixco2 = ",aerofixco2
493
494     write(*,*)"Radiatively active water ice?"
495     aeroh2o=.false.     ! default value
496     call getin_p("aeroh2o",aeroh2o)
497     write(*,*)" aeroh2o = ",aeroh2o
498
499     write(*,*)"Fixed H2O aerosol distribution?"
500     aerofixh2o=.false.     ! default value
501     call getin_p("aerofixh2o",aerofixh2o)
502     write(*,*)" aerofixh2o = ",aerofixh2o
503
504     write(*,*)"Radiatively active sulfuric acid aerosols?"
505     aeroh2so4=.false.     ! default value
506     call getin_p("aeroh2so4",aeroh2so4)
507     write(*,*)" aeroh2so4 = ",aeroh2so4
508         
509!=================================
510
511     write(*,*)"Radiatively active two-layer aerosols?"
512     aeroback2lay=.false.     ! default value
513     call getin_p("aeroback2lay",aeroback2lay)
514     write(*,*)" aeroback2lay = ",aeroback2lay
515
516     write(*,*)"Radiatively active ammonia cloud?"
517     aeronh3=.false.     ! default value
518     call getin_p("aeronh3",aeronh3)
519     write(*,*)" aeronh3 = ",aeronh3
520
521     write(*,*)"Radiatively active auroral aerosols?"
522     aeroaurora=.false.     ! default value
523     call getin_p("aeroaurora",aeroaurora)
524     write(*,*)" aeroaurora = ",aeroaurora
525
526     write(*,*)"TWOLAY AEROSOL: total optical depth ", &
527                    "in the tropospheric layer (visible)"
528     obs_tau_col_tropo=8.D0
529     call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
530     write(*,*)" obs_tau_col_tropo = ",obs_tau_col_tropo
531
532     write(*,*)"TWOLAY AEROSOL: total optical depth ", &
533                    "in the stratospheric layer (visible)"
534     obs_tau_col_strato=0.08D0
535     call getin_p("obs_tau_col_strato",obs_tau_col_strato)
536     write(*,*)" obs_tau_col_strato = ",obs_tau_col_strato
537
538     write(*,*)"TWOLAY AEROSOL: pres_bottom_tropo? in pa"
539     pres_bottom_tropo=66000.0
540     call getin_p("pres_bottom_tropo",pres_bottom_tropo)
541     write(*,*)" pres_bottom_tropo = ",pres_bottom_tropo
542
543     write(*,*)"TWOLAY AEROSOL: pres_top_tropo? in pa"
544     pres_top_tropo=18000.0
545     call getin_p("pres_top_tropo",pres_top_tropo)
546     write(*,*)" pres_top_tropo = ",pres_top_tropo
547
548     write(*,*)"TWOLAY AEROSOL: pres_bottom_strato? in pa"
549     pres_bottom_strato=2000.0
550     call getin_p("pres_bottom_strato",pres_bottom_strato)
551     write(*,*)" pres_bottom_strato = ",pres_bottom_strato
552
553     write(*,*)"TWOLAY AEROSOL: pres_top_strato? in pa"
554     pres_top_strato=100.0
555     call getin_p("pres_top_strato",pres_top_strato)
556     write(*,*)" pres_top_strato = ",pres_top_strato
557
558     write(*,*)"TWOLAY AEROSOL: particle size in the ", &
559                    "tropospheric layer, in meters"
560     size_tropo=2.e-6
561     call getin_p("size_tropo",size_tropo)
562     write(*,*)" size_tropo = ",size_tropo
563
564     write(*,*)"TWOLAY AEROSOL: particle size in the ", &
565                    "stratospheric layer, in meters"
566     size_strato=1.e-7
567     call getin_p("size_strato",size_strato)
568     write(*,*)" size_strato = ",size_strato
569
570     write(*,*)"NH3 (thin) cloud: total optical depth"
571     tau_nh3_cloud=7.D0
572     call getin_p("tau_nh3_cloud",tau_nh3_cloud)
573     write(*,*)" tau_nh3_cloud = ",tau_nh3_cloud
574
575     write(*,*)"NH3 (thin) cloud pressure level"
576     pres_nh3_cloud=7.D0
577     call getin_p("pres_nh3_cloud",pres_nh3_cloud)
578     write(*,*)" pres_nh3_cloud = ",pres_nh3_cloud
579
580     write(*,*)"NH3 (thin) cloud: particle sizes"
581     size_nh3_cloud=3.e-6
582     call getin_p("size_nh3_cloud",size_nh3_cloud)
583     write(*,*)" size_nh3_cloud = ",size_nh3_cloud
584
585!=================================
586
587     write(*,*)"Cloud pressure level (with kastprof only):"
588     cloudlvl=0. ! default value
589     call getin_p("cloudlvl",cloudlvl)
590     write(*,*)" cloudlvl = ",cloudlvl
591
592     write(*,*)"Is the variable gas species radiatively active?"
593     Tstrat=167.0
594     varactive=.false.
595     call getin_p("varactive",varactive)
596     write(*,*)" varactive = ",varactive
597
598     write(*,*)"Is the variable gas species distribution set?"
599     varfixed=.false.
600     call getin_p("varfixed",varfixed)
601     write(*,*)" varfixed = ",varfixed
602
603     write(*,*)"What is the saturation % of the variable species?"
604     satval=0.8
605     call getin_p("satval",satval)
606     write(*,*)" satval = ",satval
607
608
609! Test of incompatibility:
610! if varactive, then varfixed should be false
611     if (varactive.and.varfixed) then
612       print*,'if varactive, varfixed must be OFF!'
613       stop
614     endif
615
616     write(*,*) "Gravitationnal sedimentation ?"
617     sedimentation=.false. ! default value
618     call getin_p("sedimentation",sedimentation)
619     write(*,*) " sedimentation = ",sedimentation
620
621     write(*,*) "Compute water cycle ?"
622     water=.false. ! default value
623     call getin_p("water",water)
624     write(*,*) " water = ",water
625         
626! Test of incompatibility:
627! if water is true, there should be at least a tracer
628     if (water.and.(.not.tracer)) then
629        print*,'if water is ON, tracer must be ON too!'
630        stop
631     endif
632
633     write(*,*) "Include water condensation ?"
634     watercond=.false. ! default value
635     call getin_p("watercond",watercond)
636     write(*,*) " watercond = ",watercond
637
638! Test of incompatibility:
639! if watercond is used, then water should be used too
640     if (watercond.and.(.not.water)) then
641        print*,'if watercond is used, water should be used too'
642        stop
643     endif
644
645     write(*,*) "Include water precipitation ?"
646     waterrain=.false. ! default value
647     call getin_p("waterrain",waterrain)
648     write(*,*) " waterrain = ",waterrain
649
650     write(*,*) "Include surface hydrology ?"
651     hydrology=.false. ! default value
652     call getin_p("hydrology",hydrology)
653     write(*,*) " hydrology = ",hydrology
654
655     write(*,*) "Evolve surface water sources ?"
656     sourceevol=.false. ! default value
657     call getin_p("sourceevol",sourceevol)
658     write(*,*) " sourceevol = ",sourceevol
659
660     write(*,*) "Ice evolution timestep ?"
661     icetstep=100.0 ! default value
662     call getin_p("icetstep",icetstep)
663     write(*,*) " icetstep = ",icetstep
664         
665     write(*,*) "Spectral Dependant albedo ?"
666     albedo_spectral_mode=.false. ! default value
667     call getin_p("albedo_spectral_mode",albedo_spectral_mode)
668     write(*,*) " albedo_spectral_mode = ",albedo_spectral_mode
669
670     write(*,*) "Snow albedo ?"
671     write(*,*) "If albedo_spectral_mode=.true., then this "
672     write(*,*) "corresponds to the 0.5 microns snow albedo."
673     albedosnow=0.5         ! default value
674     call getin_p("albedosnow",albedosnow)
675     write(*,*) " albedosnow = ",albedosnow
676         
677     write(*,*) "CO2 ice albedo ?"
678     albedoco2ice=0.5       ! default value
679     call getin_p("albedoco2ice",albedoco2ice)
680     write(*,*) " albedoco2ice = ",albedoco2ice
681
682     write(*,*) "Maximum ice thickness ?"
683     maxicethick=2.0         ! default value
684     call getin_p("maxicethick",maxicethick)
685     write(*,*) " maxicethick = ",maxicethick
686
687     write(*,*) "Freezing point of seawater ?"
688     Tsaldiff=-1.8          ! default value
689     call getin_p("Tsaldiff",Tsaldiff)
690     write(*,*) " Tsaldiff = ",Tsaldiff
691
692     write(*,*) "Does user want to force cpp and mugaz?"
693     force_cpp=.false. ! default value
694     call getin_p("force_cpp",force_cpp)
695     write(*,*) " force_cpp = ",force_cpp
696
697     if (force_cpp) then
698       mugaz = -99999.
699       PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
700       call getin_p("mugaz",mugaz)
701       IF (mugaz.eq.-99999.) THEN
702           PRINT *, "mugaz must be set if force_cpp = T"
703           STOP
704       ELSE
705           write(*,*) "mugaz=",mugaz
706       ENDIF
707       cpp = -99999.
708       PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
709       call getin_p("cpp",cpp)
710       IF (cpp.eq.-99999.) THEN
711           PRINT *, "cpp must be set if force_cpp = T"
712           STOP
713       ELSE
714           write(*,*) "cpp=",cpp
715       ENDIF
716!         else
717!           mugaz=8.314*1000./pr
718     endif ! of if (force_cpp)
719     call su_gases
720     call calc_cpp_mugaz
721
722     PRINT*,'--------------------------------------------'
723     PRINT*
724     PRINT*
725  ELSE
726     write(*,*)
727     write(*,*) 'Cannot read file callphys.def. Is it here ?'
728     stop
729  ENDIF ! of IF(iscallphys)
730
731  PRINT*
732  PRINT*,'inifis: daysec',daysec
733  PRINT*
734  PRINT*,'inifis: The radiative transfer is computed:'
735  PRINT*,'           each ',iradia,' physical time-step'
736  PRINT*,'        or each ',iradia*dtphys,' seconds'
737  PRINT*
738
739
740!-----------------------------------------------------------------------
741!     Some more initialization:
742!     ------------------------
743
744  ! Initializations for comgeomfi_h
745  totarea=SSUM(ngrid,parea,1)
746  call planetwide_sumval(parea,totarea_planet)
747
748  !! those are defined in comdiurn_h.F90
749  IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
750  IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
751  IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
752  IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
753
754  DO ig=1,ngrid
755     sinlat(ig)=sin(plat(ig))
756     coslat(ig)=cos(plat(ig))
757     sinlon(ig)=sin(plon(ig))
758     coslon(ig)=cos(plon(ig))
759  ENDDO
760
761  ! initialize variables in radinc_h
762  call ini_radinc_h(nlayer)
763 
764  ! allocate "radcommon_h" arrays
765  call ini_radcommon_h()
766
767  ! allocate "comsoil_h" arrays
768  call ini_comsoil_h(ngrid)
769     
770  END SUBROUTINE inifis
771
772END MODULE inifis_mod
Note: See TracBrowser for help on using the repository browser.