source: trunk/LMDZ.MARS/libf/aeronomars/inichim_readcallphys.F @ 503

Last change on this file since 503 was 321, checked in by emillour, 13 years ago

Mars GCM:

  • Corrected small bug in newstart: initracer was not always used and thus some tracer indexes (igm_co2, igcm_h2o_vap,...) were not set. This however means that we now also call inifis from newstart and that we read in flags set in 'callphys.def' (required for sanity checks in initracer). Also adapted 'inichim_readcallphys': removed some obsolescent tests on number of tracers for given combinations of options.

EM

File size: 16.6 KB
Line 
1      SUBROUTINE inichim_readcallphys
2!
3!=======================================================================
4!
5!   subject:
6!   --------
7!
8!   Initialisation for the physical parametrisations of the LMD
9!   martian atmospheric general circulation model.
10!   This routine is called by inichim_newstart (which is used by programs
11!   newstart and testphys1d but not the GCM)
12!
13!   author: Frederic Hourdin 15 / 10 /93
14!   -------
15!   modifications: Sebastien Lebonnois 11/06/2003 (new callphys.def)
16!                  10/2008 adapted to using tracer sby name. Ehouarn Millour
17!                  07/2009 use 'getin' to read callphys.def. EM
18!
19!   arguments:
20!   ----------
21!
22!   input:
23!   ------
24!
25!=======================================================================
26!
27!-----------------------------------------------------------------------
28!   declarations:
29!   -------------
30! to use  'getin'
31      USE ioipsl_getincom
32      IMPLICIT NONE
33#include "dimensions.h"
34#include "dimphys.h"
35#include "planete.h"
36#include "comcstfi.h"
37#include "comsaison.h"
38#include "comdiurn.h"
39#include "comgeomfi.h"
40#include "callkeys.h"
41#include "surfdat.h"
42
43      character*12 ch1
44      integer ierr
45      logical chem, h2o
46
47
48! --------------------------------------------------------------
49!  Reading the "callphys.def" file controlling some key options
50! --------------------------------------------------------------
51     
52      ! check that 'callphys.def' file is around
53      OPEN(99,file='callphys.def',status='old',form='formatted'
54     &     ,iostat=ierr)
55      CLOSE(99)
56     
57      IF(ierr.EQ.0) THEN
58         PRINT*
59         PRINT*
60         PRINT*,'--------------------------------------------'
61         PRINT*,' inichim_readcallphys: Parametres pour la physique',
62     &          ' (callphys.def)'
63         PRINT*,'--------------------------------------------'
64
65
66         write(*,*) "Run with or without tracer transport ?"
67         tracer=.false. ! default value
68         call getin("tracer",tracer)
69         write(*,*) " tracer = ",tracer
70
71         write(*,*) "Diurnal cycle ?"
72         write(*,*) "(if diurnal=False, diurnal averaged solar heating)"
73         diurnal=.true. ! default value
74         call getin("diurnal",diurnal)
75         write(*,*) " diurnal = ",diurnal
76
77         write(*,*) "Seasonal cycle ?"
78         write(*,*) "(if season=False, Ls stays constant, to value ",
79     &   "set in 'start'"
80         season=.true. ! default value
81         call getin("season",season)
82         write(*,*) " season = ",season
83
84         write(*,*) "Write some extra output to the screen ?"
85         lwrite=.false. ! default value
86         call getin("lwrite",lwrite)
87         write(*,*) " lwrite = ",lwrite
88
89         write(*,*) "Save statistics in file stats.nc ?"
90         callstats=.true. ! default value
91         call getin("callstats",callstats)
92         write(*,*) " callstats = ",callstats
93
94         write(*,*) "Save EOF profiles in file 'profiles' for ",
95     &              "Climate Database?"
96         calleofdump=.false. ! default value
97         call getin("calleofdump",calleofdump)
98         write(*,*) " calleofdump = ",calleofdump
99
100         write(*,*) "Dust scenario:"
101         iaervar=3 ! default value
102         call getin("iaervar",iaervar)
103         write(*,*) " iaervar = ",iaervar
104
105         write(*,*) "Dust vertical distribution:"
106         write(*,*) "(=1 Dust opt.deph read in startfi;",
107     & " =2 Viking scenario; =3 MGS scenario,",
108     & " =4 Mars Year 24 from TES assimilation)"
109         iddist=3 ! default value
110         call getin("iddist",iddist)
111         write(*,*) " iddist = ",iddist
112
113         write(*,*) "Dust top altitude (km). (Matters only if iddist=1)"
114         topdustref= 90.0 ! default value
115         call getin("topdustref",topdustref)
116         write(*,*) " topdustref = ",topdustref
117
118
119         write(*,*) "call radiative transfer ?"
120         callrad=.true. ! default value
121         call getin("callrad",callrad)
122         write(*,*) " callrad = ",callrad
123
124         write(*,*) "call NLTE radiative schemes ?",
125     &              "(matters only if callrad=T)"
126         callnlte=.false. ! default value
127         call getin("callnlte",callnlte)
128         write(*,*) " callnlte = ",callnlte
129         
130         write(*,*) "call CO2 NIR absorption ?",
131     &              "(matters only if callrad=T)"
132         callnirco2=.false. ! default value
133         call getin("callnirco2",callnirco2)
134         write(*,*) " callnirco2 = ",callnirco2
135         
136         write(*,*) "call turbulent vertical diffusion ?"
137         calldifv=.true. ! default value
138         call getin("calldifv",calldifv)
139         write(*,*) " calldifv = ",calldifv
140
141         write(*,*) "call convective adjustment ?"
142         calladj=.true. ! default value
143         call getin("calladj",calladj)
144         write(*,*) " calladj = ",calladj
145         
146
147         write(*,*) "call CO2 condensation ?"
148         callcond=.true. ! default value
149         call getin("callcond",callcond)
150         write(*,*) " callcond = ",callcond
151
152         write(*,*)"call thermal conduction in the soil ?"
153         callsoil=.true. ! default value
154         call getin("callsoil",callsoil)
155         write(*,*) " callsoil = ",callsoil
156         
157
158         write(*,*)"call Lott's gravity wave/subgrid topography ",
159     &             "scheme ?"
160         calllott=.true. ! default value
161         call getin("calllott",calllott)
162         write(*,*)" calllott = ",calllott
163
164
165         write(*,*)"rad.transfer is computed every iradia",
166     &             " physical timestep"
167         iradia=1 ! default value
168         call getin("iradia",iradia)
169         write(*,*)" iradia = ",iradia
170         
171
172         write(*,*)"Output of the exchange coefficient mattrix ?",
173     &             "(for diagnostics only)"
174         callg2d=.false. ! default value
175         call getin("callg2d",callg2d)
176         write(*,*)" callg2d = ",callg2d
177
178         write(*,*)"Rayleigh scattering : (should be .false. for now)"
179         rayleigh=.false.
180         call getin("rayleigh",rayleigh)
181         write(*,*)" rayleigh = ",rayleigh
182
183
184! TRACERS:
185
186         write(*,*)"Transported dust ? (if >0, use 'dustbin' dust bins)"
187         dustbin=0 ! default value
188         call getin("dustbin",dustbin)
189         write(*,*)" dustbin = ",dustbin
190
191         write(*,*)"Radiatively active dust ? (matters if dustbin>0)"
192         active=.false. ! default value
193         call getin("active",active)
194         write(*,*)" active = ",active
195
196! Test of incompatibility:
197! if active is used, then dustbin should be > 0
198
199         if (active.and.(dustbin.lt.1)) then
200           print*,'if active is used, then dustbin should > 0'
201           stop
202         endif
203
204         write(*,*)"use mass and number mixing ratios to predict",
205     &             " dust size ?"
206         doubleq=.false. ! default value
207         call getin("doubleq",doubleq)
208         write(*,*)" doubleq = ",doubleq
209
210! Test of incompatibility:
211! if doubleq is used, then dustbin should be at least 2
212
213         if (doubleq.and.(dustbin.lt.2)) then
214           print*,'if doubleq is used, then dustbin should be > 1'
215           stop
216         endif
217
218         write(*,*)"dust lifted by GCM surface winds ?"
219         lifting=.false. ! default value
220         call getin("lifting",lifting)
221         write(*,*)" lifting = ",lifting
222
223! Test of incompatibility:
224! if lifting is used, then dustbin should be > 0
225
226         if (lifting.and.(dustbin.lt.1)) then
227           print*,'if lifting is used, then dustbin should > 0'
228           stop
229         endif
230
231         write(*,*)" dust lifted by dust devils ?"
232         callddevil=.false. !default value
233         call getin("callddevil",callddevil)
234         write(*,*)" callddevil = ",callddevil
235         
236
237! Test of incompatibility:
238! if dustdevil is used, then dustbin should be > 0
239
240         if (callddevil.and.(dustbin.lt.1)) then
241           print*,'if dustdevil is used, then dustbin should > 0'
242           stop
243         endif
244
245         write(*,*)"Dust scavenging by CO2 snowfall ?"
246         scavenging=.false. ! default value
247         call getin("scavenging",scavenging)
248         write(*,*)" scavenging = ",scavenging
249         
250
251! Test of incompatibility:
252! if scavenging is used, then dustbin should be > 0
253
254         if (scavenging.and.(dustbin.lt.1)) then
255           print*,'if scavenging is used, then dustbin should > 0'
256           stop
257         endif
258
259         write(*,*) "Gravitationnal sedimentation ?"
260         sedimentation=.true. ! default value
261         call getin("sedimentation",sedimentation)
262         write(*,*) " sedimentation = ",sedimentation
263
264         write(*,*) "includes water ice",
265     &              "(if true, 'water' must also be .true.)"
266
267         write(*,*) "Radiatively active transported atmospheric ",
268     &              "water ice ?"
269         activice=.false. ! default value
270         call getin("activice",activice)
271         write(*,*) " activice = ",activice
272
273         write(*,*) "Compute water cycle ?"
274         water=.false. ! default value
275         call getin("water",water)
276         write(*,*) " water = ",water
277         write(*,*) "Permanent water caps at poles ?",
278     &               " .true. is RECOMMENDED"
279         write(*,*) "(with .true., North cap is a source of water ",
280     &   "and South pole is a cold trap)"
281         caps=.true. ! default value
282         call getin("caps",caps)
283         write(*,*) " caps = ",caps
284
285! Test of incompatibility:
286! if activice is used, then water should be used too
287
288         if (activice.and..not.water) then
289           print*,'if activice is used, water should be used too'
290           stop
291         endif
292
293         write(*,*) "photochemistry: include chemical species"
294         photochem=.false. ! default value
295         call getin("photochem",photochem)
296         write(*,*) " photochem = ",photochem
297
298
299! THERMOSPHERE
300
301         write(*,*) "call thermosphere ?"
302         callthermos=.false. ! default value
303         call getin("callthermos",callthermos)
304         write(*,*) " callthermos = ",callthermos
305         
306
307         write(*,*) " water included without cycle ",
308     &              "(only if water=.false.)"
309         thermoswater=.false. ! default value
310         call getin("thermoswater",thermoswater)
311         write(*,*) " thermoswater = ",thermoswater
312
313         write(*,*) "call thermal conduction ?",
314     &    " (only if callthermos=.true.)"
315         callconduct=.false. ! default value
316         call getin("callconduct",callconduct)
317         write(*,*) " callconduct = ",callconduct
318
319         write(*,*) "call EUV heating ?",
320     &   " (only if callthermos=.true.)"
321         calleuv=.false.  ! default value
322         call getin("calleuv",calleuv)
323         write(*,*) " calleuv = ",calleuv
324
325         write(*,*) "call molecular viscosity ?",
326     &   " (only if callthermos=.true.)"
327         callmolvis=.false. ! default value
328         call getin("callmolvis",callmolvis)
329         write(*,*) " callmolvis = ",callmolvis
330
331         write(*,*) "call molecular diffusion ?",
332     &   " (only if callthermos=.true.)"
333         callmoldiff=.false. ! default value
334         call getin("callmoldiff",callmoldiff)
335         write(*,*) " callmoldiff = ",callmoldiff
336         
337
338         write(*,*) "call thermospheric photochemistry ?",
339     &   " (only if callthermos=.true.)"
340         thermochem=.false. ! default value
341         call getin("thermochem",thermochem)
342         write(*,*) " thermochem = ",thermochem
343
344         write(*,*) "date for solar flux calculation:",
345     &   " (1985 < date < 2002)"
346         write(*,*) "(Solar min=1996.4 ave=1993.4 max=1990.6)"
347         solarcondate=1993.4 ! default value
348         call getin("solarcondate",solarcondate)
349         write(*,*) " solarcondate = ",solarcondate
350         
351
352         if (.not.callthermos) then
353           if (thermoswater) then
354             print*,'if thermoswater is set, callthermos must be true'
355             stop
356           endif         
357           if (callconduct) then
358             print*,'if callconduct is set, callthermos must be true'
359             stop
360           endif       
361           if (calleuv) then
362             print*,'if calleuv is set, callthermos must be true'
363             stop
364           endif         
365           if (callmolvis) then
366             print*,'if callmolvis is set, callthermos must be true'
367             stop
368           endif       
369           if (callmoldiff) then
370             print*,'if callmoldiff is set, callthermos must be true'
371             stop
372           endif         
373           if (thermochem) then
374             print*,'if thermochem is set, callthermos must be true'
375             stop
376           endif         
377        endif
378
379! Test of incompatibility:
380! if photochem is used, then water should be used too
381
382         if (photochem.and..not.water) then
383           print*,'if photochem is used, water should be used too'
384           stop
385         endif
386
387! if callthermos is used, then thermoswater should be used too
388! (if water not used already)
389
390         if (callthermos .and. .not.water) then
391           if (callthermos .and. .not.thermoswater) then
392             print*,'if callthermos is used, water or thermoswater
393     &               should be used too'
394             stop
395           endif
396         endif
397
398         PRINT*,'--------------------------------------------'
399         PRINT*
400         PRINT*
401      ELSE
402         write(*,*)
403         write(*,*) 'Cannot read file callphys.def. Is it here ?'
404         stop
405      ENDIF
406      CLOSE(99)
407
408      pi=2.*asin(1.)
409
410!     managing the tracers, and tests:
411!     -------------------------------
412
413      if(tracer) then
414
415!          when photochem is used, nqchem_min is the rank
416!          of the first chemical species
417
418! Ehouarn: nqchem_min is now meaningless and no longer used
419!       nqchem_min = 1
420       if (photochem .or. callthermos) then
421         chem = .true.
422!        if (doubleq) then
423!          nqchem_min = 3
424!        else
425!          nqchem_min = dustbin+1
426!        end if
427       end if
428
429       if (water .or. thermoswater) h2o = .true.
430
431c          TESTS
432
433       print*,'inichim_readcallphys: TRACERS:'
434
435       if ((doubleq).and.(h2o).and.
436     $     (chem)) then
437!         print*,' 1: dust ; 2: dust (doubleq)'
438!         print*,' 3 to ',nqmx-2,': chemistry'
439!         print*,nqmx-1,': water ice ; ',nqmx,': water vapor'
440         print*,' 2 dust tracers (doubleq)'
441         print*,' 1 water vapour tracer'
442         print*,' 1 water ice tracer'
443         print*,nqmx-4,' chemistry tracers'
444       endif
445
446       if ((doubleq).and.(h2o).and.
447     $     .not.(chem)) then
448!         print*,' 1: dust ; 2: dust (doubleq)'
449!         print*,nqmx-1,': water ice ; ',nqmx,': water vapor'
450         print*,' 2 dust tracers (doubleq)'
451         print*,' 1 water vapour tracer'
452         print*,' 1 water ice tracer'
453         if (nqmx.lt.4) then
454           print*,'nqmx should be at least 4 with these options.'
455                   print*,'(or check callphys.def)'
456           stop
457         endif
458       endif
459
460!       if ((doubleq).and..not.(h2o)) then
461!         print*,' 1: dust ; 2: dust (doubleq)'
462!         print*,' 2 dust tracers (doubleq)'
463!         if (nqmx.ne.2) then
464!           print*,'nqmx should be 2 with these options...'
465!                  print*,'(or check callphys.def)'
466!           stop
467!         endif
468!       endif
469
470!       if (.not.(doubleq).and.(h2o).and.
471!     $     (chem)) then
472!         if (dustbin.gt.0) then
473!           print*,' 1 to ',dustbin,': dust bins'
474!           print*,dustbin,' dust bins'
475!         endif
476!         print*,nqchem_min,' to ',nqmx-2,': chemistry'
477!         print*,nqmx-1,': water ice ; ',nqmx,': water vapor'
478!         print*,nqmx-2-dustbin,' chemistry tracers'
479!         print*,' 1 water vapour tracer'
480!         print*,' 1 water ice tracer'
481!       endif
482
483!       if (.not.(doubleq).and.(h2o).and.
484!     $     .not.(chem)) then
485!         if (dustbin.gt.0) then
486!           print*,' 1 to ',dustbin,': dust bins'
487!           print*,dustbin,' dust bins'
488!         endif
489!         print*,nqmx-1,': water ice ; ',nqmx,': water vapor'
490!         print*,' 1 water vapour tracer'
491!         print*,' 1 water ice tracer'
492!         if (nqmx.ne.(dustbin+2)) then
493!           print*,'nqmx should be ',(dustbin+2),
494!     $            ' with these options...'
495!                  print*,'(or check callphys.def)'
496!           stop
497!         endif
498!       endif
499
500!       if (.not.(doubleq).and..not.(h2o)) then
501!         if (dustbin.gt.0) then
502!           print*,' 1 to ',dustbin,': dust bins'
503!           print*,dustbin,' dust bins'
504!           if (nqmx.ne.dustbin) then
505!             print*,'nqmx should be ',dustbin,
506!     $              ' with these options...'
507!             print*,'(or check callphys.def)'
508!             stop
509!           endif
510!         else
511!           print*,'dustbin=',dustbin,
512!     $            ': tracer should be F with these options...'
513!     $           ,'UNLESS you just want to move tracers around '
514!         endif
515!       endif
516
517      endif ! of if (tracer)
518
519      RETURN
520      END
Note: See TracBrowser for help on using the repository browser.