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

Last change on this file since 1668 was 1648, checked in by jvatant, 9 years ago

Modifications to custom radiative transfer to Titan
+ Enables an altitude dependant gfrac for CIA computations

-> many radical changes in su_gases and co ..
-> read vertical CH4 profile with call_profilgases
-> Now you need a 'profile.def' that I will add in the deftank

+ Added interpolate CIA routines for CH4
+ Added temporary mean aerosol profile opacity routine (disr_haze)

File size: 16.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 radinc_h, only: ini_radinc_h, naerkind
12  use radcommon_h, only: ini_radcommon_h
13  use datafile_mod, only: datadir
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
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 datafile_mod, only: datadir
62  use ioipsl_getin_p_mod, only: getin_p
63  IMPLICIT NONE
64
65
66
67  REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep
68  INTEGER,INTENT(IN) :: nday
69  INTEGER,INTENT(IN) :: ngrid,nlayer,nq
70  REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid)
71  integer,intent(in) :: day_ini
72  INTEGER ig,ierr
73 
74  EXTERNAL iniorbit,orbite
75  EXTERNAL SSUM
76  REAL SSUM
77 
78  ! initialize constants in comcstfi_mod
79  rad=prad
80  cpp=pcpp
81  g=pg
82  r=pr
83  rcp=r/cpp
84  pi=2.*asin(1.)
85  avocado = 6.02214179e23   ! added by RW
86
87  ! Initialize some "temporal and calendar" related variables
88  CALL init_time(day_ini,pdaysec,nday,ptimestep)
89
90  ! read in some parameters from "run.def" for physics,
91  ! or shared between dynamics and physics.
92  call getin_p("ecritphy",ecritphy) ! frequency of outputs in physics,
93                                    ! in dynamical steps
94  call getin_p("day_step",day_step) ! number of dynamical steps per day
95  call getin_p("iphysiq",iphysiq) ! call physics every iphysiq dyn step
96
97! --------------------------------------------------------------
98!  Reading the "callphys.def" file controlling some key options
99! --------------------------------------------------------------
100     
101  ! check that 'callphys.def' file is around
102  OPEN(99,file='callphys.def',status='old',form='formatted',iostat=ierr)
103  CLOSE(99)
104  IF(ierr.EQ.0) iscallphys=.true. !iscallphys initialised as false in callkeys_mod module
105     
106!!!      IF(ierr.EQ.0) THEN
107  IF(iscallphys) THEN
108     PRINT*
109     PRINT*
110     PRINT*,'--------------------------------------------'
111     PRINT*,' inifis: Parametres pour la physique (callphys.def)'
112     PRINT*,'--------------------------------------------'
113
114     write(*,*) "Directory where external input files are:"
115     ! default 'datadir' is set in "datadir_mod"
116     call getin_p("datadir",datadir) ! default path
117     write(*,*) " datadir = ",trim(datadir)
118
119     write(*,*) "Run with or without tracer transport ?"
120     tracer=.false. ! default value
121     call getin_p("tracer",tracer)
122     write(*,*) " tracer = ",tracer
123
124     write(*,*) "Run with or without atm mass update ", &
125            " due to tracer evaporation/condensation?"
126     mass_redistrib=.false. ! default value
127     call getin_p("mass_redistrib",mass_redistrib)
128     write(*,*) " mass_redistrib = ",mass_redistrib
129
130     write(*,*) "Diurnal cycle ?"
131     write(*,*) "(if diurnal=false, diurnal averaged solar heating)"
132     diurnal=.true. ! default value
133     call getin_p("diurnal",diurnal)
134     write(*,*) " diurnal = ",diurnal
135
136     write(*,*) "Seasonal cycle ?"
137     write(*,*) "(if season=false, Ls stays constant, to value ", &
138         "set in 'start'"
139     season=.true. ! default value
140     call getin_p("season",season)
141     write(*,*) " season = ",season
142
143     write(*,*) "Tidally resonant rotation ?"
144     tlocked=.false. ! default value
145     call getin_p("tlocked",tlocked)
146     write(*,*) "tlocked = ",tlocked
147
148     write(*,*) "Saturn ring shadowing ?"
149     rings_shadow = .false.
150     call getin_p("rings_shadow", rings_shadow)
151     write(*,*) "rings_shadow = ", rings_shadow
152         
153     write(*,*) "Compute latitude-dependent gravity field?"
154     oblate = .false.
155     call getin_p("oblate", oblate)
156     write(*,*) "oblate = ", oblate
157
158     write(*,*) "Flattening of the planet (a-b)/a "
159     flatten = 0.0
160     call getin_p("flatten", flatten)
161     write(*,*) "flatten = ", flatten
162         
163
164     write(*,*) "Needed if oblate=.true.: J2"
165     J2 = 0.0
166     call getin_p("J2", J2)
167     write(*,*) "J2 = ", J2
168         
169     write(*,*) "Needed if oblate=.true.: Planet mass (*1e24 kg)"
170     MassPlanet = 0.0
171     call getin_p("MassPlanet", MassPlanet)
172     write(*,*) "MassPlanet = ", MassPlanet         
173
174     write(*,*) "Needed if oblate=.true.: Planet mean radius (m)"
175     Rmean = 0.0
176     call getin_p("Rmean", Rmean)
177     write(*,*) "Rmean = ", Rmean
178         
179! Test of incompatibility:
180! if tlocked, then diurnal should be false
181     if (tlocked.and.diurnal) then
182       print*,'If diurnal=true, we should turn off tlocked.'
183       stop
184     endif
185
186     write(*,*) "Tidal resonance ratio ?"
187     nres=0          ! default value
188     call getin_p("nres",nres)
189     write(*,*) "nres = ",nres
190
191     write(*,*) "Write some extra output to the screen ?"
192     lwrite=.false. ! default value
193     call getin_p("lwrite",lwrite)
194     write(*,*) " lwrite = ",lwrite
195
196     write(*,*) "Save statistics in file stats.nc ?"
197     callstats=.true. ! default value
198     call getin_p("callstats",callstats)
199     write(*,*) " callstats = ",callstats
200
201     write(*,*) "Test energy conservation of model physics ?"
202     enertest=.false. ! default value
203     call getin_p("enertest",enertest)
204     write(*,*) " enertest = ",enertest
205
206     write(*,*) "Check to see if cpp values used match gases.def ?"
207     check_cpp_match=.true. ! default value
208     call getin_p("check_cpp_match",check_cpp_match)
209     write(*,*) " check_cpp_match = ",check_cpp_match
210
211     write(*,*) "call radiative transfer ?"
212     callrad=.true. ! default value
213     call getin_p("callrad",callrad)
214     write(*,*) " callrad = ",callrad
215
216     write(*,*) "call correlated-k radiative transfer ?"
217     corrk=.true. ! default value
218     call getin_p("corrk",corrk)
219     write(*,*) " corrk = ",corrk
220
221     write(*,*) "prohibit calculations outside corrk T grid?"
222     strictboundcorrk=.true. ! default value
223     call getin_p("strictboundcorrk",strictboundcorrk)
224     write(*,*) "strictboundcorrk = ",strictboundcorrk
225
226     write(*,*) "call gaseous absorption in the visible bands?", &
227                    "(matters only if callrad=T)"
228     callgasvis=.false. ! default value
229     call getin_p("callgasvis",callgasvis)
230     write(*,*) " callgasvis = ",callgasvis
231       
232     write(*,*) "call continuum opacities in radiative transfer ?", &
233                    "(matters only if callrad=T)"
234     continuum=.true. ! default value
235     call getin_p("continuum",continuum)
236     write(*,*) " continuum = ",continuum
237 
238     write(*,*) "call turbulent vertical diffusion ?"
239     calldifv=.true. ! default value
240     call getin_p("calldifv",calldifv)
241     write(*,*) " calldifv = ",calldifv
242
243     write(*,*) "use turbdiff instead of vdifc ?"
244     UseTurbDiff=.true. ! default value
245     call getin_p("UseTurbDiff",UseTurbDiff)
246     write(*,*) " UseTurbDiff = ",UseTurbDiff
247
248     write(*,*) "call convective adjustment ?"
249     calladj=.true. ! default value
250     call getin_p("calladj",calladj)
251     write(*,*) " calladj = ",calladj
252
253     write(*,*) "Radiative timescale for Newtonian cooling ?"
254     tau_relax=30. ! default value
255     call getin_p("tau_relax",tau_relax)
256     write(*,*) " tau_relax = ",tau_relax
257     tau_relax=tau_relax*24*3600 ! convert Earth days --> seconds
258
259     write(*,*)"call thermal conduction in the soil ?"
260     callsoil=.true. ! default value
261     call getin_p("callsoil",callsoil)
262     write(*,*) " callsoil = ",callsoil
263         
264     write(*,*)"Rad transfer is computed every iradia", &
265                   " physical timestep"
266     iradia=1 ! default value
267     call getin_p("iradia",iradia)
268     write(*,*)" iradia = ",iradia
269       
270     write(*,*)"Rayleigh scattering ?"
271     rayleigh=.false.
272     call getin_p("rayleigh",rayleigh)
273     write(*,*)" rayleigh = ",rayleigh
274
275     write(*,*) "Use blackbody for stellar spectrum ?"
276     stelbbody=.false. ! default value
277     call getin_p("stelbbody",stelbbody)
278     write(*,*) " stelbbody = ",stelbbody
279
280     write(*,*) "Stellar blackbody temperature ?"
281     stelTbb=5800.0 ! default value
282     call getin_p("stelTbb",stelTbb)
283     write(*,*) " stelTbb = ",stelTbb
284
285     write(*,*)"Output mean OLR in 1D?"
286     meanOLR=.false.
287     call getin_p("meanOLR",meanOLR)
288     write(*,*)" meanOLR = ",meanOLR
289
290     write(*,*)"Output spectral OLR in 3D?"
291     specOLR=.false.
292     call getin_p("specOLR",specOLR)
293     write(*,*)" specOLR = ",specOLR
294
295     write(*,*)"Uniform absorption in radiative transfer?"
296     graybody=.false.
297     call getin_p("graybody",graybody)
298     write(*,*)" graybody = ",graybody
299
300! Soil model
301     write(*,*)"Number of sub-surface layers for soil scheme?"
302     ! default value of nsoilmx set in comsoil_h
303     call getin_p("nsoilmx",nsoilmx)
304     write(*,*)" nsoilmx=",nsoilmx
305     
306     write(*,*)"Thickness of topmost soil layer (m)?"
307     ! default value of lay1_soil set in comsoil_h
308     call getin_p("lay1_soil",lay1_soil)
309     write(*,*)" lay1_soil = ",lay1_soil
310     
311     write(*,*)"Coefficient for soil layer thickness distribution?"
312     ! default value of alpha_soil set in comsoil_h
313     call getin_p("alpha_soil",alpha_soil)
314     write(*,*)" alpha_soil = ",alpha_soil
315
316     write(*,*)"Remove lower boundary?"
317     nosurf=.false.
318     call getin_p("nosurf",nosurf)
319     write(*,*)" nosurf = ",nosurf
320
321! Tests of incompatibility:
322     if (nosurf.and.callsoil) then
323       print*,'nosurf not compatible with soil scheme!'
324       print*,'... got to make a choice!'
325       call abort
326     endif
327
328     write(*,*)"Add an internal heat flux?", &
329                   "... matters only if callsoil=F"
330     intheat=0.
331     call getin_p("intheat",intheat)
332     write(*,*)" intheat = ",intheat
333
334     write(*,*)"Use Newtonian cooling for radiative transfer?"
335     newtonian=.false.
336     call getin_p("newtonian",newtonian)
337     write(*,*)" newtonian = ",newtonian
338
339! Tests of incompatibility:
340     if (newtonian.and.corrk) then
341       print*,'newtonian not compatible with correlated-k!'
342       call abort
343     endif
344     if (newtonian.and.calladj) then
345       print*,'newtonian not compatible with adjustment!'
346       call abort
347     endif
348     if (newtonian.and.calldifv) then
349       print*,'newtonian not compatible with a boundary layer!'
350       call abort
351     endif
352
353     write(*,*)"Test physics timescale in 1D?"
354     testradtimes=.false.
355     call getin_p("testradtimes",testradtimes)
356     write(*,*)" testradtimes = ",testradtimes
357
358! Test of incompatibility:
359! if testradtimes used, we must be in 1D
360     if (testradtimes.and.(ngrid.gt.1)) then
361       print*,'testradtimes can only be used in 1D!'
362       call abort
363     endif
364
365     write(*,*)"Default planetary temperature?"
366     tplanet=215.0
367     call getin_p("tplanet",tplanet)
368     write(*,*)" tplanet = ",tplanet
369
370     write(*,*)"Which star?"
371     startype=1 ! default value = Sol
372     call getin_p("startype",startype)
373     write(*,*)" startype = ",startype
374
375     write(*,*)"Value of stellar flux at 1 AU?"
376     Fat1AU=1356.0 ! default value = Sol today
377     call getin_p("Fat1AU",Fat1AU)
378     write(*,*)" Fat1AU = ",Fat1AU
379
380
381! TRACERS:
382
383!         write(*,*)"Number of radiatively active aerosols:"
384!         naerkind=0. ! default value
385!         call getin_p("naerkind",naerkind)
386!         write(*,*)" naerkind = ",naerkind
387
388 
389!=================================
390
391     write(*,*)"Radiatively active two-layer aersols?"
392     aeroback2lay=.false.     ! default value
393     call getin_p("aeroback2lay",aeroback2lay)
394     write(*,*)" aeroback2lay = ",aeroback2lay
395
396     write(*,*)"TWOLAY AEROSOL: total optical depth ", &
397                    "in the tropospheric layer (visible)"
398     obs_tau_col_tropo=8.D0
399     call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
400     write(*,*)" obs_tau_col_tropo = ",obs_tau_col_tropo
401
402     write(*,*)"TWOLAY AEROSOL: total optical depth ", &
403                    "in the stratospheric layer (visible)"
404     obs_tau_col_strato=0.08D0
405     call getin_p("obs_tau_col_strato",obs_tau_col_strato)
406     write(*,*)" obs_tau_col_strato = ",obs_tau_col_strato
407
408     write(*,*)"TWOLAY AEROSOL: pres_bottom_tropo? in pa"
409     pres_bottom_tropo=66000.0
410     call getin_p("pres_bottom_tropo",pres_bottom_tropo)
411     write(*,*)" pres_bottom_tropo = ",pres_bottom_tropo
412
413     write(*,*)"TWOLAY AEROSOL: pres_top_tropo? in pa"
414     pres_top_tropo=18000.0
415     call getin_p("pres_top_tropo",pres_top_tropo)
416     write(*,*)" pres_top_tropo = ",pres_top_tropo
417
418     write(*,*)"TWOLAY AEROSOL: pres_bottom_strato? in pa"
419     pres_bottom_strato=2000.0
420     call getin_p("pres_bottom_strato",pres_bottom_strato)
421     write(*,*)" pres_bottom_strato = ",pres_bottom_strato
422
423     write(*,*)"TWOLAY AEROSOL: pres_top_strato? in pa"
424     pres_top_strato=100.0
425     call getin_p("pres_top_strato",pres_top_strato)
426     write(*,*)" pres_top_strato = ",pres_top_strato
427
428     write(*,*)"TWOLAY AEROSOL: particle size in the ", &
429                    "tropospheric layer, in meters"
430     size_tropo=2.e-6
431     call getin_p("size_tropo",size_tropo)
432     write(*,*)" size_tropo = ",size_tropo
433
434     write(*,*)"TWOLAY AEROSOL: particle size in the ", &
435                    "stratospheric layer, in meters"
436     size_strato=1.e-7
437     call getin_p("size_strato",size_strato)
438     write(*,*)" size_strato = ",size_strato
439
440!=================================
441
442     write(*,*) "Gravitationnal sedimentation ?"
443     sedimentation=.false. ! default value
444     call getin_p("sedimentation",sedimentation)
445     write(*,*) " sedimentation = ",sedimentation       
446
447     write(*,*) "Does user want to force cpp and mugaz?"
448     force_cpp=.false. ! default value
449     call getin_p("force_cpp",force_cpp)
450     write(*,*) " force_cpp = ",force_cpp
451
452     if (force_cpp) then
453       mugaz = -99999.
454       PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
455       call getin_p("mugaz",mugaz)
456       IF (mugaz.eq.-99999.) THEN
457           PRINT *, "mugaz must be set if force_cpp = T"
458           STOP
459       ELSE
460           write(*,*) "mugaz=",mugaz
461       ENDIF
462       cpp = -99999.
463       PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
464       call getin_p("cpp",cpp)
465       IF (cpp.eq.-99999.) THEN
466           PRINT *, "cpp must be set if force_cpp = T"
467           STOP
468       ELSE
469           write(*,*) "cpp=",cpp
470       ENDIF
471!         else
472!           mugaz=8.314*1000./pr
473     endif ! of if (force_cpp)
474     
475     
476     call su_gases(nlayer,tracer)     
477     call calc_cpp_mugaz
478     
479     
480     PRINT*,'--------------------------------------------'
481     PRINT*
482     PRINT*
483  ELSE
484     write(*,*)
485     write(*,*) 'Cannot read file callphys.def. Is it here ?'
486     stop
487  ENDIF ! of IF(iscallphys)
488
489  PRINT*
490  PRINT*,'inifis: daysec',daysec
491  PRINT*
492  PRINT*,'inifis: The radiative transfer is computed:'
493  PRINT*,'           each ',iradia,' physical time-step'
494  PRINT*,'        or each ',iradia*dtphys,' seconds'
495  PRINT*
496
497
498!-----------------------------------------------------------------------
499!     Some more initialization:
500!     ------------------------
501
502  ! Initializations for comgeomfi_h
503  totarea=SSUM(ngrid,parea,1)
504  call planetwide_sumval(parea,totarea_planet)
505
506  !! those are defined in comdiurn_h.F90
507  IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
508  IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
509  IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
510  IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
511
512  DO ig=1,ngrid
513     sinlat(ig)=sin(plat(ig))
514     coslat(ig)=cos(plat(ig))
515     sinlon(ig)=sin(plon(ig))
516     coslon(ig)=cos(plon(ig))
517  ENDDO
518
519  call ini_radinc_h(nlayer)
520 
521  ! allocate "radcommon_h" arrays
522  call ini_radcommon_h()
523
524  ! allocate "comsoil_h" arrays
525  call ini_comsoil_h(ngrid)
526     
527  END SUBROUTINE inifis
528
529END MODULE inifis_mod
Note: See TracBrowser for help on using the repository browser.