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

Last change on this file since 481 was 455, checked in by tnavarro, 14 years ago

tiny change : nuice_sed initialisation is now done in inifis. Also changed initracer and improvedclouds.

File size: 25.6 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="/u/forget/WWW/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! water ice clouds effective variance distribution for sedimentaion       
422        write(*,*) "effective variance for water ice clouds ?"
423        nuice_sed=0.45
424        call getin("nuice_sed",nuice_sed)
425        write(*,*) "water_param nueff Sedimentation:", nuice_sed
426         
427! ccn factor if no scavenging         
428        write(*,*) "water param CCN reduc. factor ?", ccn_factor
429        ccn_factor = 4.5
430        call getin("ccn_factor",ccn_factor)
431        write(*,*)" ccn_factor = ",ccn_factor
432        write(*,*)"Careful: only used when microphys=F, otherwise"
433        write(*,*)"the contact parameter is used instead;"
434
435! microphys
436         write(*,*)"Microphysical scheme for water-ice clouds?"
437         microphys=.false. ! default value
438         call getin("microphys",microphys)
439         write(*,*)" microphys = ",microphys
440
441! microphysical parameter contact       
442         write(*,*) "water contact parameter ?"
443         mteta  = 0.95
444         call getin("mteta",mteta)
445         write(*,*) "mteta = ", mteta
446
447! scavenging
448         write(*,*)"Dust scavenging by H2O/CO2 snowfall ?"
449         scavenging=.false. ! default value
450         call getin("scavenging",scavenging)
451         write(*,*)" scavenging = ",scavenging
452         
453
454! Test of incompatibility:
455! if scavenging is used, then dustbin should be > 0
456
457         if (scavenging.and.(dustbin.lt.1)) then
458           print*,'if scavenging is used, then dustbin should > 0'
459           stop
460         endif
461         if ((microphys.and..not.doubleq).or.
462     &       (microphys.and..not.active).or.
463     &       (microphys.and..not.scavenging).or.
464     &       (microphys.and..not.water)) then
465           print*,'if microphys is used, then doubleq,'
466           print*,'active, water and scavenging must all be used!'
467           stop
468         endif
469         if ((scavenging.and..not.doubleq).or.
470     &       (scavenging.and..not.active).or.
471     &       (scavenging.and..not.water).or.
472     &       (scavenging.and..not.microphys)) then
473           print*,'if scavenging is used, then doubleq,'
474           print*,'active, water and microphys must all be used!'
475           stop
476         endif
477
478! Test of incompatibility:
479
480         write(*,*) "Permanent water caps at poles ?",
481     &               " .true. is RECOMMENDED"
482         write(*,*) "(with .true., North cap is a source of water ",
483     &   "and South pole is a cold trap)"
484         caps=.true. ! default value
485         call getin("caps",caps)
486         write(*,*) " caps = ",caps
487
488! albedo_h2o_ice
489         write(*,*) "water ice albedo ?"
490         albedo_h2o_ice=0.45
491         call getin("albedo_h2o_ice",albedo_h2o_ice)
492         write(*,*) " albedo_h2o_ice = ",albedo_h2o_ice
493! inert_h2o_ice
494         write(*,*) "water ice thermal inertia ?"
495         inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
496         call getin("inert_h2o_ice",inert_h2o_ice)
497         write(*,*) " inert_h2o_ice = ",inert_h2o_ice
498! frost_albedo_threshold
499         write(*,*) "frost thickness threshold for albedo ?"
500         frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2)
501         call getin("frost_albedo_threshold",
502     &    frost_albedo_threshold)
503         write(*,*) " frost_albedo_threshold = ",
504     &            frost_albedo_threshold
505
506!!!!!!!!!!!!!!!! TEMP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
507!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     
508         write(*,*) "temp tag for water caps ?"
509         temptag=.false.
510         call getin("temptag",temptag)
511         write(*,*) " temptag = ",temptag
512!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
513!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!             
514
515         write(*,*) "photochemistry: include chemical species"
516         photochem=.false. ! default value
517         call getin("photochem",photochem)
518         write(*,*) " photochem = ",photochem
519
520
521! THERMOSPHERE
522
523         write(*,*) "call thermosphere ?"
524         callthermos=.false. ! default value
525         call getin("callthermos",callthermos)
526         write(*,*) " callthermos = ",callthermos
527         
528
529         write(*,*) " water included without cycle ",
530     &              "(only if water=.false.)"
531         thermoswater=.false. ! default value
532         call getin("thermoswater",thermoswater)
533         write(*,*) " thermoswater = ",thermoswater
534
535         write(*,*) "call thermal conduction ?",
536     &    " (only if callthermos=.true.)"
537         callconduct=.false. ! default value
538         call getin("callconduct",callconduct)
539         write(*,*) " callconduct = ",callconduct
540
541         write(*,*) "call EUV heating ?",
542     &   " (only if callthermos=.true.)"
543         calleuv=.false.  ! default value
544         call getin("calleuv",calleuv)
545         write(*,*) " calleuv = ",calleuv
546
547         write(*,*) "call molecular viscosity ?",
548     &   " (only if callthermos=.true.)"
549         callmolvis=.false. ! default value
550         call getin("callmolvis",callmolvis)
551         write(*,*) " callmolvis = ",callmolvis
552
553         write(*,*) "call molecular diffusion ?",
554     &   " (only if callthermos=.true.)"
555         callmoldiff=.false. ! default value
556         call getin("callmoldiff",callmoldiff)
557         write(*,*) " callmoldiff = ",callmoldiff
558         
559
560         write(*,*) "call thermospheric photochemistry ?",
561     &   " (only if callthermos=.true.)"
562         thermochem=.false. ! default value
563         call getin("thermochem",thermochem)
564         write(*,*) " thermochem = ",thermochem
565
566         write(*,*) "date for solar flux calculation:",
567     &   " (1985 < date < 2002)"
568         write(*,*) "(Solar min=1996.4 ave=1993.4 max=1990.6)"
569         solarcondate=1993.4 ! default value
570         call getin("solarcondate",solarcondate)
571         write(*,*) " solarcondate = ",solarcondate
572         
573
574         if (.not.callthermos) then
575           if (thermoswater) then
576             print*,'if thermoswater is set, callthermos must be true'
577             stop
578           endif         
579           if (callconduct) then
580             print*,'if callconduct is set, callthermos must be true'
581             stop
582           endif       
583           if (calleuv) then
584             print*,'if calleuv is set, callthermos must be true'
585             stop
586           endif         
587           if (callmolvis) then
588             print*,'if callmolvis is set, callthermos must be true'
589             stop
590           endif       
591           if (callmoldiff) then
592             print*,'if callmoldiff is set, callthermos must be true'
593             stop
594           endif         
595           if (thermochem) then
596             print*,'if thermochem is set, callthermos must be true'
597             stop
598           endif         
599        endif
600
601! Test of incompatibility:
602! if photochem is used, then water should be used too
603
604         if (photochem.and..not.water) then
605           print*,'if photochem is used, water should be used too'
606           stop
607         endif
608
609! if callthermos is used, then thermoswater should be used too
610! (if water not used already)
611
612         if (callthermos .and. .not.water) then
613           if (callthermos .and. .not.thermoswater) then
614             print*,'if callthermos is used, water or thermoswater
615     &               should be used too'
616             stop
617           endif
618         endif
619
620         PRINT*,'--------------------------------------------'
621         PRINT*
622         PRINT*
623      ELSE
624         write(*,*)
625         write(*,*) 'Cannot read file callphys.def. Is it here ?'
626         stop
627      ENDIF
628
6298000  FORMAT(t5,a12,l8)
6308001  FORMAT(t5,a12,i8)
631
632      PRINT*
633      PRINT*,'inifis: daysec',daysec
634      PRINT*
635      PRINT*,'inifis: The radiative transfer is computed:'
636      PRINT*,'           each ',iradia,' physical time-step'
637      PRINT*,'        or each ',iradia*dtphys,' seconds'
638      PRINT*
639! --------------------------------------------------------------
640!  Managing the Longwave radiative transfer
641! --------------------------------------------------------------
642
643!     In most cases, the run just use the following values :
644!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
645      callemis=.true.     
646!     ilwd=10*int(daysec/dtphys) ! bug before 22/10/01       
647      ilwd=1
648      ilwn=1 !2
649      ilwb=1 !2
650      linear=.true.       
651      ncouche=3
652      alphan=0.4
653      semi=0
654
655!     BUT people working hard on the LW may want to read them in 'radia.def'
656!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
657
658      OPEN(99,file='radia.def',status='old',form='formatted'
659     .     ,iostat=ierr)
660      IF(ierr.EQ.0) THEN
661         write(*,*) 'inifis: Reading radia.def !!!'
662         READ(99,fmt='(a)') ch1
663         READ(99,*) callemis
664         WRITE(*,8000) ch1,callemis
665
666         READ(99,fmt='(a)') ch1
667         READ(99,*) iradia
668         WRITE(*,8001) ch1,iradia
669
670         READ(99,fmt='(a)') ch1
671         READ(99,*) ilwd
672         WRITE(*,8001) ch1,ilwd
673
674         READ(99,fmt='(a)') ch1
675         READ(99,*) ilwn
676         WRITE(*,8001) ch1,ilwn
677
678         READ(99,fmt='(a)') ch1
679         READ(99,*) linear
680         WRITE(*,8000) ch1,linear
681
682         READ(99,fmt='(a)') ch1
683         READ(99,*) ncouche
684         WRITE(*,8001) ch1,ncouche
685
686         READ(99,fmt='(a)') ch1
687         READ(99,*) alphan
688         WRITE(*,*) ch1,alphan
689
690         READ(99,fmt='(a)') ch1
691         READ(99,*) ilwb
692         WRITE(*,8001) ch1,ilwb
693
694
695         READ(99,fmt='(a)') ch1
696         READ(99,'(l1)') callg2d
697         WRITE(*,8000) ch1,callg2d
698
699         READ(99,fmt='(a)') ch1
700         READ(99,*) semi
701         WRITE(*,*) ch1,semi
702      end if
703      CLOSE(99)
704
705!-----------------------------------------------------------------------
706!     Some more initialization:
707!     ------------------------
708
709      ! in 'comgeomfi.h'
710      CALL SCOPY(ngrid,plon,1,long,1)
711      CALL SCOPY(ngrid,plat,1,lati,1)
712      CALL SCOPY(ngrid,parea,1,area,1)
713      totarea=SSUM(ngridmx,area,1)
714
715      ! in 'comdiurn.h'
716      DO ig=1,ngrid
717         sinlat(ig)=sin(plat(ig))
718         coslat(ig)=cos(plat(ig))
719         sinlon(ig)=sin(plon(ig))
720         coslon(ig)=cos(plon(ig))
721      ENDDO
722
723      pi=2.*asin(1.) ! NB: pi is a common in comcstfi.h
724
725!     managing the tracers, and tests:
726!     -------------------------------
727!     Ehouarn: removed; as these tests are now done in initracer.F
728!      if(tracer) then
729!
730!!          when photochem is used, nqchem_min is the rank
731!!          of the first chemical species
732!
733!! Ehouarn: nqchem_min is now meaningless and no longer used
734!!       nqchem_min = 1
735!       if (photochem .or. callthermos) then
736!         chem = .true.
737!       end if
738!
739!       if (water .or. thermoswater) h2o = .true.
740!
741!!          TESTS
742!
743!       print*,'inifis: TRACERS:'
744!       write(*,*) "    chem=",chem,"    h2o=",h2o
745!!       write(*,*) "   doubleq=",doubleq
746!!       write(*,*) "   dustbin=",dustbin
747!
748!       if ((doubleq).and.(h2o).and.
749!     $     (chem)) then
750!         print*,' 2 dust tracers (doubleq)'
751!         print*,' 1 water vapour tracer'
752!         print*,' 1 water ice tracer'
753!         print*,nqmx-4,' chemistry tracers'
754!       endif
755!
756!       if ((doubleq).and.(h2o).and.
757!     $     .not.(chem)) then
758!         print*,' 2 dust tracers (doubleq)'
759!         print*,' 1 water vapour tracer'
760!         print*,' 1 water ice tracer'
761!         if (nqmx.LT.4) then
762!           print*,'nqmx should be at least equal to'
763!           print*,'4 with these options.'
764!           stop
765!         endif
766!       endif
767!
768!       if (.not.(doubleq).and.(h2o).and.
769!     $     (chem)) then
770!         if (dustbin.gt.0) then
771!           print*,dustbin,' dust bins'
772!         endif
773!         print*,nqmx-2-dustbin,' chemistry tracers'
774!         print*,' 1 water vapour tracer'
775!         print*,' 1 water ice tracer'
776!       endif
777!
778!       if (.not.(doubleq).and.(h2o).and.
779!     $     .not.(chem)) then
780!         if (dustbin.gt.0) then
781!           print*,dustbin,' dust bins'
782!         endif
783!         print*,' 1 water vapour tracer'
784!         print*,' 1 water ice tracer'
785!         if (nqmx.gt.(dustbin+2)) then
786!           print*,'nqmx should be ',(dustbin+2),
787!     $            ' with these options...'
788!                  print*,'(or check callphys.def)'
789!         endif
790!         if (nqmx.lt.(dustbin+2)) then
791!           write(*,*) "inifis: nqmx.lt.(dustbin+2)"
792!           stop
793!         endif
794!       endif
795!
796!      endif ! of if (tracer)
797!
798!      RETURN
799      END
Note: See TracBrowser for help on using the repository browser.