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

Last change on this file since 1399 was 1399, checked in by emillour, 10 years ago

Mars GCM:

  • Follow-up to cleanup in dynamics/physics interface. Add iniprint.h to be in line with what is done in LMDZ.COMMON dynamics.

EM

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