source: trunk/LMDZ.MARS/libf/phymars/inifis.F @ 420

Last change on this file since 420 was 420, checked in by tnavarro, 13 years ago

24/11/11 == TN

corrected minor bug in updatereffrad.F : reffdust was not saved

ccn_factor as not correctly used in sedimentation.

It is now initialized in inifis.F, declared in tracer.h and
used in both simpleclouds.F & callsedim.F to update ice radius.

Commented diagfi outputs in aeropacity.F & improvedclouds.F for non scavenging users.

File size: 25.3 KB
Line 
1      SUBROUTINE inifis(
2     $           ngrid,nlayer
3     $           ,day_ini,pdaysec,ptimestep
4     $           ,plat,plon,parea
5     $           ,prad,pg,pr,pcpp
6#ifdef MESOSCALE
7#include "meso_inc/meso_inc_inifisinvar.F"
8#endif
9     $           )
10!
11!=======================================================================
12!
13!   purpose:
14!   -------
15!
16!   Initialisation for the physical parametrisations of the LMD
17!   martian atmospheric general circulation modele.
18!
19!   author: Frederic Hourdin 15 / 10 /93
20!   -------
21!   modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
22!             Ehouarn Millour (oct. 2008) tracers are now identified
23!              by their names and may not be contiguously
24!              stored in the q(:,:,:,:) array
25!             E.M. (june 2009) use getin routine to load parameters
26!             adapted to the mesoscale use - Aymeric Spiga - 01/2007-07/2011
27!
28!
29!   arguments:
30!   ----------
31!
32!   input:
33!   ------
34!
35!    ngrid                 Size of the horizontal grid.
36!                          All internal loops are performed on that grid.
37!    nlayer                Number of vertical layers.
38!    pdayref               Day of reference for the simulation
39!    pday                  Number of days counted from the North. Spring
40!                          equinoxe.
41!
42!=======================================================================
43!
44!-----------------------------------------------------------------------
45!   declarations:
46!   -------------
47! to use  'getin'
48      USE ioipsl_getincom
49      IMPLICIT NONE
50#include "dimensions.h"
51#include "dimphys.h"
52#include "planete.h"
53#include "comcstfi.h"
54#include "comsaison.h"
55#include "comdiurn.h"
56#include "comgeomfi.h"
57#include "callkeys.h"
58#include "surfdat.h"
59#include "dimradmars.h"
60#include "yomaer.h"
61#include "datafile.h"
62#include "slope.h"
63#include "microphys.h"
64#include "tracer.h"
65#ifdef MESOSCALE
66#include "comsoil.h"     !!MESOSCALE -- needed to fill volcapa
67#include "meso_inc/meso_inc_inifisvar.F"
68#endif
69      REAL prad,pg,pr,pcpp,pdaysec
70
71      REAL ptimestep
72      INTEGER day_ini
73
74      INTEGER ngrid,nlayer
75      REAL plat(ngrid),plon(ngrid),parea(ngridmx)
76      INTEGER ig,ierr
77 
78!      EXTERNAL iniorbit,orbite
79      EXTERNAL SSUM
80      REAL SSUM
81 
82      CHARACTER ch1*12
83      CHARACTER ch80*80
84
85!      logical chem, h2o
86
87!      chem = .false.
88!      h2o = .false.
89
90      rad=prad
91      cpp=pcpp
92      g=pg
93      r=pr
94      rcp=r/cpp
95      daysec=pdaysec
96      dtphys=ptimestep
97#ifdef MESOSCALE
98#include "meso_inc/meso_inc_inifisini.F"
99#endif
100
101! --------------------------------------------------------
102!     The usual Tests
103!     --------------
104      IF (nlayer.NE.nlayermx) THEN
105         PRINT*,'STOP in inifis'
106         PRINT*,'Probleme de dimensions :'
107         PRINT*,'nlayer     = ',nlayer
108         PRINT*,'nlayermx   = ',nlayermx
109         STOP
110      ENDIF
111
112      IF (ngrid.NE.ngridmx) THEN
113         PRINT*,'STOP in inifis'
114         PRINT*,'Probleme de dimensions :'
115         PRINT*,'ngrid     = ',ngrid
116         PRINT*,'ngridmx   = ',ngridmx
117         STOP
118      ENDIF
119
120! --------------------------------------------------------------
121!  Reading the "callphys.def" file controlling some key options
122! --------------------------------------------------------------
123     
124      ! check that 'callphys.def' file is around
125      OPEN(99,file='callphys.def',status='old',form='formatted'
126     &     ,iostat=ierr)
127      CLOSE(99)
128     
129      IF(ierr.EQ.0) THEN
130         PRINT*
131         PRINT*
132         PRINT*,'--------------------------------------------'
133         PRINT*,' inifis: Parameters for the physics (callphys.def)'
134         PRINT*,'--------------------------------------------'
135
136         write(*,*) "Directory where external input files are:"
137         datafile="/data/fgglmd/datagcm/datafile"
138         call getin("datadir",datafile) ! default path
139         write(*,*) " datafile = ",trim(datafile)
140
141         write(*,*) "Run with or without tracer transport ?"
142         tracer=.false. ! default value
143         call getin("tracer",tracer)
144         write(*,*) " tracer = ",tracer
145
146         write(*,*) "Diurnal cycle ?"
147         write(*,*) "(if diurnal=False, diurnal averaged solar heating)"
148         diurnal=.true. ! default value
149         call getin("diurnal",diurnal)
150         write(*,*) " diurnal = ",diurnal
151
152         write(*,*) "Seasonal cycle ?"
153         write(*,*) "(if season=False, Ls stays constant, to value ",
154     &   "set in 'start'"
155         season=.true. ! default value
156         call getin("season",season)
157         write(*,*) " season = ",season
158
159         write(*,*) "Write some extra output to the screen ?"
160         lwrite=.false. ! default value
161         call getin("lwrite",lwrite)
162         write(*,*) " lwrite = ",lwrite
163
164         write(*,*) "Save statistics in file stats.nc ?"
165#ifdef MESOSCALE
166         callstats=.false. ! default value
167#else
168         callstats=.true. ! default value
169#endif
170         call getin("callstats",callstats)
171         write(*,*) " callstats = ",callstats
172
173         write(*,*) "Save EOF profiles in file 'profiles' for ",
174     &              "Climate Database?"
175         calleofdump=.false. ! default value
176         call getin("calleofdump",calleofdump)
177         write(*,*) " calleofdump = ",calleofdump
178
179         write(*,*) "Dust scenario: 1=constant dust (read from startfi",
180     &   " or set as tauvis); 2=Viking scenario; =3 MGS scenario,",
181     &   "=4 Mars Year 24 from TES assimilation, ",
182     &   "=24,25 or 26 :Mars Year 24,25 or 26 from TES assimilation"
183         iaervar=3 ! default value
184         call getin("iaervar",iaervar)
185         write(*,*) " iaervar = ",iaervar
186
187         write(*,*) "Reference (visible) dust opacity at 700 Pa ",
188     &   "(matters only if iaervar=1)"
189         ! NB: default value of tauvis is set/read in startfi.nc file
190         call getin("tauvis",tauvis)
191         write(*,*) " tauvis = ",tauvis
192
193         write(*,*) "Dust vertical distribution:"
194         write(*,*) "(=1 top set by topdustref parameter;",
195     & " =2 Viking scenario; =3 MGS scenario)"
196         iddist=3 ! default value
197         call getin("iddist",iddist)
198         write(*,*) " iddist = ",iddist
199
200         write(*,*) "Dust top altitude (km). (Matters only if iddist=1)"
201         topdustref= 90.0 ! default value
202         call getin("topdustref",topdustref)
203         write(*,*) " topdustref = ",topdustref
204
205         write(*,*) "call radiative transfer ?"
206         callrad=.true. ! default value
207         call getin("callrad",callrad)
208         write(*,*) " callrad = ",callrad
209
210         write(*,*) "call slope insolation scheme ?",
211     &              "(matters only if callrad=T)"
212#ifdef MESOSCALE
213         callslope=.true. ! default value
214#else
215         callslope=.false. ! default value (not supported yet)
216#endif
217         call getin("callslope",callslope)
218         write(*,*) " callslope = ",callslope
219
220         write(*,*) "call NLTE radiative schemes ?",
221     &              "(matters only if callrad=T)"
222         callnlte=.false. ! default value
223         call getin("callnlte",callnlte)
224         write(*,*) " callnlte = ",callnlte
225         
226         nltemodel=0    !default value
227         write(*,*) "NLTE model?"
228         write(*,*) "0 -> old model, static O"
229         write(*,*) "1 -> old model, dynamic O"
230         write(*,*) "2 -> new model"
231         write(*,*) "(matters only if callnlte=T)"
232         call getin("nltemodel",nltemodel)
233         write(*,*) " nltemodel = ",nltemodel
234
235         write(*,*) "call CO2 NIR absorption ?",
236     &              "(matters only if callrad=T)"
237         callnirco2=.false. ! default value
238         call getin("callnirco2",callnirco2)
239         write(*,*) " callnirco2 = ",callnirco2
240
241         write(*,*) "New NIR NLTE correction ?",
242     $              "0-> old model (no correction)",
243     $              "1-> new correction",
244     $              "(matters only if callnirco2=T)"
245         call getin("nircorr",nircorr)
246         write(*,*) " nircorr = ",nircorr
247
248         write(*,*) "call turbulent vertical diffusion ?"
249         calldifv=.true. ! default value
250         call getin("calldifv",calldifv)
251         write(*,*) " calldifv = ",calldifv
252
253         write(*,*) "call thermals ?"
254         calltherm=.false. ! default value
255         call getin("calltherm",calltherm)
256         write(*,*) " calltherm = ",calltherm
257
258         write(*,*) "output thermal diagnostics ?"
259         outptherm=.false. ! default value
260         call getin("outptherm",outptherm)
261         write(*,*) " outptherm = ",outptherm
262
263         write(*,*) "call convective adjustment ?"
264         calladj=.true. ! default value
265         call getin("calladj",calladj)
266         write(*,*) " calladj = ",calladj
267         
268         if (calltherm .and. (.not. calladj)) then
269          print*,'Convadj has to be activated when using thermals'
270          stop
271         endif
272
273         write(*,*) "call Richardson-based surface layer ?"
274         callrichsl=.false. ! default value
275         call getin("callrichsl",callrichsl)
276         write(*,*) " callrichsl = ",callrichsl
277
278         if (calltherm .and. .not.callrichsl) then
279          print*,'WARNING WARNING WARNING'
280          print*,'if calltherm=T we strongly advise that '
281          print*,'you use the new surface layer scheme '
282          print*,'by setting callrichsl=T '
283         endif
284
285         write(*,*) "call CO2 condensation ?"
286         callcond=.true. ! default value
287         call getin("callcond",callcond)
288         write(*,*) " callcond = ",callcond
289
290         write(*,*)"call thermal conduction in the soil ?"
291         callsoil=.true. ! default value
292         call getin("callsoil",callsoil)
293         write(*,*) " callsoil = ",callsoil
294         
295
296         write(*,*)"call Lott's gravity wave/subgrid topography ",
297     &             "scheme ?"
298         calllott=.true. ! default value
299         call getin("calllott",calllott)
300         write(*,*)" calllott = ",calllott
301
302
303         write(*,*)"rad.transfer is computed every iradia",
304     &             " physical timestep"
305         iradia=1 ! default value
306         call getin("iradia",iradia)
307         write(*,*)" iradia = ",iradia
308         
309
310         write(*,*)"Output of the exchange coefficient mattrix ?",
311     &             "(for diagnostics only)"
312         callg2d=.false. ! default value
313         call getin("callg2d",callg2d)
314         write(*,*)" callg2d = ",callg2d
315
316         write(*,*)"Rayleigh scattering : (should be .false. for now)"
317         rayleigh=.false.
318         call getin("rayleigh",rayleigh)
319         write(*,*)" rayleigh = ",rayleigh
320
321
322! TRACERS:
323
324! dustbin
325         write(*,*)"Transported dust ? (if >0, use 'dustbin' dust bins)"
326         dustbin=0 ! default value
327         call getin("dustbin",dustbin)
328         write(*,*)" dustbin = ",dustbin
329! active
330         write(*,*)"Radiatively active dust ? (matters if dustbin>0)"
331         active=.false. ! default value
332         call getin("active",active)
333         write(*,*)" active = ",active
334
335! Test of incompatibility:
336! if active is used, then dustbin should be > 0
337
338         if (active.and.(dustbin.lt.1)) then
339           print*,'if active is used, then dustbin should > 0'
340           stop
341         endif
342! doubleq
343         write(*,*)"use mass and number mixing ratios to predict",
344     &             " dust size ?"
345         doubleq=.false. ! default value
346         call getin("doubleq",doubleq)
347         write(*,*)" doubleq = ",doubleq
348! submicron
349         submicron=.false. ! default value
350         call getin("submicron",submicron)
351         write(*,*)" submicron = ",submicron
352
353! Test of incompatibility:
354! if doubleq is used, then dustbin should be 2
355
356         if (doubleq.and.(dustbin.ne.2)) then
357           print*,'if doubleq is used, then dustbin should be 2'
358           stop
359         endif
360         if (doubleq.and.submicron.and.(nqmx.LT.3)) then
361           print*,'If doubleq is used with a submicron tracer,'
362           print*,' then the number of tracers has to be'
363           print*,' larger than 3.'
364           stop
365         endif
366! lifting
367         write(*,*)"dust lifted by GCM surface winds ?"
368         lifting=.false. ! default value
369         call getin("lifting",lifting)
370         write(*,*)" lifting = ",lifting
371
372! Test of incompatibility:
373! if lifting is used, then dustbin should be > 0
374
375         if (lifting.and.(dustbin.lt.1)) then
376           print*,'if lifting is used, then dustbin should > 0'
377           stop
378         endif
379! callddevil
380         write(*,*)" dust lifted by dust devils ?"
381         callddevil=.false. !default value
382         call getin("callddevil",callddevil)
383         write(*,*)" callddevil = ",callddevil
384
385! Test of incompatibility:
386! if dustdevil is used, then dustbin should be > 0
387
388         if (callddevil.and.(dustbin.lt.1)) then
389           print*,'if dustdevil is used, then dustbin should > 0'
390           stop
391         endif
392! sedimentation
393         write(*,*) "Gravitationnal sedimentation ?"
394         sedimentation=.true. ! default value
395         call getin("sedimentation",sedimentation)
396         write(*,*) " sedimentation = ",sedimentation
397! activice
398         write(*,*) "Radiatively active transported atmospheric ",
399     &              "water ice ?"
400         activice=.false. ! default value
401         call getin("activice",activice)
402         write(*,*) " activice = ",activice
403! water
404         write(*,*) "Compute water cycle ?"
405         water=.false. ! default value
406         call getin("water",water)
407         write(*,*) " water = ",water
408
409! Test of incompatibility:
410
411         if (activice.and..not.water) then
412           print*,'if activice is used, water should be used too'
413           stop
414         endif
415
416         if (water.and..not.tracer) then
417           print*,'if water is used, tracer should be used too'
418           stop
419         endif
420         
421! ccn factor if no scavenging         
422        write(*,*) "water param CCN reduc. factor ?", ccn_factor
423        ccn_factor = 4.5
424        call getin("ccn_factor",ccn_factor)
425        write(*,*)" ccn_factor = ",ccn_factor
426        write(*,*)"Careful: only used when microphys=F, otherwise"
427        write(*,*)"the contact parameter is used instead;"
428
429! microphys
430         write(*,*)"Microphysical scheme for water-ice clouds?"
431         microphys=.false. ! default value
432         call getin("microphys",microphys)
433         write(*,*)" microphys = ",microphys
434
435! microphysical parameter contact       
436         write(*,*) "water contact parameter ?"
437         mteta  = 0.95
438         call getin("mteta",mteta)
439         write(*,*) "mteta = ", mteta
440
441! scavenging
442         write(*,*)"Dust scavenging by H2O/CO2 snowfall ?"
443         scavenging=.false. ! default value
444         call getin("scavenging",scavenging)
445         write(*,*)" scavenging = ",scavenging
446         
447
448! Test of incompatibility:
449! if scavenging is used, then dustbin should be > 0
450
451         if (scavenging.and.(dustbin.lt.1)) then
452           print*,'if scavenging is used, then dustbin should > 0'
453           stop
454         endif
455         if ((microphys.and..not.doubleq).or.
456     &       (microphys.and..not.active).or.
457     &       (microphys.and..not.scavenging).or.
458     &       (microphys.and..not.water)) then
459           print*,'if microphys is used, then doubleq,'
460           print*,'active, water and scavenging must all be used!'
461           stop
462         endif
463         if ((scavenging.and..not.doubleq).or.
464     &       (scavenging.and..not.active).or.
465     &       (scavenging.and..not.water).or.
466     &       (scavenging.and..not.microphys)) then
467           print*,'if scavenging is used, then doubleq,'
468           print*,'active, water and microphys must all be used!'
469           stop
470         endif
471
472! Test of incompatibility:
473
474         write(*,*) "Permanent water caps at poles ?",
475     &               " .true. is RECOMMENDED"
476         write(*,*) "(with .true., North cap is a source of water ",
477     &   "and South pole is a cold trap)"
478         caps=.true. ! default value
479         call getin("caps",caps)
480         write(*,*) " caps = ",caps
481
482! albedo_h2o_ice
483         write(*,*) "water ice albedo ?"
484         albedo_h2o_ice=0.45
485         call getin("albedo_h2o_ice",albedo_h2o_ice)
486         write(*,*) " albedo_h2o_ice = ",albedo_h2o_ice
487! inert_h2o_ice
488         write(*,*) "water ice thermal inertia ?"
489         inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
490         call getin("inert_h2o_ice",inert_h2o_ice)
491         write(*,*) " inert_h2o_ice = ",inert_h2o_ice
492! frost_albedo_threshold
493         write(*,*) "frost thickness threshold for albedo ?"
494         frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2)
495         call getin("frost_albedo_threshold",
496     &    frost_albedo_threshold)
497         write(*,*) " frost_albedo_threshold = ",
498     &            frost_albedo_threshold
499
500!!!!!!!!!!!!!!!! TEMP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
501!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     
502         write(*,*) "temp tag for water caps ?"
503         temptag=.false.
504         call getin("temptag",temptag)
505         write(*,*) " temptag = ",temptag
506!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
507!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!             
508
509         write(*,*) "photochemistry: include chemical species"
510         photochem=.false. ! default value
511         call getin("photochem",photochem)
512         write(*,*) " photochem = ",photochem
513
514
515! THERMOSPHERE
516
517         write(*,*) "call thermosphere ?"
518         callthermos=.false. ! default value
519         call getin("callthermos",callthermos)
520         write(*,*) " callthermos = ",callthermos
521         
522
523         write(*,*) " water included without cycle ",
524     &              "(only if water=.false.)"
525         thermoswater=.false. ! default value
526         call getin("thermoswater",thermoswater)
527         write(*,*) " thermoswater = ",thermoswater
528
529         write(*,*) "call thermal conduction ?",
530     &    " (only if callthermos=.true.)"
531         callconduct=.false. ! default value
532         call getin("callconduct",callconduct)
533         write(*,*) " callconduct = ",callconduct
534
535         write(*,*) "call EUV heating ?",
536     &   " (only if callthermos=.true.)"
537         calleuv=.false.  ! default value
538         call getin("calleuv",calleuv)
539         write(*,*) " calleuv = ",calleuv
540
541         write(*,*) "call molecular viscosity ?",
542     &   " (only if callthermos=.true.)"
543         callmolvis=.false. ! default value
544         call getin("callmolvis",callmolvis)
545         write(*,*) " callmolvis = ",callmolvis
546
547         write(*,*) "call molecular diffusion ?",
548     &   " (only if callthermos=.true.)"
549         callmoldiff=.false. ! default value
550         call getin("callmoldiff",callmoldiff)
551         write(*,*) " callmoldiff = ",callmoldiff
552         
553
554         write(*,*) "call thermospheric photochemistry ?",
555     &   " (only if callthermos=.true.)"
556         thermochem=.false. ! default value
557         call getin("thermochem",thermochem)
558         write(*,*) " thermochem = ",thermochem
559
560         write(*,*) "date for solar flux calculation:",
561     &   " (1985 < date < 2002)"
562         write(*,*) "(Solar min=1996.4 ave=1993.4 max=1990.6)"
563         solarcondate=1993.4 ! default value
564         call getin("solarcondate",solarcondate)
565         write(*,*) " solarcondate = ",solarcondate
566         
567
568         if (.not.callthermos) then
569           if (thermoswater) then
570             print*,'if thermoswater is set, callthermos must be true'
571             stop
572           endif         
573           if (callconduct) then
574             print*,'if callconduct is set, callthermos must be true'
575             stop
576           endif       
577           if (calleuv) then
578             print*,'if calleuv is set, callthermos must be true'
579             stop
580           endif         
581           if (callmolvis) then
582             print*,'if callmolvis is set, callthermos must be true'
583             stop
584           endif       
585           if (callmoldiff) then
586             print*,'if callmoldiff is set, callthermos must be true'
587             stop
588           endif         
589           if (thermochem) then
590             print*,'if thermochem is set, callthermos must be true'
591             stop
592           endif         
593        endif
594
595! Test of incompatibility:
596! if photochem is used, then water should be used too
597
598         if (photochem.and..not.water) then
599           print*,'if photochem is used, water should be used too'
600           stop
601         endif
602
603! if callthermos is used, then thermoswater should be used too
604! (if water not used already)
605
606         if (callthermos .and. .not.water) then
607           if (callthermos .and. .not.thermoswater) then
608             print*,'if callthermos is used, water or thermoswater
609     &               should be used too'
610             stop
611           endif
612         endif
613
614         PRINT*,'--------------------------------------------'
615         PRINT*
616         PRINT*
617      ELSE
618         write(*,*)
619         write(*,*) 'Cannot read file callphys.def. Is it here ?'
620         stop
621      ENDIF
622
6238000  FORMAT(t5,a12,l8)
6248001  FORMAT(t5,a12,i8)
625
626      PRINT*
627      PRINT*,'inifis: daysec',daysec
628      PRINT*
629      PRINT*,'inifis: The radiative transfer is computed:'
630      PRINT*,'           each ',iradia,' physical time-step'
631      PRINT*,'        or each ',iradia*dtphys,' seconds'
632      PRINT*
633! --------------------------------------------------------------
634!  Managing the Longwave radiative transfer
635! --------------------------------------------------------------
636
637!     In most cases, the run just use the following values :
638!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
639      callemis=.true.     
640!     ilwd=10*int(daysec/dtphys) ! bug before 22/10/01       
641      ilwd=1
642      ilwn=1 !2
643      ilwb=1 !2
644      linear=.true.       
645      ncouche=3
646      alphan=0.4
647      semi=0
648
649!     BUT people working hard on the LW may want to read them in 'radia.def'
650!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
651
652      OPEN(99,file='radia.def',status='old',form='formatted'
653     .     ,iostat=ierr)
654      IF(ierr.EQ.0) THEN
655         write(*,*) 'inifis: Reading radia.def !!!'
656         READ(99,fmt='(a)') ch1
657         READ(99,*) callemis
658         WRITE(*,8000) ch1,callemis
659
660         READ(99,fmt='(a)') ch1
661         READ(99,*) iradia
662         WRITE(*,8001) ch1,iradia
663
664         READ(99,fmt='(a)') ch1
665         READ(99,*) ilwd
666         WRITE(*,8001) ch1,ilwd
667
668         READ(99,fmt='(a)') ch1
669         READ(99,*) ilwn
670         WRITE(*,8001) ch1,ilwn
671
672         READ(99,fmt='(a)') ch1
673         READ(99,*) linear
674         WRITE(*,8000) ch1,linear
675
676         READ(99,fmt='(a)') ch1
677         READ(99,*) ncouche
678         WRITE(*,8001) ch1,ncouche
679
680         READ(99,fmt='(a)') ch1
681         READ(99,*) alphan
682         WRITE(*,*) ch1,alphan
683
684         READ(99,fmt='(a)') ch1
685         READ(99,*) ilwb
686         WRITE(*,8001) ch1,ilwb
687
688
689         READ(99,fmt='(a)') ch1
690         READ(99,'(l1)') callg2d
691         WRITE(*,8000) ch1,callg2d
692
693         READ(99,fmt='(a)') ch1
694         READ(99,*) semi
695         WRITE(*,*) ch1,semi
696      end if
697      CLOSE(99)
698
699!-----------------------------------------------------------------------
700!     Some more initialization:
701!     ------------------------
702
703      ! in 'comgeomfi.h'
704      CALL SCOPY(ngrid,plon,1,long,1)
705      CALL SCOPY(ngrid,plat,1,lati,1)
706      CALL SCOPY(ngrid,parea,1,area,1)
707      totarea=SSUM(ngridmx,area,1)
708
709      ! in 'comdiurn.h'
710      DO ig=1,ngrid
711         sinlat(ig)=sin(plat(ig))
712         coslat(ig)=cos(plat(ig))
713         sinlon(ig)=sin(plon(ig))
714         coslon(ig)=cos(plon(ig))
715      ENDDO
716
717      pi=2.*asin(1.) ! NB: pi is a common in comcstfi.h
718
719!     managing the tracers, and tests:
720!     -------------------------------
721!     Ehouarn: removed; as these tests are now done in initracer.F
722!      if(tracer) then
723!
724!!          when photochem is used, nqchem_min is the rank
725!!          of the first chemical species
726!
727!! Ehouarn: nqchem_min is now meaningless and no longer used
728!!       nqchem_min = 1
729!       if (photochem .or. callthermos) then
730!         chem = .true.
731!       end if
732!
733!       if (water .or. thermoswater) h2o = .true.
734!
735!!          TESTS
736!
737!       print*,'inifis: TRACERS:'
738!       write(*,*) "    chem=",chem,"    h2o=",h2o
739!!       write(*,*) "   doubleq=",doubleq
740!!       write(*,*) "   dustbin=",dustbin
741!
742!       if ((doubleq).and.(h2o).and.
743!     $     (chem)) then
744!         print*,' 2 dust tracers (doubleq)'
745!         print*,' 1 water vapour tracer'
746!         print*,' 1 water ice tracer'
747!         print*,nqmx-4,' chemistry tracers'
748!       endif
749!
750!       if ((doubleq).and.(h2o).and.
751!     $     .not.(chem)) then
752!         print*,' 2 dust tracers (doubleq)'
753!         print*,' 1 water vapour tracer'
754!         print*,' 1 water ice tracer'
755!         if (nqmx.LT.4) then
756!           print*,'nqmx should be at least equal to'
757!           print*,'4 with these options.'
758!           stop
759!         endif
760!       endif
761!
762!       if (.not.(doubleq).and.(h2o).and.
763!     $     (chem)) then
764!         if (dustbin.gt.0) then
765!           print*,dustbin,' dust bins'
766!         endif
767!         print*,nqmx-2-dustbin,' chemistry tracers'
768!         print*,' 1 water vapour tracer'
769!         print*,' 1 water ice tracer'
770!       endif
771!
772!       if (.not.(doubleq).and.(h2o).and.
773!     $     .not.(chem)) then
774!         if (dustbin.gt.0) then
775!           print*,dustbin,' dust bins'
776!         endif
777!         print*,' 1 water vapour tracer'
778!         print*,' 1 water ice tracer'
779!         if (nqmx.gt.(dustbin+2)) then
780!           print*,'nqmx should be ',(dustbin+2),
781!     $            ' with these options...'
782!                  print*,'(or check callphys.def)'
783!         endif
784!         if (nqmx.lt.(dustbin+2)) then
785!           write(*,*) "inifis: nqmx.lt.(dustbin+2)"
786!           stop
787!         endif
788!       endif
789!
790!      endif ! of if (tracer)
791!
792!      RETURN
793      END
Note: See TracBrowser for help on using the repository browser.