source: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/defrun_new.F @ 3533

Last change on this file since 3533 was 3326, checked in by emillour, 8 months ago

Titan PCM:
Some corrections to make start2archve and newstart work.
LR+EM

File size: 18.0 KB
Line 
1      SUBROUTINE defrun_new( tapedef, etatinit )
2#ifndef CPP_PARA
3c
4c-----------------------------------------------------------------------
5c     Auteurs :   L. Fairhead , P. Le Van  .
6c      Modif C. Hourdin F. Forget VERSION MARTIENNE
7c
8c
9c  -------------------------------------------------------------------
10c
11c                    MODIF JUIN 2000 (zoom)
12c       .........     Version  du 29/04/97       ..........
13c
14c   Nouveaux parametres nitergdiv,nitergrot,niterh,tetagdiv,tetagrot,
15c   tetatemp   ajoutes  pour la dissipation   .
16c
17c   Autre parametre ajoute en fin de liste : ** fxyhypb **
18c
19c   Si fxyhypb = .TRUE. , choix de la fonction a derivee tangente hyperb.
20c   Sinon , choix de fxynew  , a derivee sinusoidale  ..
21c
22c   ......  etatinit = . TRUE. si defrun_new  est appele dans NEWSTART
23c   ETAT0_LMD  ou  LIMIT_LMD  pour l'initialisation de start.dat (dic) et
24c   de limit.dat (dic)  ...........
25c   Sinon  etatinit = . FALSE .
26c
27c   Donc etatinit = .F.  si on veut comparer les valeurs de  alphax ,
28c   alphay,clon,clat, fxyhypb  lues sur  le fichier  start  avec
29c   celles passees  par run.def ,  au debut du gcm, apres l'appel a
30c   lectba . 
31c   Ces parametres definissant entre autres la grille et doivent etre
32c   pareils et coherents , sinon il y aura  divergence du gcm .
33c
34c
35c-----------------------------------------------------------------------
36c   Declarations :
37c   --------------
38! to use  'getin'
39      USE ioipsl_getincom
40      use sponge_mod,only: callsponge,nsponge,mode_sponge,tetasponge
41      use control_mod,only: nday, day_step, iperiod, anneeref,
42     &                      iconser, dissip_period, iphysiq, timestart
43      USE logic_mod, ONLY: hybrid,purmats,fxyhypb,ysinus,iflag_phys
44      USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy,
45     .                  alphax,alphay,taux,tauy
46      IMPLICIT NONE
47
48#include "dimensions.h"
49#include "paramet.h"
50!#include "control.h"
51#include "comdissnew.h"
52#include "iniprint.h"
53c
54c   arguments:
55c   ---------
56      LOGICAL  etatinit ! should be .false. for a call from gcm.F
57                        ! and .true. for a call from newstart.F
58      INTEGER  tapedef  ! unit number to assign to 'run.def' file
59c
60c   local variables:
61c   ---------------
62
63      CHARACTER ch1*72,ch2*72,ch3*72,ch4*8 ! to store various strings
64      integer tapeerr ! unit number for error message
65      parameter (tapeerr=0)
66
67c     REAL clonn,clatt,alphaxx,alphayy
68c     LOGICAL  fxyhypbb
69      INTEGER ierr
70      REAL clonn,clatt,grossismxx,grossismyy
71      REAL dzoomxx,dzoomyy,tauxx,tauyy,temp
72      LOGICAL  fxyhypbb, ysinuss
73
74
75c   initialisations:
76c   ----------------
77 
78!      lunout=6
79
80c-----------------------------------------------------------------------
81c  Parametres de controle du run:
82c-----------------------------------------------------------------------
83
84
85!Initialisation des parametres "terrestres", qui ne concernent pas
86!le modele martien et ne sont donc plus lues dans "run.def"
87
88        anneeref=0
89        ! Note: anneref is a common in 'control.h'
90
91      OPEN(tapedef,file='run.def',status='old',form='formatted'
92     .     ,iostat=ierr)
93      CLOSE(tapedef) ! first call to getin will open the file
94
95      !lunout: default unit for the text outputs
96      lunout=6
97      CALL getin('lunout', lunout)
98      IF (lunout /= 5 .and. lunout /= 6) THEN
99        OPEN(UNIT=lunout,FILE='lmdz.out',ACTION='write',
100     &       STATUS='unknown',FORM='formatted')
101      ENDIF
102
103      IF(ierr.EQ.0) THEN ! if file run.def is found
104        WRITE(lunout,*) "DEFRUN_NEW: reading file run.def"
105       
106        timestart=-9999 ! default value; if <0, use last stored time
107        call getin("timestart",timestart)
108
109        WRITE(lunout,*) ""
110        WRITE(lunout,*) "Number of days to run:"
111        nday=1 ! default value
112        call getin("nday",nday)
113        WRITE(lunout,*)" nday = ",nday
114
115        WRITE(lunout,*) ""
116        WRITE(lunout,*) "Number of dynamical steps per day:",
117     & "(should be a multiple of iperiod)"
118        day_step=960 ! default value
119        call getin("day_step",day_step)
120        WRITE(lunout,*)" day_step = ",day_step
121
122        WRITE(lunout,*) ""
123        WRITE(lunout,*) "periode pour le pas Matsuno (en pas)"
124        iperiod=5 ! default value
125        call getin("iperiod",iperiod)
126        WRITE(lunout,*)" iperiod = ",iperiod
127
128        WRITE(lunout,*) ""
129        WRITE(lunout,*) "periode de sortie des variables de ",
130     &  "controle (en pas)"
131        iconser=120 ! default value
132        call getin("iconser",iconser)
133        WRITE(lunout,*)" iconser = ",iconser
134
135        WRITE(lunout,*) ""
136        WRITE(lunout,*) "periode de la dissipation (en pas)"
137        dissip_period=5 ! default value
138        call getin("idissip",dissip_period)
139        ! if there is a "dissip_period" in run.def, it overrides "idissip"
140        call getin("dissip_period",dissip_period)
141        WRITE(lunout,*)" dissip_period = ",dissip_period
142
143ccc  ....   P. Le Van , modif le 29/04/97 .pour la dissipation  ...
144ccc
145        WRITE(lunout,*) ""
146        WRITE(lunout,*) "choix de l'operateur de dissipation ",
147     & "(star ou  non star )"
148        lstardis=.true. ! default value
149        call getin("lstardis",lstardis)
150        WRITE(lunout,*)" lstardis = ",lstardis
151
152        WRITE(lunout,*) ""
153        WRITE(lunout,*) "avec ou sans coordonnee hybrides"
154        hybrid=.true. ! default value
155        call getin("hybrid",hybrid)
156        WRITE(lunout,*)" hybrid = ",hybrid
157
158        WRITE(lunout,*) ""
159        WRITE(lunout,*) "nombre d'iterations de l'operateur de ",
160     & "dissipation   gradiv "
161        nitergdiv=1 ! default value
162        call getin("nitergdiv",nitergdiv)
163        WRITE(lunout,*)" nitergdiv = ",nitergdiv
164
165        WRITE(lunout,*) ""
166        WRITE(lunout,*) "nombre d'iterations de l'operateur de ",
167     & "dissipation  nxgradrot"
168        nitergrot=2 ! default value
169        call getin("nitergrot",nitergrot)
170        WRITE(lunout,*)" nitergrot = ",nitergrot
171
172        WRITE(lunout,*) ""
173        WRITE(lunout,*) "nombre d'iterations de l'operateur de ",
174     & "dissipation  divgrad"
175        niterh=2 ! default value
176        call getin("niterh",niterh)
177        WRITE(lunout,*)" niterh = ",niterh
178
179        WRITE(lunout,*) ""
180        WRITE(lunout,*) "temps de dissipation des plus petites ",
181     & "long.d ondes pour u,v (gradiv)"
182        tetagdiv=4000. ! default value
183        call getin("tetagdiv",tetagdiv)
184        WRITE(lunout,*)" tetagdiv = ",tetagdiv
185
186        WRITE(lunout,*) ""
187        WRITE(lunout,*) "temps de dissipation des plus petites ",
188     & "long.d ondes pour u,v(nxgradrot)"
189        tetagrot=5000. ! default value
190        call getin("tetagrot",tetagrot)
191        WRITE(lunout,*)" tetagrot = ",tetagrot
192
193        WRITE(lunout,*) ""
194        WRITE(lunout,*) "temps de dissipation des plus petites ",
195     & "long.d ondes pour  h ( divgrad)"
196        tetatemp=5000. ! default value
197        call getin("tetatemp",tetatemp)
198        WRITE(lunout,*)" tetatemp = ",tetatemp
199
200        WRITE(lunout,*) ""
201        WRITE(lunout,*) "coefficient pour gamdissip"
202        coefdis=0. ! default value
203        call getin("coefdis",coefdis)
204        WRITE(lunout,*)" coefdis = ",coefdis
205c
206c    ...............................................................
207
208        WRITE(lunout,*) ""
209        WRITE(lunout,*) "choix du shema d'integration temporelle ",
210     & "(true==Matsuno ou false==Matsuno-leapfrog)"
211        purmats=.false. ! default value
212        call getin("purmats",purmats)
213        WRITE(lunout,*)" purmats = ",purmats
214
215        WRITE(lunout,*) ""
216        WRITE(lunout,*) "avec ou sans physique"
217!        physic=.true. ! default value
218!        call getin("physic",physic)
219!        WRITE(lunout,*)" physic = ",physic
220        iflag_phys=1 ! default value
221        call getin("iflag_phys",iflag_phys)
222        WRITE(lunout,*)" iflag_phys = ",iflag_phys
223
224        WRITE(lunout,*) ""
225        WRITE(lunout,*) "periode de la physique (en pas)"
226        iphysiq=20 ! default value
227        call getin("iphysiq",iphysiq)
228        WRITE(lunout,*)" iphysiq = ",iphysiq
229
230!        WRITE(lunout,*) ""
231!        WRITE(lunout,*) "choix d'une grille reguliere"
232!        grireg=.true.
233!        call getin("grireg",grireg)
234!        WRITE(lunout,*)" grireg = ",grireg
235
236ccc   .... P.Le Van, ajout le 03/01/96 pour l'ecriture phys ...
237c
238!        WRITE(lunout,*) ""
239!        WRITE(lunout,*) "frequence (en pas) de l'ecriture ",
240!     & "du fichier diagfi.nc"
241!        ecritphy=240
242!        call getin("ecritphy",ecritphy)
243!        WRITE(lunout,*)" ecritphy = ",ecritphy
244
245ccc  ....   P. Le Van , ajout  le 7/03/95 .pour le zoom ...
246c     .........   (  modif  le 17/04/96 )   .........
247c
248        if (.not.etatinit ) then
249
250           clonn=63.
251           call getin("clon",clonn)
252           
253           IF( ABS(clon - clonn).GE. 0.001 )  THEN
254             PRINT *,' La valeur de clon passee par run.def est '
255     *       ,'differente de celle lue sur le fichier start '
256             STOP
257           ENDIF
258c
259c
260           clatt=0.
261           call getin("clat",clatt)
262 
263           IF( ABS(clat - clatt).GE. 0.001 )  THEN
264             PRINT *,' La valeur de clat passee par run.def est '
265     *       ,'differente de celle lue sur le fichier start '
266             STOP
267           ENDIF
268
269           grossismxx=1.0
270           call getin("grossismx",grossismxx)
271
272           if(grossismxx.eq.0) then 
273             write(*,*)
274             write(*,*)'ERREUR : dans run.def, grossismx =0'
275             write(*,*)'Attention a ne pas utiliser une version de'
276             write(*,*)'run.def avant le nouveau zoom LMDZ2.3 (06/2000)'
277             write(*,*)'(Il faut ajouter grossismx,dzoomx,etc... a la'
278             write(*,*)'place de alphax, alphay. cf. dyn3d). '
279             write(*,*)
280             stop
281           end if
282
283           IF( ABS(grossismx - grossismxx).GE. 0.001 )  THEN
284             PRINT *,' La valeur de grossismx passee par run.def est '
285     *       ,'differente de celle lue sur le fichier  start =',
286     *        grossismx
287             if (grossismx.eq.0) then
288                  write(*,*) 'OK,Normal : c est un vieux start'
289     *             , 'd avant le nouveau zoom LMDZ2.3 (06/2000)'
290                 grossismx=grossismxx
291             else
292                   STOP
293             endif
294           ENDIF
295
296           grossismyy=1.0
297           call getin("grossismy",grossismyy)
298
299           IF( ABS(grossismy - grossismyy).GE. 0.001 )  THEN
300             PRINT *,' La valeur de grossismy passee par run.def est '
301     *       ,'differente de celle lue sur le fichier  start =',
302     *        grossismy
303             if (grossismy.eq.0) then
304                  write(*,*) 'OK,Normal : c est un vieux start'
305     *             , 'd avant le nouveau zoom LMDZ2.3 (06/2000)'
306                 grossismy=grossismyy
307             else
308                   STOP
309             endif
310           ENDIF
311
312
313           IF( grossismx.LT.1. )  THEN
314             PRINT *,' ***  ATTENTION !! grossismx < 1 .   *** '
315             STOP
316           ELSE
317             alphax = 1. - 1./ grossismx
318           ENDIF
319
320           IF( grossismy.LT.1. )  THEN
321             PRINT *,' ***  ATTENTION !! grossismy < 1 .   *** '
322             STOP
323           ELSE
324             alphay = 1. - 1./ grossismy
325           ENDIF
326
327           PRINT *,' '
328           PRINT *,' --> In defrun: alphax alphay  ',alphax,alphay
329           PRINT *,' '
330c
331           fxyhypbb=.false.
332           call getin("fxyhypbb",fxyhypbb)
333 
334           IF( .NOT.fxyhypb )  THEN
335             IF( fxyhypbb )     THEN
336                PRINT *,' ********  PBS DANS  DEFRUN  ******** '
337                PRINT *,' *** fxyhypb lu sur le fichier start est F ',
338     *          'alors  qu il est  T  sur  run.def  ***'
339                STOP
340             ENDIF
341           ELSE
342             IF( .NOT.fxyhypbb )   THEN
343                PRINT *,' ********  PBS DANS  DEFRUN  ******** '
344                PRINT *,' ***  fxyhypb lu sur le fichier start est T ',
345     *         'alors  qu il est  F  sur  run.def  ****  '
346                STOP
347             ENDIF
348           ENDIF
349           dzoomxx=0.0
350           call getin("dzoomx",dzoomxx)
351
352           IF( fxyhypb )  THEN
353              IF( ABS(dzoomx - dzoomxx).GE. 0.001 )  THEN
354                PRINT *,' La valeur de dzoomx passee par run.def est '
355     *          ,'differente de celle lue sur le fichier  start '
356                STOP
357              ENDIF
358           ENDIF
359
360           dzoomyy=0.0
361           call getin("dzoomy",dzoomyy)
362
363           IF( fxyhypb )  THEN
364              IF( ABS(dzoomy - dzoomyy).GE. 0.001 )  THEN
365                PRINT *,' La valeur de dzoomy passee par run.def est '
366     *          ,'differente de celle lue sur le fichier  start '
367                STOP
368              ENDIF
369           ENDIF
370
371           tauxx=2.0
372           call getin("taux",tauxx)
373
374           tauyy=2.0
375           call getin("tauy",tauyy)
376
377           IF( fxyhypb )  THEN
378              IF( ABS(taux - tauxx).GE. 0.001 )  THEN
379                WRITE(6,*)' La valeur de taux passee par run.def est',
380     *             'differente de celle lue sur le fichier  start '
381                CALL ABORT
382              ENDIF
383
384              IF( ABS(tauy - tauyy).GE. 0.001 )  THEN
385                WRITE(6,*)' La valeur de tauy passee par run.def est',
386     *          'differente de celle lue sur le fichier  start '
387                CALL ABORT
388              ENDIF
389           ENDIF
390 
391        ELSE    ! Below, case when etainit=.true.
392
393           WRITE(lunout,*) ""
394           WRITE(lunout,*) "longitude en degres du centre du zoom"
395           clon=63. ! default value
396           call getin("clon",clon)
397           WRITE(lunout,*)" clon = ",clon
398           
399c
400           WRITE(lunout,*) ""
401           WRITE(lunout,*) "latitude en degres du centre du zoom "
402           clat=0. ! default value
403           call getin("clat",clat)
404           WRITE(lunout,*)" clat = ",clat
405
406           WRITE(lunout,*) ""
407           WRITE(lunout,*) "facteur de grossissement du zoom,",
408     & " selon longitude"
409           grossismx=1.0 ! default value
410           call getin("grossismx",grossismx)
411           WRITE(lunout,*)" grossismx = ",grossismx
412
413           WRITE(lunout,*) ""
414           WRITE(lunout,*) "facteur de grossissement du zoom ,",
415     & " selon latitude"
416           grossismy=1.0 ! default value
417           call getin("grossismy",grossismy)
418           WRITE(lunout,*)" grossismy = ",grossismy
419
420           IF( grossismx.LT.1. )  THEN
421            PRINT *,' ***  ATTENTION !! grossismx < 1 .   *** '
422            STOP
423           ELSE
424             alphax = 1. - 1./ grossismx
425           ENDIF
426
427           IF( grossismy.LT.1. )  THEN
428             PRINT *,' ***  ATTENTION !! grossismy < 1 .   *** '
429             STOP
430           ELSE
431             alphay = 1. - 1./ grossismy
432           ENDIF
433
434           PRINT *,' Defrun  alphax alphay  ',alphax,alphay
435c
436           WRITE(lunout,*) ""
437           WRITE(lunout,*) "Fonction  f(y)  hyperbolique  si = .true.",
438     &  ", sinon  sinusoidale"
439           fxyhypb=.false. ! default value
440           call getin("fxyhypb",fxyhypb)
441           WRITE(lunout,*)" fxyhypb = ",fxyhypb
442
443           WRITE(lunout,*) ""
444           WRITE(lunout,*) "extension en longitude de la zone du zoom",
445     & " (fraction de la zone totale)"
446           dzoomx=0. ! default value
447           call getin("dzoomx",dzoomx)
448           WRITE(lunout,*)" dzoomx = ",dzoomx
449
450           WRITE(lunout,*) ""
451           WRITE(lunout,*) "extension en latitude de la zone du zoom",
452     & " (fraction de la zone totale)"
453           dzoomy=0. ! default value
454           call getin("dzoomy",dzoomy)
455           WRITE(lunout,*)" dzoomy = ",dzoomy
456
457           WRITE(lunout,*) ""
458           WRITE(lunout,*) "raideur du zoom en  X"
459           taux=2. ! default value
460           call getin("taux",taux)
461           WRITE(lunout,*)" taux = ",taux
462
463           WRITE(lunout,*) ""
464           WRITE(lunout,*) "raideur du zoom en  Y"
465           tauy=2.0 ! default value
466           call getin("tauy",tauy)
467           WRITE(lunout,*)" tauy = ",tauy
468
469        END IF ! of if (.not.etatinit )
470
471        WRITE(lunout,*) ""
472        WRITE(lunout,*) "Use a sponge layer?"
473        callsponge=.true. ! default value
474        call getin("callsponge",callsponge)
475        WRITE(lunout,*)" callsponge = ",callsponge
476
477        WRITE(lunout,*) ""
478        WRITE(lunout,*) "Sponge: number of layers over which",
479     &                    " sponge extends"
480        nsponge=3 ! default value
481        call getin("nsponge",nsponge)
482        WRITE(lunout,*)" nsponge = ",nsponge
483
484        WRITE(lunout,*)""
485        WRITE(lunout,*)"Sponge mode: (forcing is towards ..."
486        WRITE(lunout,*)"  over upper nsponge layers)"
487        WRITE(lunout,*)"  0: (h=hmean,u=v=0)"
488        WRITE(lunout,*)"  1: (h=hmean,u=umean,v=0)"
489        WRITE(lunout,*)"  2: (h=hmean,u=umean,v=vmean)"
490        mode_sponge=2 ! default value
491        call getin("mode_sponge",mode_sponge)
492        WRITE(lunout,*)" mode_sponge = ",mode_sponge
493
494        WRITE(lunout,*) ""
495        WRITE(lunout,*) "Sponge: characteristic time scale tetasponge"
496        WRITE(lunout,*) "(seconds) at topmost layer (time scale then "
497        WRITE(lunout,*) " doubles with decreasing layer index)."
498        tetasponge=50000.0
499        call getin("tetasponge",tetasponge)
500        WRITE(lunout,*)" tetasponge = ",tetasponge
501
502
503      WRITE(lunout,*) '-----------------------------------------------'
504      WRITE(lunout,*) ' '
505      WRITE(lunout,*) ' '
506c
507
508c       Unlike on Earth (cf LMDZ2.2) , always a regular grid on Mars :
509        ysinus = .false. !Mars Mettre a jour
510
511
512      WRITE(lunout,*) '-----------------------------------------------'
513      WRITE(lunout,*) ' '
514      WRITE(lunout,*) ' '
515cc
516      ELSE
517        write(tapeerr,*) ' WHERE IS run.def ? WE NEED IT !!!!!!!!!!!!!!'
518        stop
519      ENDIF ! of IF(ierr.eq.0)
520
521c     Test sur le zoom
522
523      if((grossismx.eq.1).and.(grossismy.eq.1)) then 
524c        Pas de zoom :
525         write(lunout,*) 'No zoom ? -> fxyhypb set to False'
526     &   ,'           (It runs better that way)'
527         fxyhypb = .false.
528      else     
529c        Avec Zoom
530         if (.not.fxyhypb) stop 'With zoom, fxyhypb should be set to T
531     &in run.def for this version... -> STOP ! '     
532      end if
533#else
534      write(*,*) "defrun_new should not be used in parallel mode!"
535      stop
536#endif
537! of #ifndef CPP_PARA
538      END
Note: See TracBrowser for help on using the repository browser.