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

Last change on this file since 486 was 486, checked in by rwordsworth, 13 years ago

Variable ice timestep added for ice evolution algorithm.
Sedimentation improved to allow consistent dust (L. Kerber).
Bug linked to allocatable matrices removed from rcm1d.F.
Treatment of initial aerosol radii in callcorrk.F90 improved.

File size: 17.1 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
9!=======================================================================
10!
11!   purpose:
12!   -------
13!
14!   Initialisation for the physical parametrisations of the LMD
15!   Generic Model.
16!
17!   author: Frederic Hourdin 15 / 10 /93
18!   -------
19!   modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
20!             Ehouarn Millour (oct. 2008) tracers are now identified
21!              by their names and may not be contiguously
22!              stored in the q(:,:,:,:) array
23!             E.M. (june 2009) use getin routine to load parameters
24!
25!
26!   arguments:
27!   ----------
28!
29!   input:
30!   ------
31!
32!    ngrid                 Size of the horizontal grid.
33!                          All internal loops are performed on that grid.
34!    nlayer                Number of vertical layers.
35!    pdayref               Day of reference for the simulation
36!    pday                  Number of days counted from the North. Spring
37!                          equinoxe.
38!
39!=======================================================================
40!
41!-----------------------------------------------------------------------
42!   declarations:
43!   -------------
44      use datafile_mod, only: datadir
45! to use  'getin'
46      USE ioipsl_getincom
47      IMPLICIT NONE
48#include "dimensions.h"
49#include "dimphys.h"
50#include "planete.h"
51#include "comcstfi.h"
52#include "comsaison.h"
53#include "comdiurn.h"
54#include "comgeomfi.h"
55#include "callkeys.h"
56#include "surfdat.h"
57
58
59
60      REAL prad,pg,pr,pcpp,pdaysec,ptimestep
61 
62      INTEGER ngrid,nlayer
63      REAL plat(ngrid),plon(ngrid),parea(ngridmx)
64      integer day_ini
65      INTEGER ig,ierr
66 
67      EXTERNAL iniorbit,orbite
68      EXTERNAL SSUM
69      REAL SSUM
70 
71      CHARACTER ch1*12
72      CHARACTER ch80*80
73
74      logical chem, h2o
75      logical :: parameter, doubleq=.false.
76
77      real psurf,pN2 ! added by RW for Gliese 581d N2+CO2
78
79      rad=prad
80      daysec=pdaysec
81      dtphys=ptimestep
82      cpp=pcpp
83      g=pg
84      r=pr
85      rcp=r/cpp
86
87      avocado = 6.02214179e23   ! added by RW
88
89! --------------------------------------------------------
90!     The usual Tests
91!     --------------
92      IF (nlayer.NE.nlayermx) THEN
93         PRINT*,'STOP in inifis'
94         PRINT*,'Probleme de dimensions :'
95         PRINT*,'nlayer     = ',nlayer
96         PRINT*,'nlayermx   = ',nlayermx
97         STOP
98      ENDIF
99
100      IF (ngrid.NE.ngridmx) THEN
101         PRINT*,'STOP in inifis'
102         PRINT*,'Probleme de dimensions :'
103         PRINT*,'ngrid     = ',ngrid
104         PRINT*,'ngridmx   = ',ngridmx
105         STOP
106      ENDIF
107
108! --------------------------------------------------------------
109!  Reading the "callphys.def" file controlling some key options
110! --------------------------------------------------------------
111     
112      ! check that 'callphys.def' file is around
113      OPEN(99,file='callphys.def',status='old',form='formatted'
114     &     ,iostat=ierr)
115      CLOSE(99)
116     
117      IF(ierr.EQ.0) THEN
118         PRINT*
119         PRINT*
120         PRINT*,'--------------------------------------------'
121         PRINT*,' inifis: Parametres pour la physique (callphys.def)'
122         PRINT*,'--------------------------------------------'
123
124         write(*,*) "Directory where external input files are:"
125         ! default 'datadir' is set in "datadir_mod"
126         call getin("datadir",datadir) ! default path
127         write(*,*) " datadir = ",trim(datadir)
128
129         write(*,*) "Run with or without tracer transport ?"
130         tracer=.false. ! default value
131         call getin("tracer",tracer)
132         write(*,*) " tracer = ",tracer
133
134         write(*,*) "Diurnal cycle ?"
135         write(*,*) "(if diurnal=false, diurnal averaged solar heating)"
136         diurnal=.true. ! default value
137         call getin("diurnal",diurnal)
138         write(*,*) " diurnal = ",diurnal
139
140         write(*,*) "Seasonal cycle ?"
141         write(*,*) "(if season=false, Ls stays constant, to value ",
142     &   "set in 'start'"
143         season=.true. ! default value
144         call getin("season",season)
145         write(*,*) " season = ",season
146
147         write(*,*) "Tidally resonant rotation ?"
148         tlocked=.false. ! default value
149         call getin("tlocked",tlocked)
150         write(*,*) "tlocked = ",tlocked
151
152! Test of incompatibility:
153! if tlocked, then diurnal should be false
154         if (tlocked.and.diurnal) then
155           print*,'If diurnal=true, we should turn off tlocked.'
156           stop
157         endif
158
159         write(*,*) "Tidal resonance ratio ?"
160         nres=0          ! default value
161         call getin("nres",nres)
162         write(*,*) "nres = ",nres
163
164         write(*,*) "Write some extra output to the screen ?"
165         lwrite=.false. ! default value
166         call getin("lwrite",lwrite)
167         write(*,*) " lwrite = ",lwrite
168
169         write(*,*) "Save statistics in file stats.nc ?"
170         callstats=.true. ! default value
171         call getin("callstats",callstats)
172         write(*,*) " callstats = ",callstats
173
174         write(*,*) "Test energy conservation of model physics ?"
175         enertest=.false. ! default value
176         call getin("enertest",enertest)
177         write(*,*) " enertest = ",enertest
178
179         write(*,*) "call radiative transfer ?"
180         callrad=.true. ! default value
181         call getin("callrad",callrad)
182         write(*,*) " callrad = ",callrad
183
184         write(*,*) "call correlated-k radiative transfer ?"
185         corrk=.true. ! default value
186         call getin("corrk",corrk)
187         write(*,*) " corrk = ",corrk
188
189         write(*,*) "call gaseous absorption in the visible bands?",
190     &              "(matters only if callrad=T)"
191         callgasvis=.false. ! default value
192         call getin("callgasvis",callgasvis)
193         write(*,*) " callgasvis = ",callgasvis
194         
195         write(*,*) "call turbulent vertical diffusion ?"
196         calldifv=.true. ! default value
197         call getin("calldifv",calldifv)
198         write(*,*) " calldifv = ",calldifv
199
200         write(*,*) "call convective adjustment ?"
201         calladj=.true. ! default value
202         call getin("calladj",calladj)
203         write(*,*) " calladj = ",calladj
204       
205         write(*,*)"Gas is nonideal CO2 ?"
206         nonideal=.false.
207         call getin("nonideal",nonideal)
208         write(*,*)" nonideal = ",nonideal
209
210! Test of incompatibility
211         if (enertest.and.nonideal) then
212            print*,'Energy conservation calculations currently
213     &           assume ideal gas!'
214            call abort
215         endif
216
217         write(*,*) "call CO2 condensation ?"
218         co2cond=.true. ! default value
219         call getin("co2cond",co2cond)
220         write(*,*) " co2cond = ",co2cond
221
222         write(*,*) "CO2 supersaturation level ?"
223         co2supsat=1.0 ! default value
224         call getin("co2supsat",co2supsat)
225         write(*,*) " co2supsat = ",co2supsat
226
227         write(*,*) "Radiative timescale for Newtonian cooling ?"
228         tau_relax=30. ! default value
229         call getin("tau_relax",tau_relax)
230         write(*,*) " tau_relax = ",tau_relax
231         tau_relax=tau_relax*24*3600 ! convert Earth days --> seconds
232
233         write(*,*)"call thermal conduction in the soil ?"
234         callsoil=.true. ! default value
235         call getin("callsoil",callsoil)
236         write(*,*) " callsoil = ",callsoil
237         
238         write(*,*)"Rad transfer is computed every iradia",
239     &             " physical timestep"
240         iradia=1 ! default value
241         call getin("iradia",iradia)
242         write(*,*)" iradia = ",iradia
243       
244         write(*,*)"Rayleigh scattering ?"
245         rayleigh=.false.
246         call getin("rayleigh",rayleigh)
247         write(*,*)" rayleigh = ",rayleigh
248
249         write(*,*) "Use blackbody for stellar spectrum ?"
250         stelbbody=.false. ! default value
251         call getin("stelbbody",stelbbody)
252         write(*,*) " stelbbody = ",stelbbody
253
254         write(*,*) "Stellar blackbody temperature ?"
255         stelTbb=5800.0 ! default value
256         call getin("stelTbb",stelTbb)
257         write(*,*) " stelTbb = ",stelTbb
258
259         write(*,*)"Output mean OLR in 1D?"
260         meanOLR=.false.
261         call getin("meanOLR",meanOLR)
262         write(*,*)" meanOLR = ",meanOLR
263
264         write(*,*)"Output spectral OLR in 3D?"
265         specOLR=.false.
266         call getin("specOLR",specOLR)
267         write(*,*)" specOLR = ",specOLR
268
269         write(*,*)"Operate in kastprof mode?"
270         kastprof=.false.
271         call getin("kastprof",kastprof)
272         write(*,*)" kastprof = ",kastprof
273
274         write(*,*)"Uniform absorption coefficient in IR?"
275         graybody=.false.
276         call getin("graybody",graybody)
277         write(*,*)" graybody = ",graybody
278
279! Test of incompatibility:
280! if kastprof used, we must be in 1D
281         if (kastprof.and.(ngridmx.gt.1)) then
282           print*,'kastprof can only be used in 1D!'
283           call abort
284         endif
285
286         write(*,*)"Stratospheric temperature for kastprof mode?"
287         Tstrat=167.0
288         call getin("Tstrat",Tstrat)
289         write(*,*)" Tstrat = ",Tstrat
290
291         write(*,*)"Remove lower boundary?"
292         noradsurf=.false.
293         call getin("noradsurf",noradsurf)
294         write(*,*)" noradsurf = ",noradsurf
295
296! Tests of incompatibility:
297         if (noradsurf.and.callsoil) then
298           print*,'noradsurf not compatible with soil scheme!'
299           call abort
300         endif
301         if (noradsurf.and.calldifv) then
302           print*,'noradsurf not compatible with a boundary layer!'
303           call abort
304         endif
305
306
307         write(*,*)"Use Newtonian cooling for radiative transfer?"
308         newtonian=.false.
309         call getin("newtonian",newtonian)
310         write(*,*)" newtonian = ",newtonian
311
312! Tests of incompatibility:
313         if (newtonian.and.corrk) then
314           print*,'newtonian not compatible with correlated-k!'
315           call abort
316         endif
317         if (newtonian.and.calladj) then
318           print*,'newtonian not compatible with adjustment!'
319           call abort
320         endif
321         if (newtonian.and.calldifv) then
322           print*,'newtonian not compatible with a boundary layer!'
323           call abort
324         endif
325
326         write(*,*)"Test physics timescale in 1D?"
327         testradtimes=.false.
328         call getin("testradtimes",testradtimes)
329         write(*,*)" testradtimes = ",testradtimes
330
331! Test of incompatibility:
332! if testradtimes used, we must be in 1D
333         if (testradtimes.and.(ngridmx.gt.1)) then
334           print*,'testradtimes can only be used in 1D!'
335           call abort
336         endif
337
338         write(*,*)"Default planetary temperature?"
339         tplanet=215.0
340         call getin("tplanet",tplanet)
341         write(*,*)" tplanet = ",tplanet
342
343         write(*,*)"Which star?"
344         startype=1 ! default value = Sol
345         call getin("startype",startype)
346         write(*,*)" startype = ",startype
347
348         write(*,*)"Value of stellar flux at 1 AU?"
349         Fat1AU=1356.0 ! default value = Sol today
350         call getin("Fat1AU",Fat1AU)
351         write(*,*)" Fat1AU = ",Fat1AU
352
353
354!     1D solar zenith angle
355         write(*,*)"Solar zenith angle in 1D?"
356         szangle=60.0
357         call getin("szangle",szangle)
358         write(*,*)" szangle = ",szangle
359
360! TRACERS:
361
362         write(*,*)"Fixed / inactive aerosol distribution?"
363         aerofixed=.true.       ! default value
364         call getin("aerofixed",aerofixed)
365         write(*,*)" aerofixed = ",aerofixed
366
367         write(*,*)"Varying H2O cloud fraction?"
368         CLFvarying=.false.     ! default value
369         call getin("CLFvarying",CLFvarying)
370         write(*,*)" CLFvarying = ",CLFvarying
371
372         write(*,*)"Value of fixed H2O cloud fraction?"
373         CLFfixval=1.0                ! default value
374         call getin("CLFfixval",CLFfixval)
375         write(*,*)" CLFfixval = ",CLFfixval
376
377         write(*,*)"Number mixing ratio of CO2 ice particles:"
378         Nmix_co2=100000. ! default value
379         call getin("Nmix_co2",Nmix_co2)
380         write(*,*)" Nmix_co2 = ",Nmix_co2
381
382         write(*,*)"Number mixing ratio of H2O ice particles:"
383         Nmix_h2o=10000000. ! default value
384         call getin("Nmix_h2o",Nmix_h2o)
385         write(*,*)" Nmix_h2o = ",Nmix_h2o
386
387         write(*,*)"Opacity of dust (if used):"
388         dusttau=0. ! default value
389         call getin("dusttau",dusttau)
390         write(*,*)" dusttau = ",dusttau
391
392! Test of incompatibility:
393! if less than 3 aerosols, then dusttau should = 0
394         if ((naerkind.lt.3).and.dusttau.gt.0.) then
395           print*,'Need naer>2 for radiatively active dust!'
396           stop
397         endif
398
399         write(*,*)"Cloud pressure level (with kastprof only):"
400         cloudlvl=0. ! default value
401         call getin("cloudlvl",cloudlvl)
402         write(*,*)" cloudlvl = ",cloudlvl
403
404         write(*,*)"Is the variable gas species radiatively active?"
405         Tstrat=167.0
406         varactive=.false.
407         call getin("varactive",varactive)
408         write(*,*)" varactive = ",varactive
409
410         write(*,*)"Is the variable gas species distribution set?"
411         varfixed=.false.
412         call getin("varfixed",varfixed)
413         write(*,*)" varfixed = ",varfixed
414
415         write(*,*)"What is the saturation % of the variable species?"
416         satval=0.8
417         call getin("satval",satval)
418         write(*,*)" satval = ",satval
419
420! Test of incompatibility:
421! if no tracers, then aerofixed should be true
422         if ((.not.tracer).and.(.not.aerofixed)) then
423           print*,'if tracers are off, aerofixed must be ON!'
424           stop
425         endif
426
427! Test of incompatibility:
428! if varactive, then varfixed should be false
429         if (varactive.and.varfixed) then
430           print*,'if varactive, varfixed must be OFF!'
431           stop
432         endif
433
434! Test of incompatibility:
435
436         write(*,*) "Gravitationnal sedimentation ?"
437         sedimentation=.true. ! default value
438         call getin("sedimentation",sedimentation)
439         write(*,*) " sedimentation = ",sedimentation
440
441
442! Test of incompatibility:
443
444         write(*,*) "Compute water cycle ?"
445         water=.false. ! default value
446         call getin("water",water)
447         write(*,*) " water = ",water
448         
449         write(*,*) "Include water condensation ?"
450         watercond=.false. ! default value
451         call getin("watercond",watercond)
452         write(*,*) " watercond = ",watercond
453
454         write(*,*) "Include water precipitation ?"
455         waterrain=.false. ! default value
456         call getin("waterrain",waterrain)
457         write(*,*) " waterrain = ",waterrain
458
459         write(*,*) "Precipitation threshold ?"
460         rainthreshold=0.011 ! default value (Emmanuel 1997)
461         call getin("rainthreshold",rainthreshold)
462         write(*,*) " rainthreshold = ",rainthreshold
463
464         write(*,*) "Include surface hydrology ?"
465         hydrology=.false. ! default value
466         call getin("hydrology",hydrology)
467         write(*,*) " hydrology = ",hydrology
468
469         write(*,*) "Evolve surface water sources ?"
470         sourceevol=.false. ! default value
471         call getin("sourceevol",sourceevol)
472         write(*,*) " sourceevol = ",sourceevol
473
474         write(*,*) "Ice evolution timestep ?"
475         icetstep=100.0 ! default value
476         call getin("icetstep",icetstep)
477         write(*,*) " icetstep = ",icetstep
478
479         write(*,*) "Snow albedo ?"
480         albedosnow=0.5         ! default value
481         call getin("albedosnow",albedosnow)
482         write(*,*) " albedosnow = ",albedosnow
483
484         write(*,*) "Maximum ice thickness ?"
485         maxicethick=2.0         ! default value
486         call getin("maxicethick",maxicethick)
487         write(*,*) " maxicethick = ",maxicethick
488
489         write(*,*) "Freezing point of seawater ?"
490         Tsaldiff=-1.8          ! default value
491         call getin("Tsaldiff",Tsaldiff)
492         write(*,*) " Tsaldiff = ",Tsaldiff
493
494! Test of incompatibility:
495! if watercond is used, then water should be used too
496
497         if (watercond.and.(.not.watercond)) then
498           print*,'if watercond is used, water should be used too'
499           stop
500         endif
501
502         mugaz=0.
503         cpp=0.
504         call su_gases
505         call calc_cpp_mugaz
506
507         PRINT*,'--------------------------------------------'
508         PRINT*
509         PRINT*
510      ELSE
511         write(*,*)
512         write(*,*) 'Cannot read file callphys.def. Is it here ?'
513         stop
514      ENDIF
515
5168000  FORMAT(t5,a12,l8)
5178001  FORMAT(t5,a12,i8)
518
519      PRINT*
520      PRINT*,'inifis: daysec',daysec
521      PRINT*
522      PRINT*,'inifis: The radiative transfer is computed:'
523      PRINT*,'           each ',iradia,' physical time-step'
524      PRINT*,'        or each ',iradia*dtphys,' seconds'
525      PRINT*
526
527
528!-----------------------------------------------------------------------
529!     Some more initialization:
530!     ------------------------
531
532      CALL SCOPY(ngrid,plon,1,long,1)
533      CALL SCOPY(ngrid,plat,1,lati,1)
534      CALL SCOPY(ngrid,parea,1,area,1)
535      totarea=SSUM(ngridmx,area,1)
536
537      DO ig=1,ngrid
538         sinlat(ig)=sin(plat(ig))
539         coslat(ig)=cos(plat(ig))
540         sinlon(ig)=sin(plon(ig))
541         coslon(ig)=cos(plon(ig))
542      ENDDO
543
544      pi=2.*asin(1.) ! NB: pi is a common in comcstfi.h
545
546      RETURN
547      END
Note: See TracBrowser for help on using the repository browser.