source: trunk/LMDZ.MARS/libf/dyn3d/defrun_new.F @ 937

Last change on this file since 937 was 828, checked in by tnavarro, 12 years ago

new option in run.def: ndynstep to run a number of dynamical time step. More efficient thannday for fractions of days

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