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

Last change on this file since 1422 was 1422, checked in by milmd, 10 years ago

In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).

File size: 19.7 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, only: getin
40      use control_mod, only: ndynstep, day_step, iperiod, iconser,
41     &                       idissip, iphysiq, anneeref, ecritphy,
42     &                       ecritstart, timestart, nday_r
43      USE logic_mod, ONLY: hybrid,purmats,physic,grireg,fxyhypb,ysinus
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 "sponge.h"
53#include "iniprint.h"
54c
55c   arguments:
56c   ---------
57      LOGICAL  etatinit ! should be .false. for a call from gcm.F
58                        ! and .true. for a call from newstart.F
59      INTEGER  tapedef  ! unit number to assign to 'run.def' file
60c
61c   local variables:
62c   ---------------
63
64      CHARACTER ch1*72,ch2*72,ch3*72,ch4*8 ! to store various strings
65      integer tapeerr ! unit number for error message
66      parameter (tapeerr=0)
67
68c     REAL clonn,clatt,alphaxx,alphayy
69c     LOGICAL  fxyhypbb
70      INTEGER ierr
71      REAL clonn,clatt,grossismxx,grossismyy
72      REAL dzoomxx,dzoomyy,tauxx,tauyy,temp
73      LOGICAL  fxyhypbb, ysinuss
74
75
76c   initialisations:
77c   ----------------
78 
79!      lunout=6
80
81c-----------------------------------------------------------------------
82c  Parametres de controle du run:
83c-----------------------------------------------------------------------
84
85
86!Initialisation des parametres "terrestres", qui ne concernent pas
87!le modele martien et ne sont donc plus lues dans "run.def"
88
89      anneeref=0
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        WRITE(lunout,*) ""
107        WRITE(lunout,*) "Number of days to run:"
108        nday_r=1 ! default value
109        call getin("nday",nday_r)
110        WRITE(lunout,*)" nday = ",nday_r
111       
112        ndynstep=-9999 ! default value
113        call getin("ndynstep",ndynstep)
114        if (ndynstep .gt. 0) then
115          WRITE(lunout,*) ""
116          WRITE(lunout,*) "Number of dynamical steps to run:"
117          WRITE(lunout,*) " ndynstep = ",ndynstep
118          WRITE(lunout,*) " nday value is now discarded "
119        endif
120
121        WRITE(lunout,*) ""
122        WRITE(lunout,*) "Number of dynamical steps per day:",
123     & "(should be a multiple of iperiod)"
124        day_step=960 ! default value
125        call getin("day_step",day_step)
126        WRITE(lunout,*)" day_step = ",day_step
127
128        WRITE(lunout,*) ""
129        WRITE(lunout,*) "periode pour le pas Matsuno (en pas)"
130        iperiod=5 ! default value
131        call getin("iperiod",iperiod)
132        ! verify that day_step is a multiple of iperiod
133        if (((1.*day_step)/iperiod).ne.(day_step/iperiod)) then
134          write(lunout,*)" Error! iperiod must be such that",
135     &    " day_step is a multiple of iperiod, but iperiod=",
136     &    iperiod," and day_step=",day_step
137        else
138          WRITE(lunout,*)" iperiod = ",iperiod
139        endif
140       
141        WRITE(lunout,*) ""
142        WRITE(lunout,*) "periode de sortie des variables de ",
143     &  "controle (en pas)"
144        iconser=120 ! default value
145        call getin("iconser",iconser)
146        WRITE(lunout,*)" iconser = ",iconser
147
148        WRITE(lunout,*) ""
149        WRITE(lunout,*) "periode de la dissipation (en pas)"
150        idissip=5 ! default value
151        call getin("idissip",idissip)
152        WRITE(lunout,*)" idissip = ",idissip
153
154ccc  ....   P. Le Van , modif le 29/04/97 .pour la dissipation  ...
155ccc
156        WRITE(lunout,*) ""
157        WRITE(lunout,*) "choix de l'operateur de dissipation ",
158     & "(star ou  non star )"
159        lstardis=.true. ! default value
160        call getin("lstardis",lstardis)
161        WRITE(lunout,*)" lstardis = ",lstardis
162
163        WRITE(lunout,*) ""
164        WRITE(lunout,*) "avec ou sans coordonnee hybrides"
165        hybrid=.true. ! default value
166        call getin("hybrid",hybrid)
167        WRITE(lunout,*)" hybrid = ",hybrid
168
169        WRITE(lunout,*) ""
170        WRITE(lunout,*) "nombre d'iterations de l'operateur de ",
171     & "dissipation   gradiv "
172        nitergdiv=1 ! default value
173        call getin("nitergdiv",nitergdiv)
174        WRITE(lunout,*)" nitergdiv = ",nitergdiv
175
176        WRITE(lunout,*) ""
177        WRITE(lunout,*) "nombre d'iterations de l'operateur de ",
178     & "dissipation  nxgradrot"
179        nitergrot=2 ! default value
180        call getin("nitergrot",nitergrot)
181        WRITE(lunout,*)" nitergrot = ",nitergrot
182
183        WRITE(lunout,*) ""
184        WRITE(lunout,*) "nombre d'iterations de l'operateur de ",
185     & "dissipation  divgrad"
186        niterh=2 ! default value
187        call getin("niterh",niterh)
188        WRITE(lunout,*)" niterh = ",niterh
189
190        WRITE(lunout,*) ""
191        WRITE(lunout,*) "temps de dissipation des plus petites ",
192     & "long.d ondes pour u,v (gradiv)"
193        tetagdiv=4000. ! default value
194        call getin("tetagdiv",tetagdiv)
195        WRITE(lunout,*)" tetagdiv = ",tetagdiv
196
197        WRITE(lunout,*) ""
198        WRITE(lunout,*) "temps de dissipation des plus petites ",
199     & "long.d ondes pour u,v(nxgradrot)"
200        tetagrot=5000. ! default value
201        call getin("tetagrot",tetagrot)
202        WRITE(lunout,*)" tetagrot = ",tetagrot
203
204        WRITE(lunout,*) ""
205        WRITE(lunout,*) "temps de dissipation des plus petites ",
206     & "long.d ondes pour  h ( divgrad)"
207        tetatemp=5000. ! default value
208        call getin("tetatemp",tetatemp)
209        WRITE(lunout,*)" tetatemp = ",tetatemp
210
211        WRITE(lunout,*) ""
212        WRITE(lunout,*) "coefficient pour gamdissip"
213        coefdis=0. ! default value
214        call getin("coefdis",coefdis)
215        WRITE(lunout,*)" coefdis = ",coefdis
216c
217c    ...............................................................
218
219        WRITE(lunout,*) ""
220        WRITE(lunout,*) "choix du shema d'integration temporelle ",
221     & "(true==Matsuno ou false==Matsuno-leapfrog)"
222        purmats=.false. ! default value
223        call getin("purmats",purmats)
224        WRITE(lunout,*)" purmats = ",purmats
225
226        WRITE(lunout,*) ""
227        WRITE(lunout,*) "avec ou sans physique"
228        physic=.true. ! default value
229        call getin("physic",physic)
230        WRITE(lunout,*)" physic = ",physic
231
232        WRITE(lunout,*) ""
233        WRITE(lunout,*) "periode de la physique (en pas)"
234        iphysiq=20 ! default value
235        call getin("iphysiq",iphysiq)
236        ! verify that day_step is a multiple of iphysiq
237        if (((1.*day_step)/iphysiq).ne.(day_step/iphysiq)) then
238          write(lunout,*)" Error! iphysiq must be such that",
239     &    " day_step is a multiple of iphysiq, but iphysiq=",
240     &    iphysiq," and day_step=",day_step
241        else
242          WRITE(lunout,*)" iphysiq = ",iphysiq
243        endif
244
245        WRITE(lunout,*) ""
246        WRITE(lunout,*) "choix d'une grille reguliere"
247        grireg=.true.
248        call getin("grireg",grireg)
249        WRITE(lunout,*)" grireg = ",grireg
250
251ccc   .... P.Le Van, ajout le 03/01/96 pour l'ecriture phys ...
252c
253        WRITE(lunout,*) ""
254        WRITE(lunout,*) "frequence (en pas) de l'ecriture ",
255     & "du fichier diagfi.nc"
256        ecritphy=240
257        call getin("ecritphy",ecritphy)
258        ! verify that ecriphy is indeed a multiple of iphysiq
259        if (((1.*ecritphy)/iphysiq).ne.(ecritphy/iphysiq)) then
260          write(lunout,*)" Error! ecritphy must be a multiple",
261     &    " of iphysiq, but ecritphy=",ecritphy," and iphysiq=",
262     &    iphysiq
263        else
264          WRITE(lunout,*)" ecritphy = ",ecritphy
265        endif
266       
267ccc   .... T.Navarro, read start time (06/2013) ...
268c
269        WRITE(lunout,*) ""
270        WRITE(lunout,*) "date de debut dans le fichier start.nc"
271        timestart=-9999
272        call getin("timestart",timestart)
273        WRITE(lunout,*)" timestart = ",timestart
274
275ccc   .... T.Navarro, start output (01/2013) ...
276c
277        WRITE(lunout,*) ""
278        WRITE(lunout,*) "frequence (en pas) de l'ecriture ",
279     & "du fichier start.nc"
280        ecritstart=0
281        call getin("ecritstart",ecritstart)
282        ! verify that ecritstart is indeed a multiple of iphysiq
283        if ( ((real(ecritstart))/iphysiq)
284     &       .ne.(real(ecritstart)/iphysiq) ) then
285          write(lunout,*)" Error! ecritstart must be a multiple",
286     &    " of iphysiq, but ecritstart=",ecritstart," and iphysiq=",
287     &    iphysiq
288        else
289          WRITE(lunout,*)" ecritstart = ",ecritstart
290        endif
291
292ccc  ....   P. Le Van , ajout  le 7/03/95 .pour le zoom ...
293c     .........   (  modif  le 17/04/96 )   .........
294c
295        if (.not.etatinit ) then
296
297           clonn=63.
298           call getin("clon",clonn)
299           
300           IF( ABS(clon - clonn).GE. 0.001 )  THEN
301             PRINT *,' La valeur de clon passee par run.def est '
302     *       ,'differente de celle lue sur le fichier start '
303             STOP
304           ENDIF
305c
306c
307           clatt=0.
308           call getin("clat",clatt)
309 
310           IF( ABS(clat - clatt).GE. 0.001 )  THEN
311             PRINT *,' La valeur de clat passee par run.def est '
312     *       ,'differente de celle lue sur le fichier start '
313             STOP
314           ENDIF
315
316           grossismxx=1.0
317           call getin("grossismx",grossismxx)
318
319           if(grossismxx.eq.0) then 
320             write(*,*)
321             write(*,*)'ERREUR : dans run.def, grossismx =0'
322             write(*,*)'Attention a ne pas utiliser une version de'
323             write(*,*)'run.def avant le nouveau zoom LMDZ2.3 (06/2000)'
324             write(*,*)'(Il faut ajouter grossismx,dzoomx,etc... a la'
325             write(*,*)'place de alphax, alphay. cf. dyn3d). '
326             write(*,*)
327             stop
328           end if
329
330           IF( ABS(grossismx - grossismxx).GE. 0.001 )  THEN
331             PRINT *,' La valeur de grossismx passee par run.def est '
332     *       ,'differente de celle lue sur le fichier  start =',
333     *        grossismx
334             if (grossismx.eq.0) then
335                  write(*,*) 'OK,Normal : c est un vieux start'
336     *             , 'd avant le nouveau zoom LMDZ2.3 (06/2000)'
337                 grossismx=grossismxx
338             else
339                   STOP
340             endif
341           ENDIF
342
343           grossismyy=1.0
344           call getin("grossismy",grossismyy)
345
346           IF( ABS(grossismy - grossismyy).GE. 0.001 )  THEN
347             PRINT *,' La valeur de grossismy passee par run.def est '
348     *       ,'differente de celle lue sur le fichier  start =',
349     *        grossismy
350             if (grossismy.eq.0) then
351                  write(*,*) 'OK,Normal : c est un vieux start'
352     *             , 'd avant le nouveau zoom LMDZ2.3 (06/2000)'
353                 grossismy=grossismyy
354             else
355                   STOP
356             endif
357           ENDIF
358
359
360           IF( grossismx.LT.1. )  THEN
361             PRINT *,' ***  ATTENTION !! grossismx < 1 .   *** '
362             STOP
363           ELSE
364             alphax = 1. - 1./ grossismx
365           ENDIF
366
367           IF( grossismy.LT.1. )  THEN
368             PRINT *,' ***  ATTENTION !! grossismy < 1 .   *** '
369             STOP
370           ELSE
371             alphay = 1. - 1./ grossismy
372           ENDIF
373
374           PRINT *,' '
375           PRINT *,' --> In defrun: alphax alphay  ',alphax,alphay
376           PRINT *,' '
377c
378           fxyhypbb=.false.
379           call getin("fxyhypbb",fxyhypbb)
380 
381           IF( .NOT.fxyhypb )  THEN
382             IF( fxyhypbb )     THEN
383                PRINT *,' ********  PBS DANS  DEFRUN  ******** '
384                PRINT *,' *** fxyhypb lu sur le fichier start est F ',
385     *          'alors  qu il est  T  sur  run.def  ***'
386                STOP
387             ENDIF
388           ELSE
389             IF( .NOT.fxyhypbb )   THEN
390                PRINT *,' ********  PBS DANS  DEFRUN  ******** '
391                PRINT *,' ***  fxyhypb lu sur le fichier start est T ',
392     *         'alors  qu il est  F  sur  run.def  ****  '
393                STOP
394             ENDIF
395           ENDIF
396           dzoomxx=0.0
397           call getin("dzoomx",dzoomxx)
398
399           IF( fxyhypb )  THEN
400              IF( ABS(dzoomx - dzoomxx).GE. 0.001 )  THEN
401                PRINT *,' La valeur de dzoomx passee par run.def est '
402     *          ,'differente de celle lue sur le fichier  start '
403                STOP
404              ENDIF
405           ENDIF
406
407           dzoomyy=0.0
408           call getin("dzoomy",dzoomyy)
409
410           IF( fxyhypb )  THEN
411              IF( ABS(dzoomy - dzoomyy).GE. 0.001 )  THEN
412                PRINT *,' La valeur de dzoomy passee par run.def est '
413     *          ,'differente de celle lue sur le fichier  start '
414                STOP
415              ENDIF
416           ENDIF
417
418           tauxx=2.0
419           call getin("taux",tauxx)
420
421           tauyy=2.0
422           call getin("tauy",tauyy)
423
424           IF( fxyhypb )  THEN
425              IF( ABS(taux - tauxx).GE. 0.001 )  THEN
426                WRITE(6,*)' La valeur de taux passee par run.def est',
427     *             'differente de celle lue sur le fichier  start '
428                CALL ABORT
429              ENDIF
430
431              IF( ABS(tauy - tauyy).GE. 0.001 )  THEN
432                WRITE(6,*)' La valeur de tauy passee par run.def est',
433     *          'differente de celle lue sur le fichier  start '
434                CALL ABORT
435              ENDIF
436           ENDIF
437 
438        ELSE    ! Below, case when etainit=.true.
439
440           WRITE(lunout,*) ""
441           WRITE(lunout,*) "longitude en degres du centre du zoom"
442           clon=63. ! default value
443           call getin("clon",clon)
444           WRITE(lunout,*)" clon = ",clon
445           
446c
447           WRITE(lunout,*) ""
448           WRITE(lunout,*) "latitude en degres du centre du zoom "
449           clat=0. ! default value
450           call getin("clat",clat)
451           WRITE(lunout,*)" clat = ",clat
452
453           WRITE(lunout,*) ""
454           WRITE(lunout,*) "facteur de grossissement du zoom,",
455     & " selon longitude"
456           grossismx=1.0 ! default value
457           call getin("grossismx",grossismx)
458           WRITE(lunout,*)" grossismx = ",grossismx
459
460           WRITE(lunout,*) ""
461           WRITE(lunout,*) "facteur de grossissement du zoom ,",
462     & " selon latitude"
463           grossismy=1.0 ! default value
464           call getin("grossismy",grossismy)
465           WRITE(lunout,*)" grossismy = ",grossismy
466
467           IF( grossismx.LT.1. )  THEN
468            PRINT *,' ***  ATTENTION !! grossismx < 1 .   *** '
469            STOP
470           ELSE
471             alphax = 1. - 1./ grossismx
472           ENDIF
473
474           IF( grossismy.LT.1. )  THEN
475             PRINT *,' ***  ATTENTION !! grossismy < 1 .   *** '
476             STOP
477           ELSE
478             alphay = 1. - 1./ grossismy
479           ENDIF
480
481           PRINT *,' Defrun  alphax alphay  ',alphax,alphay
482c
483           WRITE(lunout,*) ""
484           WRITE(lunout,*) "Fonction  f(y)  hyperbolique  si = .true.",
485     &  ", sinon  sinusoidale"
486           fxyhypb=.false. ! default value
487           call getin("fxyhypb",fxyhypb)
488           WRITE(lunout,*)" fxyhypb = ",fxyhypb
489
490           WRITE(lunout,*) ""
491           WRITE(lunout,*) "extension en longitude de la zone du zoom",
492     & " (fraction de la zone totale)"
493           dzoomx=0. ! default value
494           call getin("dzoomx",dzoomx)
495           WRITE(lunout,*)" dzoomx = ",dzoomx
496
497           WRITE(lunout,*) ""
498           WRITE(lunout,*) "extension en latitude de la zone du zoom",
499     & " (fraction de la zone totale)"
500           dzoomy=0. ! default value
501           call getin("dzoomy",dzoomy)
502           WRITE(lunout,*)" dzoomy = ",dzoomy
503
504           WRITE(lunout,*) ""
505           WRITE(lunout,*) "raideur du zoom en  X"
506           taux=2. ! default value
507           call getin("taux",taux)
508           WRITE(lunout,*)" taux = ",taux
509
510           WRITE(lunout,*) ""
511           WRITE(lunout,*) "raideur du zoom en  Y"
512           tauy=2.0 ! default value
513           call getin("tauy",tauy)
514           WRITE(lunout,*)" tauy = ",tauy
515
516        END IF ! of if (.not.etatinit )
517
518        WRITE(lunout,*) ""
519        WRITE(lunout,*) "Avec sponge layer"
520        callsponge=.true. ! default value
521        call getin("callsponge",callsponge)
522        WRITE(lunout,*)" callsponge = ",callsponge
523
524        WRITE(lunout,*) ""
525        WRITE(lunout,*) "Sponge: number of layers over which",
526     &                    " sponge extends"
527        nsponge=3 ! default value
528        call getin("nsponge",nsponge)
529        WRITE(lunout,*)" nsponge = ",nsponge
530
531        WRITE(lunout,*)""
532        WRITE(lunout,*)"Sponge mode: (forcing is towards ..."
533        WRITE(lunout,*)"  over upper nsponge layers)"
534        WRITE(lunout,*)"  0: (h=hmean,u=v=0)"
535        WRITE(lunout,*)"  1: (h=hmean,u=umean,v=0)"
536        WRITE(lunout,*)"  2: (h=hmean,u=umean,v=vmean)"
537        mode_sponge=2 ! default value
538        call getin("mode_sponge",mode_sponge)
539        WRITE(lunout,*)" mode_sponge = ",mode_sponge
540
541        WRITE(lunout,*) ""
542        WRITE(lunout,*) "Sponge: characteristice time scale tetasponge"
543        WRITE(lunout,*) "(seconds) at topmost layer (time scale then "
544        WRITE(lunout,*) " doubles with decreasing layer index)."
545        tetasponge=30000.0
546        call getin("tetasponge",tetasponge)
547        WRITE(lunout,*)" tetasponge = ",tetasponge
548
549
550      WRITE(lunout,*) '-----------------------------------------------'
551      WRITE(lunout,*) ' '
552      WRITE(lunout,*) ' '
553c
554
555c       Unlike on Earth (cf LMDZ2.2) , always a regular grid on Mars :
556        ysinus = .false. !Mars Mettre a jour
557
558
559      WRITE(lunout,*) '-----------------------------------------------'
560      WRITE(lunout,*) ' '
561      WRITE(lunout,*) ' '
562cc
563      ELSE
564        write(tapeerr,*) ' WHERE IS run.def ? WE NEED IT !!!!!!!!!!!!!!'
565        stop
566      ENDIF ! of IF(ierr.eq.0)
567
568c     Test sur le zoom
569
570      if((grossismx.eq.1).and.(grossismy.eq.1)) then 
571c        Pas de zoom :
572         write(lunout,*) 'No zoom ? -> fxyhypb set to False'
573     &   ,'           (It runs better that way)'
574         fxyhypb = .false.
575      else     
576c        Avec Zoom
577         if (.not.fxyhypb) stop 'With zoom, fxyhypb should be set to T
578     &in run.def for this version... -> STOP ! '     
579      end if
580#else
581      write(*,*) "defrun_new should not be used in parallel mode!"
582      stop
583#endif
584! of #ifndef CPP_PARA
585      END
Note: See TracBrowser for help on using the repository browser.