source: trunk/mesoscale/POUR_LIBF_COMMUN/meso_inifis.F @ 35

Last change on this file since 35 was 35, checked in by aslmd, 15 years ago

LMD_MM_MARS: travail de preparation a un dossier libf commun avec LMDZ.MARS

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