source: trunk/LMDZ.GENERIC/libf/phystd/inifis.F @ 1123

Last change on this file since 1123 was 1026, checked in by sglmd, 12 years ago

Added a flexible, 2-layer aerosol scenario (initially for Saturn, but can be simply tuned). Called aeroback2lay.

File size: 19.8 KB
Line 
1      SUBROUTINE inifis(ngrid,nlayer,
2     $           day_ini,pdaysec,ptimestep,
3     $           plat,plon,parea,
4     $           prad,pg,pr,pcpp)
5
6      use radinc_h, only : naerkind
7      use datafile_mod, only: datadir
8      use comdiurn_h
9
10      !! to be conservative, but why this include here?
11      use surfdat_h
12      use comsaison_h
13
14      USE comgeomfi_h
15
16!=======================================================================
17!
18!   purpose:
19!   -------
20!
21!   Initialisation for the physical parametrisations of the LMD
22!   Generic Model.
23!
24!   author: Frederic Hourdin 15 / 10 /93
25!   -------
26!   modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
27!             Ehouarn Millour (oct. 2008) tracers are now identified
28!              by their names and may not be contiguously
29!              stored in the q(:,:,:,:) array
30!             E.M. (june 2009) use getin routine to load parameters
31!
32!
33!   arguments:
34!   ----------
35!
36!   input:
37!   ------
38!
39!    ngrid                 Size of the horizontal grid.
40!                          All internal loops are performed on that grid.
41!    nlayer                Number of vertical layers.
42!    pdayref               Day of reference for the simulation
43!    pday                  Number of days counted from the North. Spring
44!                          equinoxe.
45!
46!=======================================================================
47!
48!-----------------------------------------------------------------------
49!   declarations:
50!   -------------
51      use datafile_mod, only: datadir
52! to use  'getin'
53      USE ioipsl_getincom
54      IMPLICIT NONE
55#include "dimensions.h"
56#include "dimphys.h"
57#include "planete.h"
58#include "comcstfi.h"
59#include "callkeys.h"
60
61
62
63      REAL prad,pg,pr,pcpp,pdaysec,ptimestep
64 
65      INTEGER ngrid,nlayer
66      REAL plat(ngrid),plon(ngrid),parea(ngrid)
67      integer day_ini
68      INTEGER ig,ierr
69 
70      EXTERNAL iniorbit,orbite
71      EXTERNAL SSUM
72      REAL SSUM
73 
74      CHARACTER ch1*12
75      CHARACTER ch80*80
76
77      logical chem, h2o
78      logical :: parameter, doubleq=.false.
79
80      real psurf,pN2 ! added by RW for Gliese 581d N2+CO2
81
82      rad=prad
83      daysec=pdaysec
84      dtphys=ptimestep
85      cpp=pcpp
86      g=pg
87      r=pr
88      rcp=r/cpp
89
90      avocado = 6.02214179e23   ! added by RW
91
92      cursor = 1 ! added by AS in dimphys. 1 for sequential runs.
93
94! --------------------------------------------------------
95!     The usual Tests
96!     --------------
97      IF (nlayer.NE.nlayermx) THEN
98         PRINT*,'STOP in inifis'
99         PRINT*,'Probleme de dimensions :'
100         PRINT*,'nlayer     = ',nlayer
101         PRINT*,'nlayermx   = ',nlayermx
102         STOP
103      ENDIF
104
105! --------------------------------------------------------------
106!  Reading the "callphys.def" file controlling some key options
107! --------------------------------------------------------------
108     
109      ! check that 'callphys.def' file is around
110      OPEN(99,file='callphys.def',status='old',form='formatted'
111     &     ,iostat=ierr)
112      CLOSE(99)
113     
114      IF(ierr.EQ.0) 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("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("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("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("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("season",season)
148         write(*,*) " season = ",season
149
150         write(*,*) "Tidally resonant rotation ?"
151         tlocked=.false. ! default value
152         call getin("tlocked",tlocked)
153         write(*,*) "tlocked = ",tlocked
154
155! Test of incompatibility:
156! if tlocked, then diurnal should be false
157         if (tlocked.and.diurnal) then
158           print*,'If diurnal=true, we should turn off tlocked.'
159           stop
160         endif
161
162         write(*,*) "Tidal resonance ratio ?"
163         nres=0          ! default value
164         call getin("nres",nres)
165         write(*,*) "nres = ",nres
166
167         write(*,*) "Write some extra output to the screen ?"
168         lwrite=.false. ! default value
169         call getin("lwrite",lwrite)
170         write(*,*) " lwrite = ",lwrite
171
172         write(*,*) "Save statistics in file stats.nc ?"
173         callstats=.true. ! default value
174         call getin("callstats",callstats)
175         write(*,*) " callstats = ",callstats
176
177         write(*,*) "Test energy conservation of model physics ?"
178         enertest=.false. ! default value
179         call getin("enertest",enertest)
180         write(*,*) " enertest = ",enertest
181
182         write(*,*) "Check to see if cpp values used match gases.def ?"
183         check_cpp_match=.true. ! default value
184         call getin("check_cpp_match",check_cpp_match)
185         write(*,*) " check_cpp_match = ",check_cpp_match
186
187         write(*,*) "call radiative transfer ?"
188         callrad=.true. ! default value
189         call getin("callrad",callrad)
190         write(*,*) " callrad = ",callrad
191
192         write(*,*) "call correlated-k radiative transfer ?"
193         corrk=.true. ! default value
194         call getin("corrk",corrk)
195         write(*,*) " corrk = ",corrk
196
197         write(*,*) "call gaseous absorption in the visible bands?",
198     &              "(matters only if callrad=T)"
199         callgasvis=.false. ! default value
200         call getin("callgasvis",callgasvis)
201         write(*,*) " callgasvis = ",callgasvis
202       
203         write(*,*) "call continuum opacities in radiative transfer ?",
204     &              "(matters only if callrad=T)"
205         continuum=.true. ! default value
206         call getin("continuum",continuum)
207         write(*,*) " continuum = ",continuum
208
209         write(*,*) "use analytic function for H2O continuum ?"
210         H2Ocont_simple=.false. ! default value
211         call getin("H2Ocont_simple",H2Ocont_simple)
212         write(*,*) " H2Ocont_simple = ",H2Ocont_simple
213 
214         write(*,*) "call turbulent vertical diffusion ?"
215         calldifv=.true. ! default value
216         call getin("calldifv",calldifv)
217         write(*,*) " calldifv = ",calldifv
218
219         write(*,*) "use turbdiff instead of vdifc ?"
220         UseTurbDiff=.true. ! default value
221         call getin("UseTurbDiff",UseTurbDiff)
222         write(*,*) " UseTurbDiff = ",UseTurbDiff
223
224         write(*,*) "call convective adjustment ?"
225         calladj=.true. ! default value
226         call getin("calladj",calladj)
227         write(*,*) " calladj = ",calladj
228
229         write(*,*) "call CO2 condensation ?"
230         co2cond=.false. ! default value
231         call getin("co2cond",co2cond)
232         write(*,*) " co2cond = ",co2cond
233! Test of incompatibility
234         if (co2cond.and.(.not.tracer)) then
235            print*,'We need a CO2 ice tracer to condense CO2'
236            call abort
237         endif
238 
239         write(*,*) "CO2 supersaturation level ?"
240         co2supsat=1.0 ! default value
241         call getin("co2supsat",co2supsat)
242         write(*,*) " co2supsat = ",co2supsat
243
244         write(*,*) "Radiative timescale for Newtonian cooling ?"
245         tau_relax=30. ! default value
246         call getin("tau_relax",tau_relax)
247         write(*,*) " tau_relax = ",tau_relax
248         tau_relax=tau_relax*24*3600 ! convert Earth days --> seconds
249
250         write(*,*)"call thermal conduction in the soil ?"
251         callsoil=.true. ! default value
252         call getin("callsoil",callsoil)
253         write(*,*) " callsoil = ",callsoil
254         
255         write(*,*)"Rad transfer is computed every iradia",
256     &             " physical timestep"
257         iradia=1 ! default value
258         call getin("iradia",iradia)
259         write(*,*)" iradia = ",iradia
260       
261         write(*,*)"Rayleigh scattering ?"
262         rayleigh=.false.
263         call getin("rayleigh",rayleigh)
264         write(*,*)" rayleigh = ",rayleigh
265
266         write(*,*) "Use blackbody for stellar spectrum ?"
267         stelbbody=.false. ! default value
268         call getin("stelbbody",stelbbody)
269         write(*,*) " stelbbody = ",stelbbody
270
271         write(*,*) "Stellar blackbody temperature ?"
272         stelTbb=5800.0 ! default value
273         call getin("stelTbb",stelTbb)
274         write(*,*) " stelTbb = ",stelTbb
275
276         write(*,*)"Output mean OLR in 1D?"
277         meanOLR=.false.
278         call getin("meanOLR",meanOLR)
279         write(*,*)" meanOLR = ",meanOLR
280
281         write(*,*)"Output spectral OLR in 3D?"
282         specOLR=.false.
283         call getin("specOLR",specOLR)
284         write(*,*)" specOLR = ",specOLR
285
286         write(*,*)"Operate in kastprof mode?"
287         kastprof=.false.
288         call getin("kastprof",kastprof)
289         write(*,*)" kastprof = ",kastprof
290
291         write(*,*)"Uniform absorption in radiative transfer?"
292         graybody=.false.
293         call getin("graybody",graybody)
294         write(*,*)" graybody = ",graybody
295
296! Test of incompatibility:
297! if kastprof used, we must be in 1D
298         if (kastprof.and.(ngrid.gt.1)) then
299           print*,'kastprof can only be used in 1D!'
300           call abort
301         endif
302
303         write(*,*)"Stratospheric temperature for kastprof mode?"
304         Tstrat=167.0
305         call getin("Tstrat",Tstrat)
306         write(*,*)" Tstrat = ",Tstrat
307
308         write(*,*)"Remove lower boundary?"
309         nosurf=.false.
310         call getin("nosurf",nosurf)
311         write(*,*)" nosurf = ",nosurf
312
313! Tests of incompatibility:
314         if (nosurf.and.callsoil) then
315           print*,'nosurf not compatible with soil scheme!'
316           print*,'... got to make a choice!'
317           call abort
318         endif
319
320         write(*,*)"Add an internal heat flux?",
321     .             "... matters only if callsoil=F"
322         intheat=0.
323         call getin("intheat",intheat)
324         write(*,*)" intheat = ",intheat
325
326         write(*,*)"Use Newtonian cooling for radiative transfer?"
327         newtonian=.false.
328         call getin("newtonian",newtonian)
329         write(*,*)" newtonian = ",newtonian
330
331! Tests of incompatibility:
332         if (newtonian.and.corrk) then
333           print*,'newtonian not compatible with correlated-k!'
334           call abort
335         endif
336         if (newtonian.and.calladj) then
337           print*,'newtonian not compatible with adjustment!'
338           call abort
339         endif
340         if (newtonian.and.calldifv) then
341           print*,'newtonian not compatible with a boundary layer!'
342           call abort
343         endif
344
345         write(*,*)"Test physics timescale in 1D?"
346         testradtimes=.false.
347         call getin("testradtimes",testradtimes)
348         write(*,*)" testradtimes = ",testradtimes
349
350! Test of incompatibility:
351! if testradtimes used, we must be in 1D
352         if (testradtimes.and.(ngrid.gt.1)) then
353           print*,'testradtimes can only be used in 1D!'
354           call abort
355         endif
356
357         write(*,*)"Default planetary temperature?"
358         tplanet=215.0
359         call getin("tplanet",tplanet)
360         write(*,*)" tplanet = ",tplanet
361
362         write(*,*)"Which star?"
363         startype=1 ! default value = Sol
364         call getin("startype",startype)
365         write(*,*)" startype = ",startype
366
367         write(*,*)"Value of stellar flux at 1 AU?"
368         Fat1AU=1356.0 ! default value = Sol today
369         call getin("Fat1AU",Fat1AU)
370         write(*,*)" Fat1AU = ",Fat1AU
371
372
373! TRACERS:
374
375         write(*,*)"Varying H2O cloud fraction?"
376         CLFvarying=.false.     ! default value
377         call getin("CLFvarying",CLFvarying)
378         write(*,*)" CLFvarying = ",CLFvarying
379
380         write(*,*)"Value of fixed H2O cloud fraction?"
381         CLFfixval=1.0                ! default value
382         call getin("CLFfixval",CLFfixval)
383         write(*,*)" CLFfixval = ",CLFfixval
384
385         write(*,*)"fixed radii for Cloud particles?"
386         radfixed=.false. ! default value
387         call getin("radfixed",radfixed)
388         write(*,*)" radfixed = ",radfixed
389
390         if(kastprof)then
391            radfixed=.true.
392         endif 
393
394         write(*,*)"Number mixing ratio of CO2 ice particles:"
395         Nmix_co2=1.e6 ! default value
396         call getin("Nmix_co2",Nmix_co2)
397         write(*,*)" Nmix_co2 = ",Nmix_co2
398
399!         write(*,*)"Number of radiatively active aerosols:"
400!         naerkind=0. ! default value
401!         call getin("naerkind",naerkind)
402!         write(*,*)" naerkind = ",naerkind
403
404         write(*,*)"Opacity of dust (if used):"
405         dusttau=0. ! default value
406         call getin("dusttau",dusttau)
407         write(*,*)" dusttau = ",dusttau
408
409         write(*,*)"Radiatively active CO2 aerosols?"
410         aeroco2=.false.     ! default value
411         call getin("aeroco2",aeroco2)
412         write(*,*)" aeroco2 = ",aeroco2
413
414         write(*,*)"Fixed CO2 aerosol distribution?"
415         aerofixco2=.false.     ! default value
416         call getin("aerofixco2",aerofixco2)
417         write(*,*)" aerofixco2 = ",aerofixco2
418
419         write(*,*)"Radiatively active water ice?"
420         aeroh2o=.false.     ! default value
421         call getin("aeroh2o",aeroh2o)
422         write(*,*)" aeroh2o = ",aeroh2o
423
424         write(*,*)"Fixed H2O aerosol distribution?"
425         aerofixh2o=.false.     ! default value
426         call getin("aerofixh2o",aerofixh2o)
427         write(*,*)" aerofixh2o = ",aerofixh2o
428
429         write(*,*)"Radiatively active sulfuric acid aersols?"
430         aeroh2so4=.false.     ! default value
431         call getin("aeroh2so4",aeroh2so4)
432         write(*,*)" aeroh2so4 = ",aeroh2so4
433         
434         write(*,*)"Radiatively active two-layer aersols?"
435         aeroback2lay=.false.     ! default value
436         call getin("aeroback2lay",aeroback2lay)
437         write(*,*)" aeroback2lay = ",aeroback2lay
438
439
440         write(*,*)"Cloud pressure level (with kastprof only):"
441         cloudlvl=0. ! default value
442         call getin("cloudlvl",cloudlvl)
443         write(*,*)" cloudlvl = ",cloudlvl
444
445         write(*,*)"Is the variable gas species radiatively active?"
446         Tstrat=167.0
447         varactive=.false.
448         call getin("varactive",varactive)
449         write(*,*)" varactive = ",varactive
450
451         write(*,*)"Is the variable gas species distribution set?"
452         varfixed=.false.
453         call getin("varfixed",varfixed)
454         write(*,*)" varfixed = ",varfixed
455
456         write(*,*)"What is the saturation % of the variable species?"
457         satval=0.8
458         call getin("satval",satval)
459         write(*,*)" satval = ",satval
460
461
462! Test of incompatibility:
463! if varactive, then varfixed should be false
464         if (varactive.and.varfixed) then
465           print*,'if varactive, varfixed must be OFF!'
466           stop
467         endif
468
469         write(*,*) "Gravitationnal sedimentation ?"
470         sedimentation=.false. ! default value
471         call getin("sedimentation",sedimentation)
472         write(*,*) " sedimentation = ",sedimentation
473
474         write(*,*) "Compute water cycle ?"
475         water=.false. ! default value
476         call getin("water",water)
477         write(*,*) " water = ",water
478         
479! Test of incompatibility:
480! if water is true, there should be at least a tracer
481         if (water.and.(.not.tracer)) then
482           print*,'if water is ON, tracer must be ON too!'
483           stop
484         endif
485
486         write(*,*) "Include water condensation ?"
487         watercond=.false. ! default value
488         call getin("watercond",watercond)
489         write(*,*) " watercond = ",watercond
490
491! Test of incompatibility:
492! if watercond is used, then water should be used too
493         if (watercond.and.(.not.water)) then
494           print*,'if watercond is used, water should be used too'
495           stop
496         endif
497
498         write(*,*) "Include water precipitation ?"
499         waterrain=.false. ! default value
500         call getin("waterrain",waterrain)
501         write(*,*) " waterrain = ",waterrain
502
503         write(*,*) "Include surface hydrology ?"
504         hydrology=.false. ! default value
505         call getin("hydrology",hydrology)
506         write(*,*) " hydrology = ",hydrology
507
508         write(*,*) "Evolve surface water sources ?"
509         sourceevol=.false. ! default value
510         call getin("sourceevol",sourceevol)
511         write(*,*) " sourceevol = ",sourceevol
512
513         write(*,*) "Ice evolution timestep ?"
514         icetstep=100.0 ! default value
515         call getin("icetstep",icetstep)
516         write(*,*) " icetstep = ",icetstep
517
518         write(*,*) "Snow albedo ?"
519         albedosnow=0.5         ! default value
520         call getin("albedosnow",albedosnow)
521         write(*,*) " albedosnow = ",albedosnow
522
523         write(*,*) "Maximum ice thickness ?"
524         maxicethick=2.0         ! default value
525         call getin("maxicethick",maxicethick)
526         write(*,*) " maxicethick = ",maxicethick
527
528         write(*,*) "Freezing point of seawater ?"
529         Tsaldiff=-1.8          ! default value
530         call getin("Tsaldiff",Tsaldiff)
531         write(*,*) " Tsaldiff = ",Tsaldiff
532
533         write(*,*) "Does user want to force cpp and mugaz?"
534         force_cpp=.false. ! default value
535         call getin("force_cpp",force_cpp)
536         write(*,*) " force_cpp = ",force_cpp
537
538         if (force_cpp) then
539           mugaz = -99999.
540           PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
541           call getin("mugaz",mugaz)
542           IF (mugaz.eq.-99999.) THEN
543               PRINT *, "mugaz must be set if force_cpp = T"
544               STOP
545           ELSE
546               write(*,*) "mugaz=",mugaz
547           ENDIF
548           cpp = -99999.
549           PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
550           call getin("cpp",cpp)
551           IF (cpp.eq.-99999.) THEN
552               PRINT *, "cpp must be set if force_cpp = T"
553               STOP
554           ELSE
555               write(*,*) "cpp=",cpp
556           ENDIF
557!         else
558!           mugaz=8.314*1000./pr
559         endif
560         call su_gases
561         call calc_cpp_mugaz
562
563         PRINT*,'--------------------------------------------'
564         PRINT*
565         PRINT*
566      ELSE
567         write(*,*)
568         write(*,*) 'Cannot read file callphys.def. Is it here ?'
569         stop
570      ENDIF
571
5728000  FORMAT(t5,a12,l8)
5738001  FORMAT(t5,a12,i8)
574
575      PRINT*
576      PRINT*,'inifis: daysec',daysec
577      PRINT*
578      PRINT*,'inifis: The radiative transfer is computed:'
579      PRINT*,'           each ',iradia,' physical time-step'
580      PRINT*,'        or each ',iradia*dtphys,' seconds'
581      PRINT*
582
583
584!-----------------------------------------------------------------------
585!     Some more initialization:
586!     ------------------------
587
588      ! ALLOCATE ARRAYS IN comgeomfi_h
589      IF (.not. ALLOCATED(lati)) ALLOCATE(lati(ngrid))
590      IF (.not. ALLOCATED(long)) ALLOCATE(long(ngrid))
591      IF (.not. ALLOCATED(area)) ALLOCATE(area(ngrid))
592
593      CALL SCOPY(ngrid,plon,1,long,1)
594      CALL SCOPY(ngrid,plat,1,lati,1)
595      CALL SCOPY(ngrid,parea,1,area,1)
596      totarea=SSUM(ngrid,area,1)
597
598      !! those are defined in comdiurn_h.F90
599      IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
600      IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
601      IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
602      IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
603
604      DO ig=1,ngrid
605         sinlat(ig)=sin(plat(ig))
606         coslat(ig)=cos(plat(ig))
607         sinlon(ig)=sin(plon(ig))
608         coslon(ig)=cos(plon(ig))
609      ENDDO
610
611      pi=2.*asin(1.) ! NB: pi is a common in comcstfi.h
612
613      RETURN
614      END
Note: See TracBrowser for help on using the repository browser.