source: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/defrun_new.F @ 2943

Last change on this file since 2943 was 2167, checked in by emillour, 5 years ago

Mars GCM:
Big cleanup: remove obsolete compilation scripts (makgcm_*) and dynamical core
(since the one in LMDZ.COMMON should be used instead).
EM

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