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

Last change on this file since 999 was 999, checked in by tnavarro, 11 years ago

Possibility to store multiple initial states in one start/startfi. This is RETROCOMPATIBLE. New option ecrithist in run.def to write data in start/startfi every ecrithist dynamical timestep. New option timestart in run.def to initialize the GCM with the time timestart stored in start

File size: 19.2 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   .... T.Navarro, read start time (06/2013) ...
257c
258        WRITE(tapeout,*) ""
259        WRITE(tapeout,*) "date de debut dans le fichier start.nc"
260        timestart=-9999
261        call getin("timestart",timestart)
262        WRITE(tapeout,*)" timestart = ",timestart
263
264ccc   .... T.Navarro, start output (01/2013) ...
265c
266        WRITE(tapeout,*) ""
267        WRITE(tapeout,*) "frequence (en pas) de l'ecriture ",
268     & "du fichier start.nc"
269        ecritstart=0
270        call getin("ecritstart",ecritstart)
271        ! verify that ecritstart is indeed a multiple of iphysiq
272        if ( ((real(ecritstart))/iphysiq)
273     &       .ne.(real(ecritstart)/iphysiq) ) then
274          write(tapeout,*)" Error! ecritstart must be a multiple",
275     &    " of iphysiq, but ecritstart=",ecritstart," and iphysiq=",
276     &    iphysiq
277        else
278          WRITE(tapeout,*)" ecritstart = ",ecritstart
279        endif
280
281ccc  ....   P. Le Van , ajout  le 7/03/95 .pour le zoom ...
282c     .........   (  modif  le 17/04/96 )   .........
283c
284        if (.not.etatinit ) then
285
286           clonn=63.
287           call getin("clon",clonn)
288           
289           IF( ABS(clon - clonn).GE. 0.001 )  THEN
290             PRINT *,' La valeur de clon passee par run.def est '
291     *       ,'differente de celle lue sur le fichier start '
292             STOP
293           ENDIF
294c
295c
296           clatt=0.
297           call getin("clat",clatt)
298 
299           IF( ABS(clat - clatt).GE. 0.001 )  THEN
300             PRINT *,' La valeur de clat passee par run.def est '
301     *       ,'differente de celle lue sur le fichier start '
302             STOP
303           ENDIF
304
305           grossismxx=1.0
306           call getin("grossismx",grossismxx)
307
308           if(grossismxx.eq.0) then 
309             write(*,*)
310             write(*,*)'ERREUR : dans run.def, grossismx =0'
311             write(*,*)'Attention a ne pas utiliser une version de'
312             write(*,*)'run.def avant le nouveau zoom LMDZ2.3 (06/2000)'
313             write(*,*)'(Il faut ajouter grossismx,dzoomx,etc... a la'
314             write(*,*)'place de alphax, alphay. cf. dyn3d). '
315             write(*,*)
316             stop
317           end if
318
319           IF( ABS(grossismx - grossismxx).GE. 0.001 )  THEN
320             PRINT *,' La valeur de grossismx passee par run.def est '
321     *       ,'differente de celle lue sur le fichier  start =',
322     *        grossismx
323             if (grossismx.eq.0) then
324                  write(*,*) 'OK,Normal : c est un vieux start'
325     *             , 'd avant le nouveau zoom LMDZ2.3 (06/2000)'
326                 grossismx=grossismxx
327             else
328                   STOP
329             endif
330           ENDIF
331
332           grossismyy=1.0
333           call getin("grossismy",grossismyy)
334
335           IF( ABS(grossismy - grossismyy).GE. 0.001 )  THEN
336             PRINT *,' La valeur de grossismy passee par run.def est '
337     *       ,'differente de celle lue sur le fichier  start =',
338     *        grossismy
339             if (grossismy.eq.0) then
340                  write(*,*) 'OK,Normal : c est un vieux start'
341     *             , 'd avant le nouveau zoom LMDZ2.3 (06/2000)'
342                 grossismy=grossismyy
343             else
344                   STOP
345             endif
346           ENDIF
347
348
349           IF( grossismx.LT.1. )  THEN
350             PRINT *,' ***  ATTENTION !! grossismx < 1 .   *** '
351             STOP
352           ELSE
353             alphax = 1. - 1./ grossismx
354           ENDIF
355
356           IF( grossismy.LT.1. )  THEN
357             PRINT *,' ***  ATTENTION !! grossismy < 1 .   *** '
358             STOP
359           ELSE
360             alphay = 1. - 1./ grossismy
361           ENDIF
362
363           PRINT *,' '
364           PRINT *,' --> In defrun: alphax alphay  ',alphax,alphay
365           PRINT *,' '
366c
367           fxyhypbb=.false.
368           call getin("fxyhypbb",fxyhypbb)
369 
370           IF( .NOT.fxyhypb )  THEN
371             IF( fxyhypbb )     THEN
372                PRINT *,' ********  PBS DANS  DEFRUN  ******** '
373                PRINT *,' *** fxyhypb lu sur le fichier start est F ',
374     *          'alors  qu il est  T  sur  run.def  ***'
375                STOP
376             ENDIF
377           ELSE
378             IF( .NOT.fxyhypbb )   THEN
379                PRINT *,' ********  PBS DANS  DEFRUN  ******** '
380                PRINT *,' ***  fxyhypb lu sur le fichier start est T ',
381     *         'alors  qu il est  F  sur  run.def  ****  '
382                STOP
383             ENDIF
384           ENDIF
385           dzoomxx=0.0
386           call getin("dzoomx",dzoomxx)
387
388           IF( fxyhypb )  THEN
389              IF( ABS(dzoomx - dzoomxx).GE. 0.001 )  THEN
390                PRINT *,' La valeur de dzoomx passee par run.def est '
391     *          ,'differente de celle lue sur le fichier  start '
392                STOP
393              ENDIF
394           ENDIF
395
396           dzoomyy=0.0
397           call getin("dzoomy",dzoomyy)
398
399           IF( fxyhypb )  THEN
400              IF( ABS(dzoomy - dzoomyy).GE. 0.001 )  THEN
401                PRINT *,' La valeur de dzoomy passee par run.def est '
402     *          ,'differente de celle lue sur le fichier  start '
403                STOP
404              ENDIF
405           ENDIF
406
407           tauxx=2.0
408           call getin("taux",tauxx)
409
410           tauyy=2.0
411           call getin("tauy",tauyy)
412
413           IF( fxyhypb )  THEN
414              IF( ABS(taux - tauxx).GE. 0.001 )  THEN
415                WRITE(6,*)' La valeur de taux passee par run.def est',
416     *             'differente de celle lue sur le fichier  start '
417                CALL ABORT
418              ENDIF
419
420              IF( ABS(tauy - tauyy).GE. 0.001 )  THEN
421                WRITE(6,*)' La valeur de tauy passee par run.def est',
422     *          'differente de celle lue sur le fichier  start '
423                CALL ABORT
424              ENDIF
425           ENDIF
426 
427        ELSE    ! Below, case when etainit=.true.
428
429           WRITE(tapeout,*) ""
430           WRITE(tapeout,*) "longitude en degres du centre du zoom"
431           clon=63. ! default value
432           call getin("clon",clon)
433           WRITE(tapeout,*)" clon = ",clon
434           
435c
436           WRITE(tapeout,*) ""
437           WRITE(tapeout,*) "latitude en degres du centre du zoom "
438           clat=0. ! default value
439           call getin("clat",clat)
440           WRITE(tapeout,*)" clat = ",clat
441
442           WRITE(tapeout,*) ""
443           WRITE(tapeout,*) "facteur de grossissement du zoom,",
444     & " selon longitude"
445           grossismx=1.0 ! default value
446           call getin("grossismx",grossismx)
447           WRITE(tapeout,*)" grossismx = ",grossismx
448
449           WRITE(tapeout,*) ""
450           WRITE(tapeout,*) "facteur de grossissement du zoom ,",
451     & " selon latitude"
452           grossismy=1.0 ! default value
453           call getin("grossismy",grossismy)
454           WRITE(tapeout,*)" grossismy = ",grossismy
455
456           IF( grossismx.LT.1. )  THEN
457            PRINT *,' ***  ATTENTION !! grossismx < 1 .   *** '
458            STOP
459           ELSE
460             alphax = 1. - 1./ grossismx
461           ENDIF
462
463           IF( grossismy.LT.1. )  THEN
464             PRINT *,' ***  ATTENTION !! grossismy < 1 .   *** '
465             STOP
466           ELSE
467             alphay = 1. - 1./ grossismy
468           ENDIF
469
470           PRINT *,' Defrun  alphax alphay  ',alphax,alphay
471c
472           WRITE(tapeout,*) ""
473           WRITE(tapeout,*) "Fonction  f(y)  hyperbolique  si = .true.",
474     &  ", sinon  sinusoidale"
475           fxyhypb=.false. ! default value
476           call getin("fxyhypb",fxyhypb)
477           WRITE(tapeout,*)" fxyhypb = ",fxyhypb
478
479           WRITE(tapeout,*) ""
480           WRITE(tapeout,*) "extension en longitude de la zone du zoom",
481     & " (fraction de la zone totale)"
482           dzoomx=0. ! default value
483           call getin("dzoomx",dzoomx)
484           WRITE(tapeout,*)" dzoomx = ",dzoomx
485
486           WRITE(tapeout,*) ""
487           WRITE(tapeout,*) "extension en latitude de la zone du zoom",
488     & " (fraction de la zone totale)"
489           dzoomy=0. ! default value
490           call getin("dzoomy",dzoomy)
491           WRITE(tapeout,*)" dzoomy = ",dzoomy
492
493           WRITE(tapeout,*) ""
494           WRITE(tapeout,*) "raideur du zoom en  X"
495           taux=2. ! default value
496           call getin("taux",taux)
497           WRITE(tapeout,*)" taux = ",taux
498
499           WRITE(tapeout,*) ""
500           WRITE(tapeout,*) "raideur du zoom en  Y"
501           tauy=2.0 ! default value
502           call getin("tauy",tauy)
503           WRITE(tapeout,*)" tauy = ",tauy
504
505        END IF ! of if (.not.etatinit )
506
507        WRITE(tapeout,*) ""
508        WRITE(tapeout,*) "Avec sponge layer"
509        callsponge=.true. ! default value
510        call getin("callsponge",callsponge)
511        WRITE(tapeout,*)" callsponge = ",callsponge
512
513        WRITE(tapeout,*) ""
514        WRITE(tapeout,*) "Sponge: number of layers over which",
515     &                    " sponge extends"
516        nsponge=3 ! default value
517        call getin("nsponge",nsponge)
518        WRITE(tapeout,*)" nsponge = ",nsponge
519
520        WRITE(tapeout,*)""
521        WRITE(tapeout,*)"Sponge mode: (forcing is towards ..."
522        WRITE(tapeout,*)"  over upper nsponge layers)"
523        WRITE(tapeout,*)"  0: (h=hmean,u=v=0)"
524        WRITE(tapeout,*)"  1: (h=hmean,u=umean,v=0)"
525        WRITE(tapeout,*)"  2: (h=hmean,u=umean,v=vmean)"
526        mode_sponge=2 ! default value
527        call getin("mode_sponge",mode_sponge)
528        WRITE(tapeout,*)" mode_sponge = ",mode_sponge
529
530        WRITE(tapeout,*) ""
531        WRITE(tapeout,*) "Sponge: characteristice time scale tetasponge"
532        WRITE(tapeout,*) "(seconds) at topmost layer (time scale then "
533        WRITE(tapeout,*) " doubles with decreasing layer index)."
534        tetasponge=30000.0
535        call getin("tetasponge",tetasponge)
536        WRITE(tapeout,*)" tetasponge = ",tetasponge
537
538
539      WRITE(tapeout,*) '-----------------------------------------------'
540      WRITE(tapeout,*) ' '
541      WRITE(tapeout,*) ' '
542c
543
544c       Unlike on Earth (cf LMDZ2.2) , always a regular grid on Mars :
545        ysinus = .false. !Mars Mettre a jour
546
547
548      WRITE(tapeout,*) '-----------------------------------------------'
549      WRITE(tapeout,*) ' '
550      WRITE(tapeout,*) ' '
551cc
552      ELSE
553        write(tapeerr,*) ' WHERE IS run.def ? WE NEED IT !!!!!!!!!!!!!!'
554        stop
555      ENDIF ! of IF(ierr.eq.0)
556
557c     Test sur le zoom
558
559      if((grossismx.eq.1).and.(grossismy.eq.1)) then 
560c        Pas de zoom :
561         write(tapeout,*) 'No zoom ? -> fxyhypb set to False'
562     &   ,'           (It runs better that way)'
563         fxyhypb = .false.
564      else     
565c        Avec Zoom
566         if (.not.fxyhypb) stop 'With zoom, fxyhypb should be set to T
567     &in run.def for this version... -> STOP ! '     
568      end if
569
570      RETURN
571      END
Note: See TracBrowser for help on using the repository browser.