source: trunk/mesoscale/LMDZ.MARS/libf_gcm/phymars/inifis.F @ 113

Last change on this file since 113 was 57, checked in by aslmd, 14 years ago

mineur LMD_MM_MARS: ajout du GCM ancienne physique, systeme maintenant complet sur SVN (ne manque que la base de donnees d'etats initiaux)
File size: 16.4 KB
RevLine 
[57]1      SUBROUTINE inifis(ngrid,nlayer,
2     $           day_ini,pdaysec,ptimestep,
3     $           plat,plon,parea,
4     $           prad,pg,pr,pcpp)
5      IMPLICIT NONE
6c
7c=======================================================================
8c
9c   subject:
10c   --------
11c
12c   Initialisation for the physical parametrisations of the LMD
13c   martian atmospheric general circulation modele.
14c
15c   author: Frederic Hourdin 15 / 10 /93
16c   -------
17c   modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
18c
19c
20c   arguments:
21c   ----------
22c
23c   input:
24c   ------
25c
26c    ngrid                 Size of the horizontal grid.
27c                          All internal loops are performed on that grid.
28c    nlayer                Number of vertical layers.
29c    pdayref               Day of reference for the simulation
30c    firstcall             True at the first call
31c    lastcall              True at the last call
32c    pday                  Number of days counted from the North. Spring
33c                          equinoxe.
34c
35c=======================================================================
36c
37c-----------------------------------------------------------------------
38c   declarations:
39c   -------------
40 
41#include "dimensions.h"
42#include "dimphys.h"
43#include "planete.h"
44#include "comcstfi.h"
45#include "comsaison.h"
46#include "comdiurn.h"
47#include "comgeomfi.h"
48#include "callkeys.h"
49#include "surfdat.h"
50
51
52      REAL prad,pg,pr,pcpp,pdaysec,ptimestep
53 
54      INTEGER ngrid,nlayer
55      REAL plat(ngrid),plon(ngrid),parea(ngridmx)
56      integer day_ini
57      INTEGER ig,ierr
58 
59      EXTERNAL iniorbit,orbite
60      EXTERNAL SSUM
61      REAL SSUM
62 
63      CHARACTER ch1*12
64      CHARACTER ch80*80
65
66      logical chem, h2o
67
68      rad=prad
69      daysec=pdaysec
70      dtphys=ptimestep
71      cpp=pcpp
72      g=pg
73      r=pr
74      rcp=r/cpp
75
76c --------------------------------------------------------
77c     The usual Tests
78c     --------------
79      IF (nlayer.NE.nlayermx) THEN
80         PRINT*,'STOP in inifis'
81         PRINT*,'Probleme de dimensions :'
82         PRINT*,'nlayer     = ',nlayer
83         PRINT*,'nlayermx   = ',nlayermx
84         STOP
85      ENDIF
86
87      IF (ngrid.NE.ngridmx) THEN
88         PRINT*,'STOP in inifis'
89         PRINT*,'Probleme de dimensions :'
90         PRINT*,'ngrid     = ',ngrid
91         PRINT*,'ngridmx   = ',ngridmx
92         STOP
93      ENDIF
94
95c --------------------------------------------------------------
96c  Reading the "callphys.def" file controlling some key options
97c --------------------------------------------------------------
98
99      callrad=.true.
100      calldifv=.true.
101      calladj=.true.
102      callcond=.true.
103      callsoil=.true.
104      season=.true.
105      diurnal=.false.
106      lwrite=.false.
107      calllott=.true.
108      iaervar=2
109      iddist=3
110      topdustref=55.
111      OPEN(99,file='callphys.def',status='old',form='formatted'
112     .     ,iostat=ierr)
113      IF(ierr.EQ.0) THEN
114         PRINT*
115         PRINT*
116         PRINT*,'--------------------------------------------'
117         PRINT*,' Parametres pour la physique (callphys.def)'
118         PRINT*,'--------------------------------------------'
119
120         READ(99,*)
121         READ(99,*)
122
123         READ(99,fmt='(a)') ch1
124         READ(99,*) tracer
125         WRITE(*,8000) ch1,tracer
126
127         READ(99,fmt='(a)') ch1
128         READ(99,'(l1)') diurnal
129         WRITE(*,8000) ch1,diurnal
130
131         READ(99,fmt='(a)') ch1
132         READ(99,'(l1)') season
133         WRITE(*,8000) ch1,season
134
135         READ(99,fmt='(a)') ch1
136         READ(99,'(l1)') lwrite
137         WRITE(*,8000) ch1,lwrite
138
139         READ(99,fmt='(a)') ch1
140         READ(99,'(l1)') callstats
141         WRITE(*,8000) ch1,callstats
142
143         READ(99,fmt='(a)') ch1
144         READ(99,'(l1)') calleofdump
145         WRITE(*,8000) ch1,calleofdump
146
147         READ(99,*)
148         READ(99,*)
149
150         READ(99,fmt='(a)') ch1
151         READ(99,*,iostat=ierr) iaervar
152         if(ierr.ne.0) stop'Can t read iaervar in callphys.def (old?)'
153         WRITE(*,8001) ch1,iaervar
154
155         READ(99,fmt='(a)') ch1
156         READ(99,*) iddist
157         WRITE(*,8001) ch1,iddist
158
159         READ(99,fmt='(a)') ch1
160         READ(99,*) topdustref
161         WRITE(*,8002) ch1,topdustref
162
163         READ(99,*)
164         READ(99,*)
165
166         READ(99,fmt='(a)') ch1
167         READ(99,'(l1)') callrad
168         WRITE(*,8000) ch1,callrad
169
170         READ(99,fmt='(a)') ch1
171         READ(99,'(l1)') callnlte
172         WRITE(*,8000) ch1,callnlte
173         
174         READ(99,fmt='(a)') ch1
175         READ(99,'(l1)') callnirco2
176         WRITE(*,8000) ch1,callnirco2
177
178         READ(99,fmt='(a)') ch1
179         READ(99,'(l1)') calldifv
180         WRITE(*,8000) ch1,calldifv
181
182         READ(99,fmt='(a)') ch1
183         READ(99,'(l1)') calladj
184         WRITE(*,8000) ch1,calladj
185
186         READ(99,fmt='(a)') ch1
187         READ(99,'(l1)') callcond
188         WRITE(*,8000) ch1,callcond
189
190         READ(99,fmt='(a)') ch1
191         READ(99,'(l1)') callsoil
192         WRITE(*,8000) ch1,callsoil
193
194         READ(99,fmt='(a)') ch1
195         READ(99,'(l1)') calllott
196         WRITE(*,8000) ch1,calllott
197
198         READ(99,*)
199         READ(99,*)
200
201         READ(99,fmt='(a)') ch1
202         READ(99,*) iradia
203         WRITE(*,8001) ch1,iradia
204
205         READ(99,fmt='(a)') ch1
206         READ(99,'(l1)') callg2d
207         WRITE(*,8000) ch1,callg2d
208
209         READ(99,fmt='(a)') ch1
210         READ(99,*) rayleigh
211         WRITE(*,8000) ch1,rayleigh
212
213         READ(99,*)
214         READ(99,*)
215
216c TRACERS:
217
218         READ(99,fmt='(a)') ch1
219         READ(99,*) dustbin
220         WRITE(*,8001) ch1,dustbin
221
222         READ(99,fmt='(a)') ch1
223         READ(99,*) active
224         WRITE(*,8000) ch1,active
225
226c Test of incompatibility:
227c if active is used, then dustbin should be > 0
228
229         if (active.and.(dustbin.lt.1)) then
230           print*,'if active is used, then dustbin should > 0'
231           stop
232         endif
233
234         READ(99,fmt='(a)') ch1
235         READ(99,*) doubleq
236         WRITE(*,8000) ch1,doubleq
237
238c Test of incompatibility:
239c if doubleq is used, then dustbin should be 1
240
241         if (doubleq.and.(dustbin.ne.1)) then
242           print*,'if doubleq is used, then dustbin should be 1'
243           stop
244         endif
245
246         READ(99,fmt='(a)') ch1
247         READ(99,*) lifting
248         WRITE(*,8000) ch1,lifting
249
250c Test of incompatibility:
251c if lifting is used, then dustbin should be > 0
252
253         if (lifting.and.(dustbin.lt.1)) then
254           print*,'if lifting is used, then dustbin should > 0'
255           stop
256         endif
257
258         READ(99,fmt='(a)') ch1
259         READ(99,*) callddevil
260         WRITE(*,8000) ch1,callddevil
261
262c Test of incompatibility:
263c if dustdevil is used, then dustbin should be > 0
264
265         if (callddevil.and.(dustbin.lt.1)) then
266           print*,'if dustdevil is used, then dustbin should > 0'
267           stop
268         endif
269
270         READ(99,fmt='(a)') ch1
271         READ(99,*) scavenging
272         WRITE(*,8000) ch1,scavenging
273
274c Test of incompatibility:
275c if scavenging is used, then dustbin should be > 0
276
277         if (scavenging.and.(dustbin.lt.1)) then
278           print*,'if scavenging is used, then dustbin should > 0'
279           stop
280         endif
281
282         READ(99,fmt='(a)') ch1
283         READ(99,*) sedimentation
284         WRITE(*,8000) ch1,sedimentation
285
286         READ(99,fmt='(a)') ch1
287         READ(99,*) iceparty
288         WRITE(*,8000) ch1,iceparty
289
290         READ(99,fmt='(a)') ch1
291         READ(99,*) activice
292         WRITE(*,8000) ch1,activice
293
294c Test of incompatibility:
295c if activice is used, then iceparty should be used too
296
297         if (activice.and..not.iceparty) then
298           print*,'if activice is used, iceparty should be used too'
299           stop
300         endif
301
302         READ(99,fmt='(a)') ch1
303         READ(99,*) water
304         WRITE(*,8000) ch1,water
305
306c Test of incompatibility:
307c if iceparty is used, then water should be used too
308
309         if (.not.water) then
310            iceparty = .false.
311         endif
312
313         READ(99,fmt='(a)') ch1
314         READ(99,*) caps
315         WRITE(*,8000) ch1,caps
316
317         READ(99,fmt='(a)') ch1
318         READ(99,*) photochem
319         WRITE(*,8000) ch1,photochem
320
321         READ(99,*)
322         READ(99,*)
323
324c THERMOSPHERE
325
326         READ(99,fmt='(a)') ch1
327         READ(99,'(l1)') callthermos
328         WRITE(*,8000) ch1,callthermos
329
330         READ(99,fmt='(a)') ch1
331         READ(99,'(l1)') thermoswater
332         WRITE(*,8000) ch1,thermoswater
333
334         READ(99,fmt='(a)') ch1
335         READ(99,'(l1)') callconduct
336         WRITE(*,8000) ch1,callconduct
337
338         READ(99,fmt='(a)') ch1
339         READ(99,'(l1)') calleuv
340         WRITE(*,8000) ch1,calleuv
341
342         READ(99,fmt='(a)') ch1
343         READ(99,'(l1)') callmolvis
344         WRITE(*,8000) ch1,callmolvis
345
346         READ(99,fmt='(a)') ch1
347         READ(99,'(l1)') callmoldiff
348         WRITE(*,8000) ch1,callmoldiff
349
350         READ(99,fmt='(a)') ch1
351         READ(99,'(l1)') thermochem
352         WRITE(*,8000) ch1,thermochem
353
354         READ(99,fmt='(a)') ch1
355         READ(99,*) solarcondate
356         WRITE(*,*) ch1,solarcondate
357
358         if (.not.callthermos) then
359                 thermoswater=.false.         
360                 callconduct=.false.         
361                 calleuv=.false.         
362                 callmolvis=.false.         
363                 callmoldiff=.false.         
364                 thermochem=.false.         
365        end if
366
367c Test of incompatibility:
368c if photochem is used, then water should be used too
369
370         if (photochem.and..not.water) then
371           print*,'if photochem is used, water should be used too'
372           stop
373         endif
374
375c if callthermos is used, then thermoswater should be used too
376c (if water not used already)
377
378         if (callthermos .and. .not.water) then
379           if (callthermos .and. .not.thermoswater) then
380             print*,'if callthermos is used, water or thermoswater
381     &               should be used too'
382             stop
383           endif
384         endif
385
386         PRINT*,'--------------------------------------------'
387         PRINT*
388         PRINT*
389      ELSE
390         write(*,*)
391         write(*,*) 'Cannot read file callphys.def. Is it here ?'
392         stop
393      ENDIF
394      CLOSE(99)
395
3968000  FORMAT(t5,a12,l8)
3978001  FORMAT(t5,a12,i8)
3988002  FORMAT(t5,a12,f8.1)
399
400      PRINT*
401      PRINT*,'daysec',daysec
402      PRINT*
403      PRINT*,'The radiative transfer is computed:'
404      PRINT*,'           each ',iradia,' physical time-step'
405      PRINT*,'        or each ',iradia*dtphys,' seconds'
406      PRINT*
407c --------------------------------------------------------------
408c  Managing the Longwave radiative transfer
409c --------------------------------------------------------------
410
411c     In most cases, the run just use the following values :
412c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
413      callemis=.true.     
414c     ilwd=10*int(daysec/dtphys) ! bug before 22/10/01       
415      ilwd=10*int(daysec/dtphys)
416      ilwn=2               
417      linear=.true.       
418      ncouche=3
419      alphan=0.4
420      ilwb=2
421      semi=0
422
423c     BUT people working hard on the LW may want to read them in 'radia.def'
424c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
425
426      OPEN(99,file='radia.def',status='old',form='formatted'
427     .     ,iostat=ierr)
428      IF(ierr.EQ.0) THEN
429         write(*,*) 'Reading radia.def !!!'
430         READ(99,fmt='(a)') ch1
431         READ(99,*) callemis
432         WRITE(*,8000) ch1,callemis
433
434         READ(99,fmt='(a)') ch1
435         READ(99,*) iradia
436         WRITE(*,8001) ch1,iradia
437
438         READ(99,fmt='(a)') ch1
439         READ(99,*) ilwd
440         WRITE(*,8001) ch1,ilwd
441
442         READ(99,fmt='(a)') ch1
443         READ(99,*) ilwn
444         WRITE(*,8001) ch1,ilwn
445
446         READ(99,fmt='(a)') ch1
447         READ(99,*) linear
448         WRITE(*,8000) ch1,linear
449
450         READ(99,fmt='(a)') ch1
451         READ(99,*) ncouche
452         WRITE(*,8001) ch1,ncouche
453
454         READ(99,fmt='(a)') ch1
455         READ(99,*) alphan
456         WRITE(*,*) ch1,alphan
457
458         READ(99,fmt='(a)') ch1
459         READ(99,*) ilwb
460         WRITE(*,8001) ch1,ilwb
461
462
463         READ(99,fmt='(a)') ch1
464         READ(99,'(l1)') callg2d
465         WRITE(*,8000) ch1,callg2d
466
467         READ(99,fmt='(a)') ch1
468         READ(99,*) semi
469         WRITE(*,*) ch1,semi
470      end if
471
472c-----------------------------------------------------------------------
473c     Some more initialization:
474c     ------------------------
475
476      CALL SCOPY(ngrid,plon,1,long,1)
477      CALL SCOPY(ngrid,plat,1,lati,1)
478      CALL SCOPY(ngrid,parea,1,area,1)
479      totarea=SSUM(ngridmx,area,1)
480
481      DO ig=1,ngrid
482         sinlat(ig)=sin(plat(ig))
483         coslat(ig)=cos(plat(ig))
484         sinlon(ig)=sin(plon(ig))
485         coslon(ig)=cos(plon(ig))
486      ENDDO
487
488      pi=2.*asin(1.)
489
490c     managing the tracers, and tests:
491c     -------------------------------
492
493      if(tracer) then
494
495c          when photochem is used, nqchem_min is the rank
496c          of the first chemical species
497
498       nqchem_min = 1
499       if (photochem .or. callthermos) then
500         chem = .true.
501        if (doubleq) then
502          nqchem_min = 3
503        else
504          nqchem_min = dustbin+1
505        end if
506       end if
507
508       if (water .or. thermoswater) h2o = .true.
509
510c          TESTS
511
512       print*,'TRACERS:'
513
514       if ((doubleq).and.(h2o).and.
515     $     (chem).and.(iceparty)) then
516         print*,' 1: dust ; 2: dust (doubleq)'
517         print*,' 3 to ',nqmx-2,': chemistry'
518         print*,nqmx-1,': water ice ; ',nqmx,': water vapor'
519       endif
520
521       if ((doubleq).and.(h2o).and.
522     $     (chem).and..not.(iceparty)) then
523         print*,' 1: dust ; 2: dust (doubleq)'
524         print*,' 3 to ',nqmx-1,': chemistry'
525         print*,nqmx,': water vapor'
526       endif
527
528       if ((doubleq).and.(h2o).and.
529     $     .not.(chem).and.(iceparty)) then
530         print*,' 1: dust ; 2: dust (doubleq)'
531         print*,nqmx-1,': water ice ; ',nqmx,': water vapor'
532         if (nqmx.ne.4) then
533           print*,'nqmx should be 4 with these options.'
534                   print*,'(or check callphys.def)'
535           stop
536         endif
537       endif
538
539       if ((doubleq).and.(h2o).and.
540     $     .not.(chem).and..not.(iceparty)) then
541         print*,' 1: dust ; 2: dust (doubleq)'
542         print*,nqmx,': water vapor'
543         if (nqmx.ne.3) then
544           print*,'nqmx should be 3 with these options...'
545                   print*,'(or check callphys.def)'
546           stop
547         endif
548       endif
549
550       if ((doubleq).and..not.(h2o)) then
551         print*,' 1: dust ; 2: dust (doubleq)'
552         if (nqmx.ne.2) then
553           print*,'nqmx should be 2 with these options...'
554                   print*,'(or check callphys.def)'
555           stop
556         endif
557       endif
558
559       if (.not.(doubleq).and.(h2o).and.
560     $     (chem).and.(iceparty)) then
561         if (dustbin.gt.0) then
562           print*,' 1 to ',dustbin,': dust bins'
563         endif
564         print*,nqchem_min,' to ',nqmx-2,': chemistry'
565         print*,nqmx-1,': water ice ; ',nqmx,': water vapor'
566       endif
567       if (.not.(doubleq).and.(h2o).and.
568     $     (chem).and..not.(iceparty)) then
569         if (dustbin.gt.0) then
570           print*,' 1 to ',dustbin,': dust bins'
571         endif
572         print*,nqchem_min,' to ',nqmx-1,': chemistry'
573         print*,nqmx,': water vapor'
574       endif
575       if (.not.(doubleq).and.(h2o).and.
576     $     .not.(chem).and.(iceparty)) then
577         if (dustbin.gt.0) then
578           print*,' 1 to ',dustbin,': dust bins'
579         endif
580         print*,nqmx-1,': water ice ; ',nqmx,': water vapor'
581         if (nqmx.ne.(dustbin+2)) then
582           print*,'nqmx should be ',(dustbin+2),
583     $            ' with these options...'
584                   print*,'(or check callphys.def)'
585           stop
586         endif
587       endif
588       if (.not.(doubleq).and.(h2o).and.
589     $     .not.(chem).and..not.(iceparty)) then
590         if (dustbin.gt.0) then
591           print*,' 1 to ',dustbin,': dust bins'
592         endif
593         print*,nqmx,': water vapor'
594         if (nqmx.ne.(dustbin+1)) then
595           print*,'nqmx should be ',(dustbin+1),
596     $            ' with these options...'
597                   print*,'(or check callphys.def)'
598           stop
599         endif
600       endif
601       if (.not.(doubleq).and..not.(h2o)) then
602         if (dustbin.gt.0) then
603           print*,' 1 to ',dustbin,': dust bins'
604           if (nqmx.ne.dustbin) then
605             print*,'nqmx should be ',dustbin,
606     $              ' with these options...'
607             print*,'(or check callphys.def)'
608             stop
609           endif
610         else
611           print*,'dustbin=',dustbin,
612     $            ': tracer should be F with these options...'
613     $           ,'UNLESS you just want to move tracers around '
614         endif
615       endif
616
617      endif
618
619
620      RETURN
621      END
Note: See TracBrowser for help on using the repository browser.