source: LMDZ5/branches/LMDZ5_AR5/libf/dyn3d/gcm.F @ 4116

Last change on this file since 4116 was 1529, checked in by Laurent Fairhead, 13 years ago

Necessary modifications to use the model in an aquaplanet or terraplanet configuration


Modifications nécessaires pour l'utilisation du moèle en configuration aquaplanète
ou terraplanète

Dans la dynamique:
Changement pour le paramètre iflag_phys:

auparavant : 0 pas de physique, 1 phyique, 2 rappel
ici on ajoute iflag_phys>=100 pour les aqua et terra

Du coup, dans leapfrog.F et gcm.F on appelle la physique si
iflag_phys=1.or.iflag_phys>=100
Dans iniacademic, on initialise si iflag_phys>=2 au lieu de =2
Dans gcm.F, on appelle en plus iniaqua (sous une clef CPP_EARTH)
Dans iniacademic, on met ps=108080 pour les aqua et terra pour répondre
à une specification CFMIP.

Dans la physique:
On ajoute phyaqua.F qui contient :
iniaqua : initialise les startphy.nc et limit.nc pour la phyique
zenang_an : calcule un ensolleillement moyen sur l'année en fonction de

la latiude (A. Campoy).

profil_sst : qui calcule différents profils latitudinaux de SSTs.
writelim : écrit le fichier limit.nc

