source: trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90 @ 2050

Last change on this file since 2050 was 2050, checked in by jvatant, 6 years ago

Major radiative transfer contribution : Add the 'corrk_recombin' option that allows to use
correlated-k for single species instead of pre-mix and enables more flexiblity for variable species.
-> Algorithm inspired from Lacis and Oinas 1991 and Amundsen et al 2016

+ Added 'recombin_corrk.F90' - Important improvements in 'sugas_corrk.90' and 'callcork.F90'

  • Must have the desired variable species as tracers -> TBD : Enable to force composition even if no tracers
  • To have decent CPU time recombining is not done on all gridpoints and wavelenghts but we calculate a gasi/v_recomb variable on the reference corrk-k T,P grid (only for T,P values who match the atmospheric conditions ) which is then processed as a standard pre-mix in optci/v routines, but updated every time tracers on the ref P grid have varied > 1%.


READ CAREFULY :

  • In case of 'corrk_recombin', the variable L_NREFVAR doesn't have the same meaning as before and doesn't stand for the different mixing ratios but the different species.
  • Input corr-k should be found in corrkdir within 'corrk_gcm_IR/VI_XXX.dat' and can contain a 'fixed' specie ( compulsory if you include self-broadening ) that MUST have been created with correct mixing ratios, or a variable specie for which mixing ratio MUST have been set to 1 ( no self-broadening then, assume it's a trace sepecie ) -> You can't neither have CIA of variable species included upstream in the corr-k
File size: 19.1 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
13  use datafile_mod
14  use comdiurn_h, only: sinlat, coslat, sinlon, coslon
15  use comgeomfi_h, only: totarea, totarea_planet
16  use comsoil_h, only: ini_comsoil_h, nsoilmx, lay1_soil, alpha_soil
17  use time_phylmdz_mod, only: ecritphy,day_step,iphysiq, &
18                              init_time, daysec, dtphys
19  use comcstfi_mod, only: rad, cpp, g, r, rcp, &
20                          mugaz, pi, avocado, kbol
21  use planete_mod, only: nres
22  use planetwide_mod, only: planetwide_sumval
23  use callkeys_mod
24  use mod_phys_lmdz_para, only : is_parallel
25
26!=======================================================================
27!
28!   purpose:
29!   -------
30!
31!   Initialisation for the physical parametrisations of the LMD
32!   Generic Model.
33!
34!   author: Frederic Hourdin 15 / 10 /93
35!   -------
36!   modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
37!             Ehouarn Millour (oct. 2008) tracers are now identified
38!              by their names and may not be contiguously
39!              stored in the q(:,:,:,:) array
40!             E.M. (june 2009) use getin routine to load parameters
41!
42!
43!   arguments:
44!   ----------
45!
46!   input:
47!   ------
48!
49!    ngrid                 Size of the horizontal grid.
50!                          All internal loops are performed on that grid.
51!    nlayer                Number of vertical layers.
52!    pdayref               Day of reference for the simulation
53!    pday                  Number of days counted from the North. Spring
54!                          equinoxe.
55!
56!=======================================================================
57!
58!-----------------------------------------------------------------------
59!   declarations:
60!   -------------
61  use ioipsl_getin_p_mod, only: getin_p
62  IMPLICIT NONE
63
64
65
66  REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep
67  INTEGER,INTENT(IN) :: nday
68  INTEGER,INTENT(IN) :: ngrid,nlayer,nq
69  REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid)
70  integer,intent(in) :: day_ini
71  INTEGER ig,ierr
72 
73  EXTERNAL iniorbit,orbite
74  EXTERNAL SSUM
75  REAL SSUM
76 
77  ! Initialize flags lunout, prt_level, debug (in print_control_mod)
78  CALL init_print_control
79
80  ! initialize constants in comcstfi_mod
81  rad=prad
82  cpp=pcpp
83  g=pg
84  r=pr
85  rcp=r/cpp
86  pi=2.*asin(1.)
87  avocado = 6.02214179e23   ! added by RW
88  kbol = 1.38064852e-23  ! added by JVO for Titan chem
89
90  ! Initialize some "temporal and calendar" related variables
91  CALL init_time(day_ini,pdaysec,nday,ptimestep)
92
93  ! read in some parameters from "run.def" for physics,
94  ! or shared between dynamics and physics.
95  call getin_p("ecritphy",ecritphy) ! frequency of outputs in physics,
96                                    ! in dynamical steps
97  call getin_p("day_step",day_step) ! number of dynamical steps per day
98  call getin_p("iphysiq",iphysiq) ! call physics every iphysiq dyn step
99
100  ! do we read a startphy.nc file? (default: .true.)
101  call getin_p("startphy_file",startphy_file)
102 
103! --------------------------------------------------------------
104!  Reading the "callphys.def" file controlling some key options
105! --------------------------------------------------------------
106     
107  ! check that 'callphys.def' file is around
108  OPEN(99,file='callphys.def',status='old',form='formatted',iostat=ierr)
109  CLOSE(99)
110  IF(ierr.EQ.0) iscallphys=.true. !iscallphys initialised as false in callkeys_mod module
111     
112!!!      IF(ierr.EQ.0) THEN
113  IF(iscallphys) THEN
114     PRINT*
115     PRINT*
116     PRINT*,'--------------------------------------------'
117     PRINT*,' inifis: Parametres pour la physique (callphys.def)'
118     PRINT*,'--------------------------------------------'
119
120     write(*,*) "Directory where external input files are:"
121     ! default 'datadir' is set in "datadir_mod"
122     call getin_p("datadir",datadir) ! default path
123     write(*,*) " datadir = ",trim(datadir)
124
125     write(*,*) "Run with or without tracer transport ?"
126     tracer=.false. ! default value
127     call getin_p("tracer",tracer)
128     write(*,*) " tracer = ",tracer
129
130     write(*,*) "Run with or without atm mass update ", &
131            " due to tracer evaporation/condensation?"
132     mass_redistrib=.false. ! default value
133     call getin_p("mass_redistrib",mass_redistrib)
134     write(*,*) " mass_redistrib = ",mass_redistrib
135
136     write(*,*) "Diurnal cycle ?"
137     write(*,*) "(if diurnal=false, diurnal averaged solar heating)"
138     diurnal=.true. ! default value
139     call getin_p("diurnal",diurnal)
140     write(*,*) " diurnal = ",diurnal
141
142     write(*,*) "Seasonal cycle ?"
143     write(*,*) "(if season=false, Ls stays constant, to value ", &
144         "set in 'start'"
145     season=.true. ! default value
146     call getin_p("season",season)
147     write(*,*) " season = ",season
148
149     write(*,*) "Tidally resonant rotation ?"
150     tlocked=.false. ! default value
151     call getin_p("tlocked",tlocked)
152     write(*,*) "tlocked = ",tlocked
153
154     write(*,*) "Saturn ring shadowing ?"
155     rings_shadow = .false.
156     call getin_p("rings_shadow", rings_shadow)
157     write(*,*) "rings_shadow = ", rings_shadow
158         
159     write(*,*) "Compute latitude-dependent gravity field?"
160     oblate = .false.
161     call getin_p("oblate", oblate)
162     write(*,*) "oblate = ", oblate
163
164     write(*,*) "Flattening of the planet (a-b)/a "
165     flatten = 0.0
166     call getin_p("flatten", flatten)
167     write(*,*) "flatten = ", flatten         
168
169     write(*,*) "Needed if oblate=.true.: J2"
170     J2 = 0.0
171     call getin_p("J2", J2)
172     write(*,*) "J2 = ", J2
173         
174     write(*,*) "Needed if oblate=.true.: Planet mass (*1e24 kg)"
175     MassPlanet = 0.0
176     call getin_p("MassPlanet", MassPlanet)
177     write(*,*) "MassPlanet = ", MassPlanet         
178
179     write(*,*) "Needed if oblate=.true.: Planet mean radius (m)"
180     Rmean = 0.0
181     call getin_p("Rmean", Rmean)
182     write(*,*) "Rmean = ", Rmean
183     
184     write(*,*) "Compute effective altitude-dependent gravity field?"
185     eff_gz = .false.
186     call getin_p("eff_gz", eff_gz)
187     write(*,*) "eff_gz = ", eff_gz
188         
189! Test of incompatibility:
190! if tlocked, then diurnal should be false
191     if (tlocked.and.diurnal) then
192       print*,'If diurnal=true, we should turn off tlocked.'
193       stop
194     endif
195
196     write(*,*) "Tidal resonance ratio ?"
197     nres=0          ! default value
198     call getin_p("nres",nres)
199     write(*,*) "nres = ",nres
200
201     write(*,*) "Write some extra output to the screen ?"
202     lwrite=.false. ! default value
203     call getin_p("lwrite",lwrite)
204     write(*,*) " lwrite = ",lwrite
205
206     write(*,*) "Save statistics in file stats.nc ?"
207     callstats=.true. ! default value
208     call getin_p("callstats",callstats)
209     write(*,*) " callstats = ",callstats
210
211     write(*,*) "Test energy conservation of model physics ?"
212     enertest=.false. ! default value
213     call getin_p("enertest",enertest)
214     write(*,*) " enertest = ",enertest
215
216     write(*,*) "Check to see if cpp values used match gases.def ?"
217     check_cpp_match=.true. ! default value
218     call getin_p("check_cpp_match",check_cpp_match)
219     write(*,*) " check_cpp_match = ",check_cpp_match
220
221     write(*,*) "call radiative transfer ?"
222     callrad=.true. ! default value
223     call getin_p("callrad",callrad)
224     write(*,*) " callrad = ",callrad
225
226     write(*,*) "call correlated-k radiative transfer ?"
227     corrk=.true. ! default value
228     call getin_p("corrk",corrk)
229     write(*,*) " corrk = ",corrk
230     
231     if (corrk) then
232       ! default path is set in datadir
233       write(*,*) "callcorrk: Correlated-k data base folder:",trim(datadir)
234       call getin_p("corrkdir",corrkdir)
235       write(*,*) " corrkdir = ",corrkdir
236       
237       write(*,*) "use correlated-k recombination instead of pre-mixed values ?"
238       corrk_recombin=.true. ! default value
239       call getin_p("corrk_recombin",corrk_recombin)
240       write(*,*) " corrk_recombin = ",corrk_recombin
241     endif
242     
243     if (corrk .and. ngrid.eq.1) then
244       write(*,*) "simulate global averaged conditions ?"
245       global1d = .false. ! default value
246       call getin_p("global1d",global1d)
247       write(*,*) " global1d = ",global1d
248       
249       ! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
250       if (global1d.and.diurnal) then
251          write(*,*) "if global1d is true, diurnal must be set to false"
252          stop
253       endif
254
255       if (global1d) then
256          write(*,*) "Solar Zenith angle (deg.) ?"
257          write(*,*) "(assumed for averaged solar flux S/4)"
258          szangle=60.0  ! default value
259          call getin_p("szangle",szangle)
260          write(*,*) " szangle = ",szangle
261       endif
262     endif
263
264     write(*,*) "prohibit calculations outside corrk T grid?"
265     strictboundcorrk=.true. ! default value
266     call getin_p("strictboundcorrk",strictboundcorrk)
267     write(*,*) "strictboundcorrk = ",strictboundcorrk
268
269     write(*,*) "call gaseous absorption in the visible bands?", &
270                    "(matters only if callrad=T)"
271     callgasvis=.false. ! default value
272     call getin_p("callgasvis",callgasvis)
273     write(*,*) " callgasvis = ",callgasvis
274       
275     write(*,*) "call continuum opacities in radiative transfer ?", &
276                    "(matters only if callrad=T)"
277     continuum=.true. ! default value
278     call getin_p("continuum",continuum)
279     write(*,*) " continuum = ",continuum
280 
281     write(*,*) "call turbulent vertical diffusion ?"
282     calldifv=.true. ! default value
283     call getin_p("calldifv",calldifv)
284     write(*,*) " calldifv = ",calldifv
285
286     write(*,*) "use turbdiff instead of vdifc ?"
287     UseTurbDiff=.true. ! default value
288     call getin_p("UseTurbDiff",UseTurbDiff)
289     write(*,*) " UseTurbDiff = ",UseTurbDiff
290
291     write(*,*) "call convective adjustment ?"
292     calladj=.true. ! default value
293     call getin_p("calladj",calladj)
294     write(*,*) " calladj = ",calladj
295
296     write(*,*) "Radiative timescale for Newtonian cooling ?"
297     tau_relax=30. ! default value
298     call getin_p("tau_relax",tau_relax)
299     write(*,*) " tau_relax = ",tau_relax
300     tau_relax=tau_relax*24*3600 ! convert Earth days --> seconds
301
302     write(*,*)"call thermal conduction in the soil ?"
303     callsoil=.true. ! default value
304     call getin_p("callsoil",callsoil)
305     write(*,*) " callsoil = ",callsoil
306         
307     write(*,*)"Rad transfer is computed every iradia", &
308                   " physical timestep"
309     iradia=1 ! default value
310     call getin_p("iradia",iradia)
311     write(*,*)" iradia = ",iradia
312       
313     write(*,*)"Rayleigh scattering ?"
314     rayleigh=.false.
315     call getin_p("rayleigh",rayleigh)
316     write(*,*)" rayleigh = ",rayleigh
317
318     write(*,*) "Use blackbody for stellar spectrum ?"
319     stelbbody=.false. ! default value
320     call getin_p("stelbbody",stelbbody)
321     write(*,*) " stelbbody = ",stelbbody
322
323     write(*,*) "Stellar blackbody temperature ?"
324     stelTbb=5800.0 ! default value
325     call getin_p("stelTbb",stelTbb)
326     write(*,*) " stelTbb = ",stelTbb
327
328     write(*,*)"Output mean OLR in 1D?"
329     meanOLR=.false.
330     call getin_p("meanOLR",meanOLR)
331     write(*,*)" meanOLR = ",meanOLR
332
333     write(*,*)"Output spectral OLR in 3D?"
334     specOLR=.false.
335     call getin_p("specOLR",specOLR)
336     write(*,*)" specOLR = ",specOLR
337
338     write(*,*)"Uniform absorption in radiative transfer?"
339     graybody=.false.
340     call getin_p("graybody",graybody)
341     write(*,*)" graybody = ",graybody
342
343! Chemistry
344
345     write(*,*) "Run with or without chemistry?"
346     callchim=.false. ! default value
347     call getin_p("callchim",callchim)
348     write(*,*) " callchim = ",callchim
349
350     ! sanity check
351     if (callchim.and.(.not.tracer)) then
352       print*,"You are running chemistry without tracer"
353       print*,"Please start again with tracer =.true."
354       stop
355     endif
356     
357     ! sanity check warning
358     if (callchim.and.(.not.eff_gz)) then
359       print*,"WARNING : You run chemistry without effective altitude-dependent gravity field !!"
360       print*,"You will have wrong vertical photodissociations rates, that's crazy !!"
361       print*,"I let you continue but you should rather set eff_gz =.true. ..."
362       !stop
363     endif
364
365
366     write(*,*)"Chemistry is computed every ichim", &
367                   " physical timestep"
368     ichim=1 ! default value
369     call getin_p("ichim",ichim)
370     write(*,*)" ichim = ",ichim
371
372! Microphysics
373
374     write(*,*) "Use haze seasonal model from Karkoschka 2016 ?"
375     seashaze=.true. ! default value
376     call getin_p("seashaze",seashaze)
377     write(*,*)" seashaze = ",seashaze
378
379     write(*,*) "Run with or without microphysics?"
380     callmufi=.false. ! default value
381     call getin_p("callmufi",callmufi)
382     write(*,*)" callmufi = ",callmufi
383
384     ! sanity check
385     if (callmufi.and.(.not.tracer)) then
386       print*,"You are running microphysics without tracer"
387       print*,"Please start again with tracer =.true."
388       stop
389     endif
390
391     write(*,*) "Compute clouds?"
392     callclouds=.false. ! default value
393     call getin_p("callclouds",callclouds)
394     write(*,*)" callclouds = ",callclouds
395
396     ! sanity check
397     if (callclouds.and.(.not.callmufi)) then
398       print*,"You are trying to make clouds without microphysics !"
399       print*,"Please start again with callmufi =.true."
400       stop
401     endif
402
403     write(*,*) "Disable the coupling of microphysics within rad. transf. ?"
404     write(*,*) "If disabled we will assume a planetwide vert. profile of extinction ..."
405     uncoupl_optic_haze=.true. ! default value - true as long as the microphysics is bugged
406     call getin_p("uncoupl_optic_haze",uncoupl_optic_haze)
407     write(*,*)" uncoupl_optic_haze = ",uncoupl_optic_haze
408
409     write(*,*) "Pressure level of aer. production (Pa) ?"
410     p_prod=1.0 ! default value
411     call getin_p("p_prod",p_prod)
412     write(*,*)" p_prod = ",p_prod
413     
414     write(*,*) "Aerosol production rate (kg.m-2.s-1) ?"
415     tx_prod=3.5e-13 ! default value
416     call getin_p("tx_prod",tx_prod)
417     write(*,*)" tx_prod = ",tx_prod
418
419     write(*,*) "Equivalent radius production (m) ?"
420     rc_prod=2.0e-8 ! default value
421     call getin_p("rc_prod",rc_prod)
422     write(*,*)" rhc_prod = ",rc_prod
423
424     write(*,*) "Radius of air (nitrogen) molecule (m) ?"
425     air_rad=1.75e-10 ! default value
426     call getin_p("air_rad",air_rad)
427     write(*,*)" air_rad = ",air_rad
428
429     write(*,*) "Path to microphys. config file ?"
430     config_mufi='datagcm/microphysics/config.cfg' ! default value
431     call getin_p("config_mufi",config_mufi)
432     write(*,*)" config_mufi = ",config_mufi
433
434! Soil model
435     write(*,*)"Number of sub-surface layers for soil scheme?"
436     ! default value of nsoilmx set in comsoil_h
437     call getin_p("nsoilmx",nsoilmx)
438     write(*,*)" nsoilmx=",nsoilmx
439     
440     write(*,*)"Thickness of topmost soil layer (m)?"
441     ! default value of lay1_soil set in comsoil_h
442     call getin_p("lay1_soil",lay1_soil)
443     write(*,*)" lay1_soil = ",lay1_soil
444     
445     write(*,*)"Coefficient for soil layer thickness distribution?"
446     ! default value of alpha_soil set in comsoil_h
447     call getin_p("alpha_soil",alpha_soil)
448     write(*,*)" alpha_soil = ",alpha_soil
449
450     write(*,*)"Remove lower boundary?"
451     nosurf=.false.
452     call getin_p("nosurf",nosurf)
453     write(*,*)" nosurf = ",nosurf
454
455! Tests of incompatibility:
456     if (nosurf.and.callsoil) then
457       print*,'nosurf not compatible with soil scheme!'
458       print*,'... got to make a choice!'
459       call abort
460     endif
461
462     write(*,*)"Add an internal heat flux?", &
463                   "... matters only if callsoil=F"
464     intheat=0.
465     call getin_p("intheat",intheat)
466     write(*,*)" intheat = ",intheat
467
468     write(*,*)"Use Newtonian cooling for radiative transfer?"
469     newtonian=.false.
470     call getin_p("newtonian",newtonian)
471     write(*,*)" newtonian = ",newtonian
472
473! Tests of incompatibility:
474     if (newtonian.and.corrk) then
475       print*,'newtonian not compatible with correlated-k!'
476       call abort
477     endif
478     if (newtonian.and.calladj) then
479       print*,'newtonian not compatible with adjustment!'
480       call abort
481     endif
482     if (newtonian.and.calldifv) then
483       print*,'newtonian not compatible with a boundary layer!'
484       call abort
485     endif
486
487     write(*,*)"Test physics timescale in 1D?"
488     testradtimes=.false.
489     call getin_p("testradtimes",testradtimes)
490     write(*,*)" testradtimes = ",testradtimes
491
492! Test of incompatibility:
493! if testradtimes used, we must be in 1D
494     if (testradtimes.and.(ngrid.gt.1)) then
495       print*,'testradtimes can only be used in 1D!'
496       call abort
497     endif
498
499     write(*,*)"Which star?"
500     startype=1 ! default value = Sol
501     call getin_p("startype",startype)
502     write(*,*)" startype = ",startype
503
504     write(*,*)"Value of stellar flux at 1 AU?"
505     Fat1AU=1356.0 ! default value = Sol today
506     call getin_p("Fat1AU",Fat1AU)
507     write(*,*)" Fat1AU = ",Fat1AU
508
509     write(*,*) "Does user want to force cpp and mugaz?"
510     force_cpp=.false. ! default value
511     call getin_p("force_cpp",force_cpp)
512     write(*,*) " force_cpp = ",force_cpp
513
514     if (force_cpp) then
515       mugaz = -99999.
516       PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
517       call getin_p("mugaz",mugaz)
518       IF (mugaz.eq.-99999.) THEN
519           PRINT *, "mugaz must be set if force_cpp = T"
520           STOP
521       ELSE
522           write(*,*) "mugaz=",mugaz
523       ENDIF
524       cpp = -99999.
525       PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
526       call getin_p("cpp",cpp)
527       IF (cpp.eq.-99999.) THEN
528           PRINT *, "cpp must be set if force_cpp = T"
529           STOP
530       ELSE
531           write(*,*) "cpp=",cpp
532       ENDIF
533!         else
534!           mugaz=8.314*1000./pr
535     endif ! of if (force_cpp)
536     
537     
538     call su_gases(nlayer,tracer)     
539     call calc_cpp_mugaz
540     
541     
542     PRINT*,'--------------------------------------------'
543     PRINT*
544     PRINT*
545  ELSE
546     write(*,*)
547     write(*,*) 'Cannot read file callphys.def. Is it here ?'
548     stop
549  ENDIF ! of IF(iscallphys)
550
551  PRINT*
552  PRINT*,'inifis: daysec',daysec
553  PRINT*
554  PRINT*,'inifis: The radiative transfer is computed:'
555  PRINT*,'           each ',iradia,' physical time-step'
556  PRINT*,'        or each ',iradia*dtphys,' seconds'
557  PRINT*
558  PRINT*,'inifis: Chemistry is computed:'
559  PRINT*,'           each ',ichim,' physical time-step'
560  PRINT*,'        or each ',ichim*dtphys,' seconds'
561  PRINT*
562
563!-----------------------------------------------------------------------
564!     Some more initialization:
565!     ------------------------
566
567  ! Initializations for comgeomfi_h
568  totarea=SSUM(ngrid,parea,1)
569  call planetwide_sumval(parea,totarea_planet)
570
571  !! those are defined in comdiurn_h.F90
572  IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
573  IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
574  IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
575  IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
576
577  DO ig=1,ngrid
578     sinlat(ig)=sin(plat(ig))
579     coslat(ig)=cos(plat(ig))
580     sinlon(ig)=sin(plon(ig))
581     coslon(ig)=cos(plon(ig))
582  ENDDO
583
584  ! initialize variables in radinc_h
585  call ini_radinc_h(nlayer)
586 
587  ! allocate "comsoil_h" arrays
588  call ini_comsoil_h(ngrid)
589     
590  END SUBROUTINE inifis
591
592END MODULE inifis_mod
Note: See TracBrowser for help on using the repository browser.