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

Last change on this file since 1130 was 1130, checked in by emillour, 11 years ago

Mars GCM:
Series of changes to enable running in parallel (using LMDZ.COMMON dynamics);
Current LMDZ.MARS can still notheless be compiled and run in serial mode
"as previously".
Summary of main changes:

  • Main programs (newstart, start2archive, xvik) that used to be in dyn3d have been moved to phymars.
  • dyn3d/control.h is now module control_mod.F90
  • rearanged input/outputs routines everywhere to handle serial/MPI cases. physdem.F => phyredem.F90 , phyetat0.F => phyetat0.F90 ; all read/write routines for startfi files are gathered in module iostart.F90
  • added parallelism related routines init_phys_lmdz.F90, comgeomphy.F90, dimphy.F90, iniphysiq.F90, mod_grid_phy_lmdz.F90, mod_phys_lmdz_mpi_data.F90, mod_phys_lmdz_mpi_transfert.F90, mod_phys_lmdz_omp_data.F90, mod_phys_lmdz_omp_transfert.F90, mod_phys_lmdz_para.F90, mod_phys_lmdz_transfert_para.F90 in phymars and mod_const_mpi.F90 in dyn3d (for compliance with parallel case)
  • created generic routines 'planetwide_maxval' and 'planetwide_minval', in module "planetwide_mod", that enable obtaining the min and max of a field over the whole planet.

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