Dans physiq.F
On ajoute la possibilité d'appeller un calcul d'ensoleillement moyen
sur l'année quand solarlong0=1000.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.3 KB
Line 
1!
2! $Id: gcm.F 1529 2011-05-26 15:17:33Z idelkadi $
3!
4c
5c
6      PROGRAM gcm
7
8#ifdef CPP_IOIPSL
9      USE IOIPSL
10#else
11! if not using IOIPSL, we still need to use (a local version of) getin
12      USE ioipsl_getincom
13#endif
14
15      USE filtreg_mod
16      USE infotrac
17      USE control_mod
18
19!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
20! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
21! A nettoyer. On ne veut qu'une ou deux routines d'interface
22! dynamique -> physique pour l'initialisation
23! Ehouarn: for now these only apply to Earth:
24#ifdef CPP_EARTH
25      USE dimphy
26      USE comgeomphy
27      USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
28#endif
29!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
30
31      IMPLICIT NONE
32
33c      ......   Version  du 10/01/98    ..........
34
35c             avec  coordonnees  verticales hybrides
36c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
37
38c=======================================================================
39c
40c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
41c   -------
42c
43c   Objet:
44c   ------
45c
46c   GCM LMD nouvelle grille
47c
48c=======================================================================
49c
50c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
51c      et possibilite d'appeler une fonction f(y)  a derivee tangente
52c      hyperbolique a la  place de la fonction a derivee sinusoidale.
53c  ... Possibilite de choisir le schema pour l'advection de
54c        q  , en modifiant iadv dans traceur.def  (MAF,10/02) .
55c
56c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
57c      Pour Van-Leer iadv=10
58c
59c-----------------------------------------------------------------------
60c   Declarations:
61c   -------------
62
63#include "dimensions.h"
64#include "paramet.h"
65#include "comconst.h"
66#include "comdissnew.h"
67#include "comvert.h"
68#include "comgeom.h"
69#include "logic.h"
70#include "temps.h"
71!!!!!!!!!!!#include "control.h"
72#include "ener.h"
73#include "description.h"
74#include "serre.h"
75!#include "com_io_dyn.h"
76#include "iniprint.h"
77#include "tracstoke.h"
78#ifdef INCA
79! Only INCA needs these informations (from the Earth's physics)
80#include "indicesol.h"
81#endif
82      INTEGER         longcles
83      PARAMETER     ( longcles = 20 )
84      REAL  clesphy0( longcles )
85      SAVE  clesphy0
86
87
88
89      REAL zdtvr
90      INTEGER nbetatmoy, nbetatdem,nbetat
91
92c   variables dynamiques
93      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
94      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
95      REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
96      REAL ps(ip1jmp1)                       ! pression  au sol
97      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
98      REAL pks(ip1jmp1)                      ! exner au  sol
99      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
100      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
101      REAL masse(ip1jmp1,llm)                ! masse d'air
102      REAL phis(ip1jmp1)                     ! geopotentiel au sol
103      REAL phi(ip1jmp1,llm)                  ! geopotentiel
104      REAL w(ip1jmp1,llm)                    ! vitesse verticale
105
106c variables dynamiques intermediaire pour le transport
107
108c   variables pour le fichier histoire
109      REAL dtav      ! intervalle de temps elementaire
110
111      REAL time_0
112
113      LOGICAL lafin
114      INTEGER ij,iq,l,i,j
115
116
117      real time_step, t_wrt, t_ops
118
119      LOGICAL first
120
121      LOGICAL call_iniphys
122      data call_iniphys/.true./
123
124      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
125c+jld variables test conservation energie
126c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
127C     Tendance de la temp. potentiel d (theta)/ d t due a la
128C     tansformation d'energie cinetique en energie thermique
129C     cree par la dissipation
130      REAL dhecdt(ip1jmp1,llm)
131c      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
132c      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
133      CHARACTER (len=15) :: ztit
134c-jld
135
136
137      character (len=80) :: dynhist_file, dynhistave_file
138      character (len=20) :: modname
139      character (len=80) :: abort_message
140! locales pour gestion du temps
141      INTEGER :: an, mois, jour
142      REAL :: heure
143
144
145c-----------------------------------------------------------------------
146c    variables pour l'initialisation de la physique :
147c    ------------------------------------------------
148      INTEGER ngridmx
149      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
150      REAL zcufi(ngridmx),zcvfi(ngridmx)
151      REAL latfi(ngridmx),lonfi(ngridmx)
152      REAL airefi(ngridmx)
153      SAVE latfi, lonfi, airefi
154
155c-----------------------------------------------------------------------
156c   Initialisations:
157c   ----------------
158
159      abort_message = 'last timestep reached'
160      modname = 'gcm'
161      descript = 'Run GCM LMDZ'
162      lafin    = .FALSE.
163      dynhist_file = 'dyn_hist.nc'
164      dynhistave_file = 'dyn_hist_ave.nc'
165
166
167
168c----------------------------------------------------------------------
169c  lecture des fichiers gcm.def ou run.def
170c  ---------------------------------------
171c
172! Ehouarn: dump possibility of using defrun
173!#ifdef CPP_IOIPSL
174      CALL conf_gcm( 99, .TRUE. , clesphy0 )
175!#else
176!      CALL defrun( 99, .TRUE. , clesphy0 )
177!#endif
178
179!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
180! FH 2008/05/02
181! A nettoyer. On ne veut qu'une ou deux routines d'interface
182! dynamique -> physique pour l'initialisation
183! Ehouarn : temporarily (?) keep this only for Earth
184      if (planet_type.eq."earth") then
185#ifdef CPP_EARTH
186      CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
187      call InitComgeomphy
188#endif
189      endif
190!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
191c-----------------------------------------------------------------------
192c   Choix du calendrier
193c   -------------------
194
195c      calend = 'earth_365d'
196
197#ifdef CPP_IOIPSL
198      if (calend == 'earth_360d') then
199        call ioconf_calendar('360d')
200        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
201      else if (calend == 'earth_365d') then
202        call ioconf_calendar('noleap')
203        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
204      else if (calend == 'earth_366d') then
205        call ioconf_calendar('gregorian')
206        write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
207      else
208        abort_message = 'Mauvais choix de calendrier'
209        call abort_gcm(modname,abort_message,1)
210      endif
211#endif
212c-----------------------------------------------------------------------
213
214      IF (config_inca /= 'none') THEN
215#ifdef INCA
216      call init_const_lmdz(nbtr,anneeref,dayref,iphysiq,day_step,nday,
217     $        nbsrf, is_oce,is_sic,is_ter,is_lic)
218      call init_inca_para(iim,jjm+1,klon,1,klon_mpi_para_nb,0)
219#endif
220      END IF
221c
222c
223c------------------------------------
224c   Initialisation partie parallele
225c------------------------------------
226
227c
228c
229c-----------------------------------------------------------------------
230c   Initialisation des traceurs
231c   ---------------------------
232c  Choix du nombre de traceurs et du schema pour l'advection
233c  dans fichier traceur.def, par default ou via INCA
234      call infotrac_init
235
236c Allocation de la tableau q : champs advectes   
237      allocate(q(ip1jmp1,llm,nqtot))
238
239c-----------------------------------------------------------------------
240c   Lecture de l'etat initial :
241c   ---------------------------
242
243c  lecture du fichier start.nc
244      if (read_start) then
245      ! we still need to run iniacademic to initialize some
246      ! constants & fields, if we run the 'newtonian' or 'SW' cases:
247        if (iflag_phys.ne.1) then
248          CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
249        endif
250
251!        if (planet_type.eq."earth") then
252! Load an Earth-format start file
253         CALL dynetat0("start.nc",vcov,ucov,
254     &              teta,q,masse,ps,phis, time_0)
255!        endif ! of if (planet_type.eq."earth")
256       
257c       write(73,*) 'ucov',ucov
258c       write(74,*) 'vcov',vcov
259c       write(75,*) 'teta',teta
260c       write(76,*) 'ps',ps
261c       write(77,*) 'q',q
262
263      endif ! of if (read_start)
264
265      IF (config_inca /= 'none') THEN
266#ifdef INCA
267         call init_inca_dim(klon,llm,iim,jjm,
268     $        rlonu,rlatu,rlonv,rlatv)
269#endif
270      END IF
271
272
273c le cas echeant, creation d un etat initial
274      IF (prt_level > 9) WRITE(lunout,*)
275     .              'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
276      if (.not.read_start) then
277         CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
278      endif
279
280
281c-----------------------------------------------------------------------
282c   Lecture des parametres de controle pour la simulation :
283c   -------------------------------------------------------
284c  on recalcule eventuellement le pas de temps
285
286      IF(MOD(day_step,iperiod).NE.0) THEN
287        abort_message =
288     .  'Il faut choisir un nb de pas par jour multiple de iperiod'
289        call abort_gcm(modname,abort_message,1)
290      ENDIF
291
292      IF(MOD(day_step,iphysiq).NE.0) THEN
293        abort_message =
294     * 'Il faut choisir un nb de pas par jour multiple de iphysiq'
295        call abort_gcm(modname,abort_message,1)
296      ENDIF
297
298      zdtvr    = daysec/REAL(day_step)
299        IF(dtvr.NE.zdtvr) THEN
300         WRITE(lunout,*)
301     .    'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
302        ENDIF
303
304C
305C on remet le calendrier à zero si demande
306c
307      IF (raz_date == 1) THEN
308        annee_ref = anneeref
309        day_ref = dayref
310        day_ini = dayref
311        itau_dyn = 0
312        itau_phy = 0
313        time_0 = 0.
314        write(lunout,*)
315     .   'GCM: On reinitialise a la date lue dans gcm.def'
316      ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN
317        write(lunout,*)
318     .  'GCM: Attention les dates initiales lues dans le fichier'
319        write(lunout,*)
320     .  ' restart ne correspondent pas a celles lues dans '
321        write(lunout,*)' gcm.def'
322        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
323        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
324        write(lunout,*)' Pas de remise a zero'
325      ENDIF
326
327c      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
328c        write(lunout,*)
329c     .  'GCM: Attention les dates initiales lues dans le fichier'
330c        write(lunout,*)
331c     .  ' restart ne correspondent pas a celles lues dans '
332c        write(lunout,*)' gcm.def'
333c        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
334c        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
335c        if (raz_date .ne. 1) then
336c          write(lunout,*)
337c     .    'GCM: On garde les dates du fichier restart'
338c        else
339c          annee_ref = anneeref
340c          day_ref = dayref
341c          day_ini = dayref
342c          itau_dyn = 0
343c          itau_phy = 0
344c          time_0 = 0.
345c          write(lunout,*)
346c     .   'GCM: On reinitialise a la date lue dans gcm.def'
347c        endif
348c      ELSE
349c        raz_date = 0
350c      endif
351
352#ifdef CPP_IOIPSL
353      mois = 1
354      heure = 0.
355      call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
356      jH_ref = jD_ref - int(jD_ref)
357      jD_ref = int(jD_ref)
358
359      call ioconf_startdate(INT(jD_ref), jH_ref)
360
361      write(lunout,*)'DEBUG'
362      write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
363      write(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
364      call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
365      write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
366      write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
367#else
368! Ehouarn: we still need to define JD_ref and JH_ref
369! and since we don't know how many days there are in a year
370! we set JD_ref to 0 (this should be improved ...)
371      jD_ref=0
372      jH_ref=0
373#endif
374
375c  nombre d'etats dans les fichiers demarrage et histoire
376      nbetatdem = nday / iecri
377      nbetatmoy = nday / periodav + 1
378
379      if (iflag_phys.eq.1) then
380      ! these initialisations have already been done (via iniacademic)
381      ! if running in SW or Newtonian mode
382c-----------------------------------------------------------------------
383c   Initialisation des constantes dynamiques :
384c   ------------------------------------------
385        dtvr = zdtvr
386        CALL iniconst
387
388c-----------------------------------------------------------------------
389c   Initialisation de la geometrie :
390c   --------------------------------
391        CALL inigeom
392
393c-----------------------------------------------------------------------
394c   Initialisation du filtre :
395c   --------------------------
396        CALL inifilr
397      endif ! of if (iflag_phys.eq.1)
398c
399c-----------------------------------------------------------------------
400c   Initialisation de la dissipation :
401c   ----------------------------------
402
403      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
404     *                tetagdiv, tetagrot , tetatemp              )
405
406c-----------------------------------------------------------------------
407c   Initialisation de la physique :
408c   -------------------------------
409
410      IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN
411         latfi(1)=rlatu(1)
412         lonfi(1)=0.
413         zcufi(1) = cu(1)
414         zcvfi(1) = cv(1)
415         DO j=2,jjm
416            DO i=1,iim
417               latfi((j-2)*iim+1+i)= rlatu(j)
418               lonfi((j-2)*iim+1+i)= rlonv(i)
419               zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i)
420               zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i)
421            ENDDO
422         ENDDO
423         latfi(ngridmx)= rlatu(jjp1)
424         lonfi(ngridmx)= 0.
425         zcufi(ngridmx) = cu(ip1jm+1)
426         zcvfi(ngridmx) = cv(ip1jm-iim)
427         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
428         WRITE(lunout,*)
429     .       'GCM: WARNING!!! vitesse verticale nulle dans la physique'
430! Earth:
431         if (planet_type.eq."earth") then
432#ifdef CPP_EARTH
433         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys ,
434     ,                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp     )
435#endif
436         endif ! of if (planet_type.eq."earth")
437         call_iniphys=.false.
438      ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
439!#endif
440
441c  numero de stockage pour les fichiers de redemarrage:
442
443c-----------------------------------------------------------------------
444c   Initialisation des I/O :
445c   ------------------------
446
447
448      day_end = day_ini + nday
449      WRITE(lunout,300)day_ini,day_end
450 300  FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
451
452#ifdef CPP_IOIPSL
453      call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
454      write (lunout,301)jour, mois, an
455      call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
456      write (lunout,302)jour, mois, an
457 301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
458 302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
459#endif
460
461#ifdef CPP_EARTH
462! Create start file (startphy.nc) and boundary conditions (limit.nc)
463! for the Earth verstion
464       if (iflag_phys>=100) then
465          call iniaqua(ngridmx,latfi,lonfi,iflag_phys)
466       endif
467#endif
468
469!      if (planet_type.eq."earth") then
470! Write an Earth-format restart file
471
472        CALL dynredem0("restart.nc", day_end, phis)
473!      endif
474
475      ecripar = .TRUE.
476
477#ifdef CPP_IOIPSL
478      time_step = zdtvr
479      if (ok_dyn_ins) then
480        ! initialize output file for instantaneous outputs
481        ! t_ops = iecri * daysec ! do operations every t_ops
482        t_ops =((1.0*iecri)/day_step) * daysec 
483        t_wrt = daysec ! iecri * daysec ! write output every t_wrt
484        CALL inithist(day_ref,annee_ref,time_step,
485     &              t_ops,t_wrt)
486      endif
487
488      IF (ok_dyn_ave) THEN
489        ! initialize output file for averaged outputs
490        t_ops = iperiod * time_step ! do operations every t_ops
491        t_wrt = periodav * daysec   ! write output every t_wrt
492        CALL initdynav(day_ref,annee_ref,time_step,
493     &       t_ops,t_wrt)
494      END IF
495      dtav = iperiod*dtvr/daysec
496#endif
497! #endif of #ifdef CPP_IOIPSL
498
499c  Choix des frequences de stokage pour le offline
500c      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
501c      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
502      istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
503      istphy=istdyn/iphysiq     
504
505
506c
507c-----------------------------------------------------------------------
508c   Integration temporelle du modele :
509c   ----------------------------------
510
511c       write(78,*) 'ucov',ucov
512c       write(78,*) 'vcov',vcov
513c       write(78,*) 'teta',teta
514c       write(78,*) 'ps',ps
515c       write(78,*) 'q',q
516
517
518      CALL leapfrog(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
519     .              time_0)
520
521      END
522
Note: See TracBrowser for help on using the repository browser.