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

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

Update newstart and start2archive programs of LMDZ.GENERIC and LMDZ.MARS to the new organization.

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