Ignore:
Timestamp:
Sep 11, 2024, 6:03:07 PM (12 days ago)
Author:
abarral
Message:

Encapsulate files in modules

File:
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/lmdz_leapfrog.f90

    r5185 r5186  
    1 ! $Id$
    2 
    3 
    4 
    5 SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, q, time_0)
    6 
    7 
    8   !IM : pour sortir les param. du modele dans un fis. netcdf 110106
    9   USE IOIPSL
    10   USE lmdz_infotrac, ONLY: nqtot, isoCheck
    11   USE guide_mod, ONLY: guide_main
    12   USE lmdz_write_field, ONLY: writefield
    13   USE control_mod, ONLY: nday, day_step, planet_type, offline, &
    14           iconser, iphysiq, iperiod, dissip_period, &
    15           iecri, ip_ebil_dyn, ok_dynzon, ok_dyn_ins, &
    16           periodav, ok_dyn_ave, output_grads_dyn
    17   USE exner_hyb_m, ONLY: exner_hyb
    18   USE exner_milieu_m, ONLY: exner_milieu
    19   USE comvert_mod, ONLY: ap, bp, pressure_exner, presnivs
    20   USE comconst_mod, ONLY: cpp, dtphys, dtvr, pi, ihf
    21   USE logic_mod, ONLY: iflag_phys, ok_guide, forward, leapf, apphys, &
    22           statcl, conser, apdiss, purmats, ok_strato
    23   USE temps_mod, ONLY: jD_ref, jH_ref, itaufin, day_ini, day_ref, &
    24           start_time, dt
    25   USE lmdz_strings, ONLY: msg
    26   USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_PHYS
    27   USE lmdz_description, ONLY: descript
    28   USE lmdz_iniprint, ONLY: lunout, prt_level
    29   USE lmdz_ssum_scopy, ONLY: scopy, ssum
    30   USE lmdz_academic, ONLY: tetarappel, knewt_t, kfrict, knewt_g, clat4
    31   USE lmdz_comdissnew, ONLY: lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
    32           tetagrot, tetatemp, coefdis, vert_prof_dissip
    33   USE lmdz_comgeom
    34 
    35   USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
    36   USE lmdz_paramet
    37   IMPLICIT NONE
    38 
    39   ! ......   Version  du 10/01/98    ..........
    40 
    41   !        avec  coordonnees  verticales hybrides
    42   !   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
    43 
    44   !=======================================================================
    45 
    46   !   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
    47   !   -------
    48 
    49   !   Objet:
    50   !   ------
    51 
    52   !   GCM LMD nouvelle grille
    53 
    54   !=======================================================================
    55 
    56   !  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
    57   !  et possibilite d'appeler une fonction f(y)  a derivee tangente
    58   !  hyperbolique a la  place de la fonction a derivee sinusoidale.
    59 
    60   !  ... Possibilite de choisir le shema pour l'advection de
    61   !    q  , en modifiant iadv dans traceur.def  (10/02) .
    62 
    63   !  Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
    64   !  Pour Van-Leer iadv=10
    65 
    66   !-----------------------------------------------------------------------
    67   !   Declarations:
    68   !   -------------
    69 
    70 
    71 
    72 
    73   REAL, INTENT(IN) :: time_0 ! not used
    74 
    75   !   dynamical variables:
    76   REAL, INTENT(INOUT) :: ucov(ip1jmp1, llm)    ! zonal covariant wind
    77   REAL, INTENT(INOUT) :: vcov(ip1jm, llm)      ! meridional covariant wind
    78   REAL, INTENT(INOUT) :: teta(ip1jmp1, llm)    ! potential temperature
    79   REAL, INTENT(INOUT) :: ps(ip1jmp1)          ! surface pressure (Pa)
    80   REAL, INTENT(INOUT) :: masse(ip1jmp1, llm)   ! air mass
    81   REAL, INTENT(INOUT) :: phis(ip1jmp1)        ! geopotentiat at the surface
    82   REAL, INTENT(INOUT) :: q(ip1jmp1, llm, nqtot) ! advected tracers
    83 
    84   REAL :: p (ip1jmp1, llmp1)               ! interlayer pressure
    85   REAL :: pks(ip1jmp1)                      ! exner at the surface
    86   REAL :: pk(ip1jmp1, llm)                   ! exner at mid-layer
    87   REAL :: pkf(ip1jmp1, llm)                  ! filtered exner at mid-layer
    88   REAL :: phi(ip1jmp1, llm)                  ! geopotential
    89   REAL :: w(ip1jmp1, llm)                    ! vertical velocity
    90 
    91   REAL :: zqmin, zqmax
    92 
    93   ! variables dynamiques intermediaire pour le transport
    94   REAL :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm) !flux de masse
    95 
    96   !   variables dynamiques au pas -1
    97   REAL :: vcovm1(ip1jm, llm), ucovm1(ip1jmp1, llm)
    98   REAL :: tetam1(ip1jmp1, llm), psm1(ip1jmp1)
    99   REAL :: massem1(ip1jmp1, llm)
    100 
    101   !   tendances dynamiques
    102   REAL :: dv(ip1jm, llm), du(ip1jmp1, llm)
    103   REAL :: dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqtot), dp(ip1jmp1)
    104 
    105   !   tendances de la dissipation
    106   REAL :: dvdis(ip1jm, llm), dudis(ip1jmp1, llm)
    107   REAL :: dtetadis(ip1jmp1, llm)
    108 
    109   !   tendances physiques
    110   REAL :: dvfi(ip1jm, llm), dufi(ip1jmp1, llm)
    111   REAL :: dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqtot), dpfi(ip1jmp1)
    112 
    113   !   variables pour le fichier histoire
    114   REAL :: dtav      ! intervalle de temps elementaire
    115 
    116   REAL :: tppn(iim), tpps(iim), tpn, tps
    117 
    118   INTEGER :: itau, itaufinp1, iav
    119   ! INTEGER  iday ! jour julien
    120   REAL :: time
    121 
    122   ! REAL finvmaold(ip1jmp1,llm)
    123 
    124   !ym      LOGICAL  lafin
    125   LOGICAL :: lafin = .FALSE.
    126   INTEGER :: ij, iq, l
    127   INTEGER :: ik
    128 
    129   REAL :: time_step, t_wrt, t_ops
    130 
    131   ! REAL rdayvrai,rdaym_ini
    132   ! jD_cur: jour julien courant
    133   ! jH_cur: heure julienne courante
    134   REAL :: jD_cur, jH_cur
    135   INTEGER :: an, mois, jour
    136   REAL :: secondes
    137 
    138   LOGICAL :: first, callinigrads
    139   !IM : pour sortir les param. du modele dans un fis. netcdf 110106
    140   save first
    141   data first/.TRUE./
    142   REAL :: dt_cum
    143   CHARACTER(LEN = 10) :: infile
    144   INTEGER :: zan, tau0, thoriid
    145   INTEGER :: nid_ctesGCM
    146   save nid_ctesGCM
    147   REAL :: degres
    148   REAL :: rlong(iip1), rlatg(jjp1)
    149   REAL :: zx_tmp_2d(iip1, jjp1)
    150   INTEGER :: ndex2d(iip1 * jjp1)
    151   LOGICAL :: ok_sync
    152   parameter (ok_sync = .TRUE.)
    153   LOGICAL :: physic
    154 
    155   data callinigrads/.TRUE./
    156   CHARACTER(LEN = 10) :: string10
    157 
    158   REAL :: flxw(ip1jmp1, llm)  ! flux de masse verticale
    159 
    160   !+jld variables test conservation energie
    161   REAL :: ecin(ip1jmp1, llm), ecin0(ip1jmp1, llm)
    162   ! Tendance de la temp. potentiel d (theta)/ d t due a la
    163   ! tansformation d'energie cinetique en energie thermique
    164   ! cree par la dissipation
    165   REAL :: dtetaecdt(ip1jmp1, llm)
    166   REAL :: vcont(ip1jm, llm), ucont(ip1jmp1, llm)
    167   REAL :: vnat(ip1jm, llm), unat(ip1jmp1, llm)
    168   REAL :: d_h_vcol, d_qt, d_qw, d_ql, d_ec
    169   CHARACTER(len = 15) :: ztit
    170   !IM   INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
    171   !IM   SAVE      ip_ebil_dyn
    172   !IM   DATA      ip_ebil_dyn/0/
    173   !-jld
    174 
    175   CHARACTER(LEN = 80) :: dynhist_file, dynhistave_file
    176   CHARACTER(LEN = *), parameter :: modname = "leapfrog"
    177   CHARACTER(LEN = 80) :: abort_message
    178 
    179   LOGICAL :: dissip_conservative
    180   save dissip_conservative
    181   data dissip_conservative/.TRUE./
    182 
    183   LOGICAL :: prem
    184   save prem
    185   DATA prem/.TRUE./
    186   INTEGER :: testita
    187   PARAMETER (testita = 9)
    188 
    189   logical, parameter :: flag_verif = .FALSE.
    190 
    191   INTEGER :: itau_w   ! pas de temps ecriture = itap + itau_phy
    192 
    193   IF (nday>=0) THEN
    194     itaufin = nday * day_step
    195   else
    196     itaufin = -nday
    197   ENDIF
    198   itaufinp1 = itaufin + 1
    199   itau = 0
    200   physic = .TRUE.
    201   IF (iflag_phys==0.OR.iflag_phys==2) physic = .FALSE.
    202 
    203   ! iday = day_ini+itau/day_step
    204   ! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
    205   !    IF(time.GT.1.) THEN
    206   !     time = time-1.
    207   !     iday = iday+1
    208   !    ENDIF
    209 
    210 
    211   !-----------------------------------------------------------------------
    212   !   On initialise la pression et la fonction d'Exner :
    213   !   --------------------------------------------------
    214 
    215   dq(:, :, :) = 0.
    216   CALL pression (ip1jmp1, ap, bp, ps, p)
    217   IF (pressure_exner) THEN
    218     CALL exner_hyb(ip1jmp1, ps, p, pks, pk, pkf)
    219   else
    220     CALL exner_milieu(ip1jmp1, ps, p, pks, pk, pkf)
    221   ENDIF
    222 
    223   !-----------------------------------------------------------------------
    224   !   Debut de l'integration temporelle:
    225   !   ----------------------------------
    226 
    227   1   CONTINUE ! Matsuno Forward step begins here
    228 
    229   !   date: (NB: date remains unchanged for Backward step)
    230   !   -----
    231 
    232   jD_cur = jD_ref + day_ini - day_ref + &
    233           (itau + 1) / day_step
    234   jH_cur = jH_ref + start_time + &
    235           mod(itau + 1, day_step) / float(day_step)
    236   jD_cur = jD_cur + int(jH_cur)
    237   jH_cur = jH_cur - int(jH_cur)
    238 
    239   CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 321')
    240 
    241   IF (ok_guide) THEN
    242     CALL guide_main(itau, ucov, vcov, teta, q, masse, ps)
    243   ENDIF
    244 
    245 
    246 
    247   ! IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
    248   !   CALL  test_period ( ucov,vcov,teta,q,p,phis )
    249   !   PRINT *,' ----   Test_period apres continue   OK ! -----', itau
    250   ! ENDIF
    251   !
    252 
    253   ! Save fields obtained at previous time step as '...m1'
    254   CALL SCOPY(ijmllm, vcov, 1, vcovm1, 1)
    255   CALL SCOPY(ijp1llm, ucov, 1, ucovm1, 1)
    256   CALL SCOPY(ijp1llm, teta, 1, tetam1, 1)
    257   CALL SCOPY(ijp1llm, masse, 1, massem1, 1)
    258   CALL SCOPY(ip1jmp1, ps, 1, psm1, 1)
    259 
    260   forward = .TRUE.
    261   leapf = .FALSE.
    262   dt = dtvr
    263 
    264   !   ...    P.Le Van .26/04/94  ....
    265   ! Ehouarn: finvmaold is actually not used
    266   ! CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
    267   ! CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
    268 
    269   CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 400')
    270 
    271   2   CONTINUE ! Matsuno backward or leapfrog step begins here
    272 
    273   !-----------------------------------------------------------------------
    274 
    275   !   date: (NB: only leapfrog step requires recomputing date)
    276   !   -----
    277 
    278   IF (leapf) THEN
     1MODULE lmdz_leapfrog
     2
     3  IMPLICIT NONE; PRIVATE
     4  PUBLIC leapfrog
     5
     6CONTAINS
     7
     8  SUBROUTINE leapfrog(ucov, vcov, teta, ps, masse, phis, q, time_0)
     9
     10    !IM : pour sortir les param. du modele dans un fis. netcdf 110106
     11    USE IOIPSL
     12    USE lmdz_infotrac, ONLY: nqtot, isoCheck
     13    USE guide_mod, ONLY: guide_main
     14    USE lmdz_write_field, ONLY: writefield
     15    USE control_mod, ONLY: nday, day_step, planet_type, offline, &
     16            iconser, iphysiq, iperiod, dissip_period, &
     17            iecri, ip_ebil_dyn, ok_dynzon, ok_dyn_ins, &
     18            periodav, ok_dyn_ave, output_grads_dyn
     19    USE exner_hyb_m, ONLY: exner_hyb
     20    USE exner_milieu_m, ONLY: exner_milieu
     21    USE comvert_mod, ONLY: ap, bp, pressure_exner, presnivs
     22    USE comconst_mod, ONLY: cpp, dtphys, dtvr, pi, ihf
     23    USE logic_mod, ONLY: iflag_phys, ok_guide, forward, leapf, apphys, &
     24            statcl, conser, apdiss, purmats, ok_strato
     25    USE temps_mod, ONLY: jD_ref, jH_ref, itaufin, day_ini, day_ref, &
     26            start_time, dt
     27    USE lmdz_strings, ONLY: msg
     28    USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_PHYS
     29    USE lmdz_description, ONLY: descript
     30    USE lmdz_iniprint, ONLY: lunout, prt_level
     31    USE lmdz_ssum_scopy, ONLY: scopy, ssum
     32    USE lmdz_academic, ONLY: tetarappel, knewt_t, kfrict, knewt_g, clat4
     33    USE lmdz_comdissnew, ONLY: lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
     34            tetagrot, tetatemp, coefdis, vert_prof_dissip
     35    USE lmdz_comgeom
     36
     37    USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
     38    USE lmdz_paramet
     39    USE lmdz_addfi, ONLY: addfi
     40    USE lmdz_caldyn, ONLY: caldyn
     41    USE lmdz_caladvtrac, ONLY: caladvtrac
     42    USE lmdz_bilan_dyn, ONLY: bilan_dyn
     43    USE lmdz_check_isotopes, ONLY: check_isotopes_seq
     44    USE lmdz_writedynav, ONLY: writedynav
     45    USE lmdz_writehist, ONLY: writehist
     46    USE lmdz_dissip, ONLY: dissip
     47    USE lmdz_dynredem, ONLY: dynredem1
     48
     49    IMPLICIT NONE
     50
     51    ! ......   Version  du 10/01/98    ..........
     52
     53    !        avec  coordonnees  verticales hybrides
     54    !   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
     55
     56    !=======================================================================
     57
     58    !   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
     59    !   -------
     60
     61    !   Objet:
     62    !   ------
     63
     64    !   GCM LMD nouvelle grille
     65
     66    !=======================================================================
     67
     68    !  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
     69    !  et possibilite d'appeler une fonction f(y)  a derivee tangente
     70    !  hyperbolique a la  place de la fonction a derivee sinusoidale.
     71
     72    !  ... Possibilite de choisir le shema pour l'advection de
     73    !    q  , en modifiant iadv dans traceur.def  (10/02) .
     74
     75    !  Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
     76    !  Pour Van-Leer iadv=10
     77
     78    !-----------------------------------------------------------------------
     79    !   Declarations:
     80    !   -------------
     81
     82    REAL, INTENT(IN) :: time_0 ! not used
     83
     84    !   dynamical variables:
     85    REAL, INTENT(INOUT) :: ucov(ip1jmp1, llm)    ! zonal covariant wind
     86    REAL, INTENT(INOUT) :: vcov(ip1jm, llm)      ! meridional covariant wind
     87    REAL, INTENT(INOUT) :: teta(ip1jmp1, llm)    ! potential temperature
     88    REAL, INTENT(INOUT) :: ps(ip1jmp1)          ! surface pressure (Pa)
     89    REAL, INTENT(INOUT) :: masse(ip1jmp1, llm)   ! air mass
     90    REAL, INTENT(INOUT) :: phis(ip1jmp1)        ! geopotentiat at the surface
     91    REAL, INTENT(INOUT) :: q(ip1jmp1, llm, nqtot) ! advected tracers
     92
     93    REAL :: p (ip1jmp1, llmp1)               ! interlayer pressure
     94    REAL :: pks(ip1jmp1)                      ! exner at the surface
     95    REAL :: pk(ip1jmp1, llm)                   ! exner at mid-layer
     96    REAL :: pkf(ip1jmp1, llm)                  ! filtered exner at mid-layer
     97    REAL :: phi(ip1jmp1, llm)                  ! geopotential
     98    REAL :: w(ip1jmp1, llm)                    ! vertical velocity
     99
     100    REAL :: zqmin, zqmax
     101
     102    ! variables dynamiques intermediaire pour le transport
     103    REAL :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm) !flux de masse
     104
     105    !   variables dynamiques au pas -1
     106    REAL :: vcovm1(ip1jm, llm), ucovm1(ip1jmp1, llm)
     107    REAL :: tetam1(ip1jmp1, llm), psm1(ip1jmp1)
     108    REAL :: massem1(ip1jmp1, llm)
     109
     110    !   tendances dynamiques
     111    REAL :: dv(ip1jm, llm), du(ip1jmp1, llm)
     112    REAL :: dteta(ip1jmp1, llm), dq(ip1jmp1, llm, nqtot), dp(ip1jmp1)
     113
     114    !   tendances de la dissipation
     115    REAL :: dvdis(ip1jm, llm), dudis(ip1jmp1, llm)
     116    REAL :: dtetadis(ip1jmp1, llm)
     117
     118    !   tendances physiques
     119    REAL :: dvfi(ip1jm, llm), dufi(ip1jmp1, llm)
     120    REAL :: dtetafi(ip1jmp1, llm), dqfi(ip1jmp1, llm, nqtot), dpfi(ip1jmp1)
     121
     122    !   variables pour le fichier histoire
     123    REAL :: dtav      ! intervalle de temps elementaire
     124
     125    REAL :: tppn(iim), tpps(iim), tpn, tps
     126
     127    INTEGER :: itau, itaufinp1, iav
     128    ! INTEGER  iday ! jour julien
     129    REAL :: time
     130
     131    ! REAL finvmaold(ip1jmp1,llm)
     132
     133    !ym      LOGICAL  lafin
     134    LOGICAL :: lafin = .FALSE.
     135    INTEGER :: ij, iq, l
     136    INTEGER :: ik
     137
     138    REAL :: time_step, t_wrt, t_ops
     139
     140    ! REAL rdayvrai,rdaym_ini
     141    ! jD_cur: jour julien courant
     142    ! jH_cur: heure julienne courante
     143    REAL :: jD_cur, jH_cur
     144    INTEGER :: an, mois, jour
     145    REAL :: secondes
     146
     147    LOGICAL :: first, callinigrads
     148    !IM : pour sortir les param. du modele dans un fis. netcdf 110106
     149    save first
     150    data first/.TRUE./
     151    REAL :: dt_cum
     152    CHARACTER(LEN = 10) :: infile
     153    INTEGER :: zan, tau0, thoriid
     154    INTEGER :: nid_ctesGCM
     155    save nid_ctesGCM
     156    REAL :: degres
     157    REAL :: rlong(iip1), rlatg(jjp1)
     158    REAL :: zx_tmp_2d(iip1, jjp1)
     159    INTEGER :: ndex2d(iip1 * jjp1)
     160    LOGICAL :: ok_sync
     161    parameter (ok_sync = .TRUE.)
     162    LOGICAL :: physic
     163
     164    data callinigrads/.TRUE./
     165    CHARACTER(LEN = 10) :: string10
     166
     167    REAL :: flxw(ip1jmp1, llm)  ! flux de masse verticale
     168
     169    !+jld variables test conservation energie
     170    REAL :: ecin(ip1jmp1, llm), ecin0(ip1jmp1, llm)
     171    ! Tendance de la temp. potentiel d (theta)/ d t due a la
     172    ! tansformation d'energie cinetique en energie thermique
     173    ! cree par la dissipation
     174    REAL :: dtetaecdt(ip1jmp1, llm)
     175    REAL :: vcont(ip1jm, llm), ucont(ip1jmp1, llm)
     176    REAL :: vnat(ip1jm, llm), unat(ip1jmp1, llm)
     177    REAL :: d_h_vcol, d_qt, d_qw, d_ql, d_ec
     178    CHARACTER(len = 15) :: ztit
     179    !IM   INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
     180    !IM   SAVE      ip_ebil_dyn
     181    !IM   DATA      ip_ebil_dyn/0/
     182    !-jld
     183
     184    CHARACTER(LEN = 80) :: dynhist_file, dynhistave_file
     185    CHARACTER(LEN = *), parameter :: modname = "leapfrog"
     186    CHARACTER(LEN = 80) :: abort_message
     187
     188    LOGICAL :: dissip_conservative
     189    save dissip_conservative
     190    data dissip_conservative/.TRUE./
     191
     192    LOGICAL :: prem
     193    save prem
     194    DATA prem/.TRUE./
     195    INTEGER :: testita
     196    PARAMETER (testita = 9)
     197
     198    logical, parameter :: flag_verif = .FALSE.
     199
     200    INTEGER :: itau_w   ! pas de temps ecriture = itap + itau_phy
     201
     202    IF (nday>=0) THEN
     203      itaufin = nday * day_step
     204    else
     205      itaufin = -nday
     206    ENDIF
     207    itaufinp1 = itaufin + 1
     208    itau = 0
     209    physic = .TRUE.
     210    IF (iflag_phys==0.OR.iflag_phys==2) physic = .FALSE.
     211
     212    ! iday = day_ini+itau/day_step
     213    ! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
     214    !    IF(time.GT.1.) THEN
     215    !     time = time-1.
     216    !     iday = iday+1
     217    !    ENDIF
     218
     219
     220    !-----------------------------------------------------------------------
     221    !   On initialise la pression et la fonction d'Exner :
     222    !   --------------------------------------------------
     223
     224    dq(:, :, :) = 0.
     225    CALL pression (ip1jmp1, ap, bp, ps, p)
     226    IF (pressure_exner) THEN
     227      CALL exner_hyb(ip1jmp1, ps, p, pks, pk, pkf)
     228    else
     229      CALL exner_milieu(ip1jmp1, ps, p, pks, pk, pkf)
     230    ENDIF
     231
     232    !-----------------------------------------------------------------------
     233    !   Debut de l'integration temporelle:
     234    !   ----------------------------------
     235
     236    1   CONTINUE ! Matsuno Forward step begins here
     237
     238    !   date: (NB: date remains unchanged for Backward step)
     239    !   -----
     240
    279241    jD_cur = jD_ref + day_ini - day_ref + &
    280242            (itau + 1) / day_step
    … …  
    283245    jD_cur = jD_cur + int(jH_cur)
    284246    jH_cur = jH_cur - int(jH_cur)
    285   ENDIF
    286 
    287 
    288   !   gestion des appels de la physique et des dissipations:
    289   !   ------------------------------------------------------
    290 
    291   !   ...    P.Le Van  ( 6/02/95 )  ....
    292 
    293   apphys = .FALSE.
    294   statcl = .FALSE.
    295   conser = .FALSE.
    296   apdiss = .FALSE.
    297 
    298   IF(purmats) THEN
    299     ! Purely Matsuno time stepping
    300     IF(MOD(itau, iconser) ==0.AND.  forward) conser = .TRUE.
    301     IF(MOD(itau, dissip_period)==0.AND..NOT.forward) &
    302             apdiss = .TRUE.
    303     IF(MOD(itau, iphysiq)==0.AND..NOT.forward &
    304             .AND. physic) apphys = .TRUE.
    305   ELSE
    306     ! Leapfrog/Matsuno time stepping
    307     IF(MOD(itau, iconser) == 0) conser = .TRUE.
    308     IF(MOD(itau + 1, dissip_period)==0 .AND. .NOT. forward) &
    309             apdiss = .TRUE.
    310     IF(MOD(itau + 1, iphysiq)==0.AND.physic) apphys = .TRUE.
    311   END IF
    312 
    313   ! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
    314   ! supress dissipation step
    315   IF (llm==1) THEN
     247
     248    CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 321')
     249
     250    IF (ok_guide) THEN
     251      CALL guide_main(itau, ucov, vcov, teta, q, masse, ps)
     252    ENDIF
     253
     254
     255
     256    ! IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
     257    !   CALL  test_period ( ucov,vcov,teta,q,p,phis )
     258    !   PRINT *,' ----   Test_period apres continue   OK ! -----', itau
     259    ! ENDIF
     260    !
     261
     262    ! Save fields obtained at previous time step as '...m1'
     263    CALL SCOPY(ijmllm, vcov, 1, vcovm1, 1)
     264    CALL SCOPY(ijp1llm, ucov, 1, ucovm1, 1)
     265    CALL SCOPY(ijp1llm, teta, 1, tetam1, 1)
     266    CALL SCOPY(ijp1llm, masse, 1, massem1, 1)
     267    CALL SCOPY(ip1jmp1, ps, 1, psm1, 1)
     268
     269    forward = .TRUE.
     270    leapf = .FALSE.
     271    dt = dtvr
     272
     273    !   ...    P.Le Van .26/04/94  ....
     274    ! Ehouarn: finvmaold is actually not used
     275    ! CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
     276    ! CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
     277
     278    CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 400')
     279
     280    2   CONTINUE ! Matsuno backward or leapfrog step begins here
     281
     282    !-----------------------------------------------------------------------
     283
     284    !   date: (NB: only leapfrog step requires recomputing date)
     285    !   -----
     286
     287    IF (leapf) THEN
     288      jD_cur = jD_ref + day_ini - day_ref + &
     289              (itau + 1) / day_step
     290      jH_cur = jH_ref + start_time + &
     291              mod(itau + 1, day_step) / float(day_step)
     292      jD_cur = jD_cur + int(jH_cur)
     293      jH_cur = jH_cur - int(jH_cur)
     294    ENDIF
     295
     296
     297    !   gestion des appels de la physique et des dissipations:
     298    !   ------------------------------------------------------
     299
     300    !   ...    P.Le Van  ( 6/02/95 )  ....
     301
     302    apphys = .FALSE.
     303    statcl = .FALSE.
     304    conser = .FALSE.
    316305    apdiss = .FALSE.
    317   ENDIF
    318 
    319   CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 589')
    320 
    321   !-----------------------------------------------------------------------
    322   !   calcul des tendances dynamiques:
    323   !   --------------------------------
    324 
    325   ! compute geopotential phi()
    326   CALL geopot  (ip1jmp1, teta, pk, pks, phis, phi)
    327 
    328   time = jD_cur + jH_cur
    329   CALL caldyn &
    330           (itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, &
    331           phi, conser, du, dv, dteta, dp, w, pbaru, pbarv, time)
    332 
    333 
    334   !-----------------------------------------------------------------------
    335   !   calcul des tendances advection des traceurs (dont l'humidite)
    336   !   -------------------------------------------------------------
    337 
    338   CALL check_isotopes_seq(q, ip1jmp1, &
    339           'leapfrog 686: avant caladvtrac')
    340 
    341   IF(forward .OR.  leapf)  THEN
    342     ! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
    343     CALL caladvtrac(q, pbaru, pbarv, &
    344             p, masse, dq, teta, &
    345             flxw, pk)
    346     !WRITE(*,*) 'caladvtrac 346'
    347 
    348     IF (offline) THEN
    349       !maf stokage du flux de masse pour traceurs OFF-LINE
    350 
    351       CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
    352               dtvr, itau)
    353 
    354     ENDIF ! of IF (offline)
    355 
    356   ENDIF ! of IF( forward .OR.  leapf )
    357 
    358 
    359   !-----------------------------------------------------------------------
    360   !   integrations dynamique et traceurs:
    361   !   ----------------------------------
    362 
    363   CALL msg('720', modname, isoCheck)
    364   CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 756')
    365 
    366   CALL integrd (nqtot, vcovm1, ucovm1, tetam1, psm1, massem1, &
    367           dv, du, dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis)
    368   ! $              finvmaold                                    )
    369 
    370   CALL msg('724', modname, isoCheck)
    371   CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 762')
    372 
    373   ! .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
    374 
    375   !-----------------------------------------------------------------------
    376   !   calcul des tendances physiques:
    377   !   -------------------------------
    378   !    ########   P.Le Van ( Modif le  6/02/95 )   ###########
    379 
    380   IF(purmats)  THEN
    381     IF(itau==itaufin.AND..NOT.forward) lafin = .TRUE.
    382   ELSE
    383     IF(itau + 1 == itaufin)              lafin = .TRUE.
    384   ENDIF
    385 
    386 
    387   IF(apphys)  THEN
    388 
    389     ! .......   Ajout   P.Le Van ( 17/04/96 )   ...........
    390     !
     306
     307    IF(purmats) THEN
     308      ! Purely Matsuno time stepping
     309      IF(MOD(itau, iconser) ==0.AND.  forward) conser = .TRUE.
     310      IF(MOD(itau, dissip_period)==0.AND..NOT.forward) &
     311              apdiss = .TRUE.
     312      IF(MOD(itau, iphysiq)==0.AND..NOT.forward &
     313              .AND. physic) apphys = .TRUE.
     314    ELSE
     315      ! Leapfrog/Matsuno time stepping
     316      IF(MOD(itau, iconser) == 0) conser = .TRUE.
     317      IF(MOD(itau + 1, dissip_period)==0 .AND. .NOT. forward) &
     318              apdiss = .TRUE.
     319      IF(MOD(itau + 1, iphysiq)==0.AND.physic) apphys = .TRUE.
     320    END IF
     321
     322    ! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
     323    ! supress dissipation step
     324    IF (llm==1) THEN
     325      apdiss = .FALSE.
     326    ENDIF
     327
     328    CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 589')
     329
     330    !-----------------------------------------------------------------------
     331    !   calcul des tendances dynamiques:
     332    !   --------------------------------
     333
     334    ! compute geopotential phi()
     335    CALL geopot  (ip1jmp1, teta, pk, pks, phis, phi)
     336
     337    time = jD_cur + jH_cur
     338    CALL caldyn &
     339            (itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, &
     340            phi, conser, du, dv, dteta, dp, w, pbaru, pbarv, time)
     341
     342
     343    !-----------------------------------------------------------------------
     344    !   calcul des tendances advection des traceurs (dont l'humidite)
     345    !   -------------------------------------------------------------
     346
     347    CALL check_isotopes_seq(q, ip1jmp1, &
     348            'leapfrog 686: avant caladvtrac')
     349
     350    IF(forward .OR.  leapf)  THEN
     351      ! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
     352      CALL caladvtrac(q, pbaru, pbarv, &
     353              p, masse, dq, teta, &
     354              flxw, pk)
     355      !WRITE(*,*) 'caladvtrac 346'
     356
     357      IF (offline) THEN
     358        !maf stokage du flux de masse pour traceurs OFF-LINE
     359
     360        CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
     361                dtvr, itau)
     362
     363      ENDIF ! of IF (offline)
     364
     365    ENDIF ! of IF( forward .OR.  leapf )
     366
     367
     368    !-----------------------------------------------------------------------
     369    !   integrations dynamique et traceurs:
     370    !   ----------------------------------
     371
     372    CALL msg('720', modname, isoCheck)
     373    CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 756')
     374
     375    CALL integrd (nqtot, vcovm1, ucovm1, tetam1, psm1, massem1, &
     376            dv, du, dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis)
     377    ! $              finvmaold                                    )
     378
     379    CALL msg('724', modname, isoCheck)
     380    CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 762')
     381
     382    ! .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
     383
     384    !-----------------------------------------------------------------------
     385    !   calcul des tendances physiques:
     386    !   -------------------------------
     387    !    ########   P.Le Van ( Modif le  6/02/95 )   ###########
     388
     389    IF(purmats)  THEN
     390      IF(itau==itaufin.AND..NOT.forward) lafin = .TRUE.
     391    ELSE
     392      IF(itau + 1 == itaufin)              lafin = .TRUE.
     393    ENDIF
     394
     395    IF(apphys)  THEN
     396
     397      ! .......   Ajout   P.Le Van ( 17/04/96 )   ...........
     398      !
     399
     400      CALL pression (ip1jmp1, ap, bp, ps, p)
     401      IF (pressure_exner) THEN
     402        CALL exner_hyb(ip1jmp1, ps, p, pks, pk, pkf)
     403      else
     404        CALL exner_milieu(ip1jmp1, ps, p, pks, pk, pkf)
     405      endif
     406
     407      ! Appel a geopot ajoute le 2014/05/08 pour garantir la convergence numerique
     408      ! avec dyn3dmem
     409      CALL geopot  (ip1jmp1, teta, pk, pks, phis, phi)
     410
     411      ! rdaym_ini  = itau * dtvr / daysec
     412      ! rdayvrai   = rdaym_ini  + day_ini
     413      ! jD_cur = jD_ref + day_ini - day_ref
     414      ! $        + int (itau * dtvr / daysec)
     415      !       jH_cur = jH_ref +                                            &
     416      ! &              (itau * dtvr / daysec - int(itau * dtvr / daysec))
     417      jD_cur = jD_ref + day_ini - day_ref + &
     418              (itau + 1) / day_step
     419
     420      IF (planet_type =="generic") THEN
     421        ! AS: we make jD_cur to be pday
     422        jD_cur = int(day_ini + itau / day_step)
     423      ENDIF
     424
     425      jH_cur = jH_ref + start_time + &
     426              mod(itau + 1, day_step) / float(day_step)
     427      jD_cur = jD_cur + int(jH_cur)
     428      jH_cur = jH_cur - int(jH_cur)
     429      ! WRITE(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
     430      ! CALL ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
     431      ! WRITE(lunout,*)'current date = ',an, mois, jour, secondes
     432
     433      ! rajout debug
     434      ! lafin = .TRUE.
     435
     436
     437      !   Inbterface avec les routines de phylmd (phymars ... )
     438      !   -----------------------------------------------------
     439
     440      !+jld
     441
     442      !  Diagnostique de conservation de l'energie : initialisation
     443      IF (ip_ebil_dyn>=1) THEN
     444        ztit = 'bil dyn'
     445        ! Ehouarn: be careful, diagedyn is Earth-specific!
     446        IF (planet_type=="earth") THEN
     447          CALL diagedyn(ztit, 2, 1, 1, dtphys &
     448                  , ucov, vcov, ps, p, pk, teta, q(:, :, 1), q(:, :, 2))
     449        ENDIF
     450      ENDIF ! of IF (ip_ebil_dyn.ge.1 )
     451      IF (CPPKEY_PHYS) THEN
     452        CALL calfis(lafin, jD_cur, jH_cur, &
     453                ucov, vcov, teta, q, masse, ps, p, pk, phis, phi, &
     454                du, dv, dteta, dq, &
     455                flxw, dufi, dvfi, dtetafi, dqfi, dpfi)
     456      END IF
     457      ! ajout des tendances physiques:
     458      ! ------------------------------
     459      CALL addfi(dtphys, leapf, forward, &
     460              ucov, vcov, teta, q, ps, &
     461              dufi, dvfi, dtetafi, dqfi, dpfi)
     462      ! since addfi updates ps(), also update p(), masse() and pk()
     463      CALL pression (ip1jmp1, ap, bp, ps, p)
     464      CALL massdair(p, masse)
     465      IF (pressure_exner) THEN
     466        CALL exner_hyb(ip1jmp1, ps, p, pks, pk, pkf)
     467      else
     468        CALL exner_milieu(ip1jmp1, ps, p, pks, pk, pkf)
     469      endif
     470
     471      IF (ok_strato) THEN
     472        CALL top_bound(vcov, ucov, teta, masse, dtphys)
     473      ENDIF
     474
     475
     476      !  Diagnostique de conservation de l'energie : difference
     477      IF (ip_ebil_dyn>=1) THEN
     478        ztit = 'bil phys'
     479        IF (planet_type=="earth") THEN
     480          CALL diagedyn(ztit, 2, 1, 1, dtphys &
     481                  , ucov, vcov, ps, p, pk, teta, q(:, :, 1), q(:, :, 2))
     482        ENDIF
     483      ENDIF ! of IF (ip_ebil_dyn.ge.1 )
     484
     485    ENDIF ! of IF( apphys )
     486
     487    IF(iflag_phys==2) THEN ! "Newtonian" case
     488      !   Academic case : Simple friction and Newtonan relaxation
     489      !   -------------------------------------------------------
     490      DO l = 1, llm
     491        DO ij = 1, ip1jmp1
     492          teta(ij, l) = teta(ij, l) - dtvr * &
     493                  (teta(ij, l) - tetarappel(ij, l)) * (knewt_g + knewt_t(l) * clat4(ij))
     494        ENDDO
     495      ENDDO ! of DO l=1,llm
     496
     497      IF (planet_type=="giant") THEN
     498        ! add an intrinsic heat flux at the base of the atmosphere
     499        teta(:, 1) = teta(:, 1) + dtvr * aire(:) * ihf / cpp / masse(:, 1)
     500      endif
     501
     502      CALL friction(ucov, vcov, dtvr)
     503
     504      ! Sponge layer (if any)
     505      IF (ok_strato) THEN
     506        ! dufi(:,:)=0.
     507        ! dvfi(:,:)=0.
     508        ! dtetafi(:,:)=0.
     509        ! dqfi(:,:,:)=0.
     510        !          dpfi(:)=0.
     511        ! CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
     512        CALL top_bound(vcov, ucov, teta, masse, dtvr)
     513        ! CALL addfi( dtvr, leapf, forward   ,
     514        ! $                  ucov, vcov, teta , q   ,ps ,
     515        ! $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
     516      ENDIF ! of IF (ok_strato)
     517    ENDIF ! of IF (iflag_phys.EQ.2)
     518
     519
     520    !-jld
    391521
    392522    CALL pression (ip1jmp1, ap, bp, ps, p)
    … …  
    395525    else
    396526      CALL exner_milieu(ip1jmp1, ps, p, pks, pk, pkf)
    397     endif
    398 
    399     ! Appel a geopot ajoute le 2014/05/08 pour garantir la convergence numerique
    400     ! avec dyn3dmem
    401     CALL geopot  (ip1jmp1, teta, pk, pks, phis, phi)
    402 
    403     ! rdaym_ini  = itau * dtvr / daysec
    404     ! rdayvrai   = rdaym_ini  + day_ini
    405     ! jD_cur = jD_ref + day_ini - day_ref
    406     ! $        + int (itau * dtvr / daysec)
    407     !       jH_cur = jH_ref +                                            &
    408     ! &              (itau * dtvr / daysec - int(itau * dtvr / daysec))
    409     jD_cur = jD_ref + day_ini - day_ref + &
    410             (itau + 1) / day_step
    411 
    412     IF (planet_type =="generic") THEN
    413       ! AS: we make jD_cur to be pday
    414       jD_cur = int(day_ini + itau / day_step)
    415527    ENDIF
    416 
    417     jH_cur = jH_ref + start_time + &
    418             mod(itau + 1, day_step) / float(day_step)
    419     jD_cur = jD_cur + int(jH_cur)
    420     jH_cur = jH_cur - int(jH_cur)
    421     ! WRITE(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
    422     ! CALL ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
    423     ! WRITE(lunout,*)'current date = ',an, mois, jour, secondes
    424 
    425     ! rajout debug
    426     ! lafin = .TRUE.
    427 
    428 
    429     !   Inbterface avec les routines de phylmd (phymars ... )
    430     !   -----------------------------------------------------
    431 
    432     !+jld
    433 
    434     !  Diagnostique de conservation de l'energie : initialisation
    435     IF (ip_ebil_dyn>=1) THEN
    436       ztit = 'bil dyn'
    437       ! Ehouarn: be careful, diagedyn is Earth-specific!
    438       IF (planet_type=="earth") THEN
    439         CALL diagedyn(ztit, 2, 1, 1, dtphys &
    440                 , ucov, vcov, ps, p, pk, teta, q(:, :, 1), q(:, :, 2))
     528    CALL massdair(p, masse)
     529
     530    CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1196')
     531
     532    !-----------------------------------------------------------------------
     533    !   dissipation horizontale et verticale  des petites echelles:
     534    !   ----------------------------------------------------------
     535
     536    IF(apdiss) THEN
     537
     538
     539      !   calcul de l'energie cinetique avant dissipation
     540      CALL covcont(llm, ucov, vcov, ucont, vcont)
     541      CALL enercin(vcov, ucov, vcont, ucont, ecin0)
     542
     543      !   dissipation
     544      CALL dissip(vcov, ucov, teta, p, dvdis, dudis, dtetadis)
     545      ucov = ucov + dudis
     546      vcov = vcov + dvdis
     547      ! teta=teta+dtetadis
     548
     549
     550      !------------------------------------------------------------------------
     551      IF (dissip_conservative) THEN
     552        ! On rajoute la tendance due a la transform. Ec -> E therm. cree
     553        ! lors de la dissipation
     554        CALL covcont(llm, ucov, vcov, ucont, vcont)
     555        CALL enercin(vcov, ucov, vcont, ucont, ecin)
     556        dtetaecdt = (ecin0 - ecin) / pk
     557        ! teta=teta+dtetaecdt
     558        dtetadis = dtetadis + dtetaecdt
     559      endif
     560      teta = teta + dtetadis
     561      !------------------------------------------------------------------------
     562
     563
     564      !    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
     565      !   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
     566      !
     567
     568      DO l = 1, llm
     569        DO ij = 1, iim
     570          tppn(ij) = aire(ij) * teta(ij, l)
     571          tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)
     572        ENDDO
     573        tpn = SSUM(iim, tppn, 1) / apoln
     574        tps = SSUM(iim, tpps, 1) / apols
     575
     576        DO ij = 1, iip1
     577          teta(ij, l) = tpn
     578          teta(ij + ip1jm, l) = tps
     579        ENDDO
     580      ENDDO
     581
     582      IF (1 == 0) THEN
     583        !!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
     584        !!!                     2) should probably not be here anyway
     585        !!! but are kept for those who would want to revert to previous behaviour
     586        DO ij = 1, iim
     587          tppn(ij) = aire(ij) * ps (ij)
     588          tpps(ij) = aire(ij + ip1jm) * ps (ij + ip1jm)
     589        ENDDO
     590        tpn = SSUM(iim, tppn, 1) / apoln
     591        tps = SSUM(iim, tpps, 1) / apols
     592
     593        DO ij = 1, iip1
     594          ps(ij) = tpn
     595          ps(ij + ip1jm) = tps
     596        ENDDO
     597      endif ! of if (1 == 0)
     598
     599    END IF ! of IF(apdiss)
     600
     601    ! ajout debug
     602    ! IF( lafin ) THEN
     603    !   abort_message = 'Simulation finished'
     604    !   CALL abort_gcm(modname,abort_message,0)
     605    ! ENDIF
     606
     607    !   ********************************************************************
     608    !   ********************************************************************
     609    !   .... fin de l'integration dynamique  et physique pour le pas itau ..
     610    !   ********************************************************************
     611    !   ********************************************************************
     612
     613    !   preparation du pas d'integration suivant  ......
     614
     615    CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1509')
     616
     617    IF (.NOT.purmats) THEN
     618      ! ........................................................
     619      ! ..............  schema matsuno + leapfrog  ..............
     620      ! ........................................................
     621
     622      IF(forward .OR. leapf) THEN
     623        itau = itau + 1
     624        ! iday= day_ini+itau/day_step
     625        ! time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
     626        !   IF(time.GT.1.) THEN
     627        !     time = time-1.
     628        !     iday = iday+1
     629        !   ENDIF
    441630      ENDIF
    442     ENDIF ! of IF (ip_ebil_dyn.ge.1 )
    443     IF (CPPKEY_PHYS) THEN
    444       CALL calfis(lafin, jD_cur, jH_cur, &
    445               ucov, vcov, teta, q, masse, ps, p, pk, phis, phi, &
    446               du, dv, dteta, dq, &
    447               flxw, dufi, dvfi, dtetafi, dqfi, dpfi)
    448     END IF
    449     ! ajout des tendances physiques:
    450     ! ------------------------------
    451     CALL addfi(dtphys, leapf, forward, &
    452             ucov, vcov, teta, q, ps, &
    453             dufi, dvfi, dtetafi, dqfi, dpfi)
    454     ! since addfi updates ps(), also update p(), masse() and pk()
    455     CALL pression (ip1jmp1, ap, bp, ps, p)
    456     CALL massdair(p, masse)
    457     IF (pressure_exner) THEN
    458       CALL exner_hyb(ip1jmp1, ps, p, pks, pk, pkf)
    459     else
    460       CALL exner_milieu(ip1jmp1, ps, p, pks, pk, pkf)
    461     endif
    462 
    463     IF (ok_strato) THEN
    464       CALL top_bound(vcov, ucov, teta, masse, dtphys)
    465     ENDIF
    466 
    467 
    468     !  Diagnostique de conservation de l'energie : difference
    469     IF (ip_ebil_dyn>=1) THEN
    470       ztit = 'bil phys'
    471       IF (planet_type=="earth") THEN
    472         CALL diagedyn(ztit, 2, 1, 1, dtphys &
    473                 , ucov, vcov, ps, p, pk, teta, q(:, :, 1), q(:, :, 2))
    474       ENDIF
    475     ENDIF ! of IF (ip_ebil_dyn.ge.1 )
    476 
    477   ENDIF ! of IF( apphys )
    478 
    479   IF(iflag_phys==2) THEN ! "Newtonian" case
    480     !   Academic case : Simple friction and Newtonan relaxation
    481     !   -------------------------------------------------------
    482     DO l = 1, llm
    483       DO ij = 1, ip1jmp1
    484         teta(ij, l) = teta(ij, l) - dtvr * &
    485                 (teta(ij, l) - tetarappel(ij, l)) * (knewt_g + knewt_t(l) * clat4(ij))
    486       ENDDO
    487     ENDDO ! of DO l=1,llm
    488 
    489     IF (planet_type=="giant") THEN
    490       ! add an intrinsic heat flux at the base of the atmosphere
    491       teta(:, 1) = teta(:, 1) + dtvr * aire(:) * ihf / cpp / masse(:, 1)
    492     endif
    493 
    494     CALL friction(ucov, vcov, dtvr)
    495 
    496     ! Sponge layer (if any)
    497     IF (ok_strato) THEN
    498       ! dufi(:,:)=0.
    499       ! dvfi(:,:)=0.
    500       ! dtetafi(:,:)=0.
    501       ! dqfi(:,:,:)=0.
    502       !          dpfi(:)=0.
    503       ! CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
    504       CALL top_bound(vcov, ucov, teta, masse, dtvr)
    505       ! CALL addfi( dtvr, leapf, forward   ,
    506       ! $                  ucov, vcov, teta , q   ,ps ,
    507       ! $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
    508     ENDIF ! of IF (ok_strato)
    509   ENDIF ! of IF (iflag_phys.EQ.2)
    510 
    511 
    512   !-jld
    513 
    514   CALL pression (ip1jmp1, ap, bp, ps, p)
    515   IF (pressure_exner) THEN
    516     CALL exner_hyb(ip1jmp1, ps, p, pks, pk, pkf)
    517   else
    518     CALL exner_milieu(ip1jmp1, ps, p, pks, pk, pkf)
    519   ENDIF
    520   CALL massdair(p, masse)
    521 
    522   CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1196')
    523 
    524   !-----------------------------------------------------------------------
    525   !   dissipation horizontale et verticale  des petites echelles:
    526   !   ----------------------------------------------------------
    527 
    528   IF(apdiss) THEN
    529 
    530 
    531     !   calcul de l'energie cinetique avant dissipation
    532     CALL covcont(llm, ucov, vcov, ucont, vcont)
    533     CALL enercin(vcov, ucov, vcont, ucont, ecin0)
    534 
    535     !   dissipation
    536     CALL dissip(vcov, ucov, teta, p, dvdis, dudis, dtetadis)
    537     ucov = ucov + dudis
    538     vcov = vcov + dvdis
    539     ! teta=teta+dtetadis
    540 
    541 
    542     !------------------------------------------------------------------------
    543     IF (dissip_conservative) THEN
    544       ! On rajoute la tendance due a la transform. Ec -> E therm. cree
    545       ! lors de la dissipation
    546       CALL covcont(llm, ucov, vcov, ucont, vcont)
    547       CALL enercin(vcov, ucov, vcont, ucont, ecin)
    548       dtetaecdt = (ecin0 - ecin) / pk
    549       ! teta=teta+dtetaecdt
    550       dtetadis = dtetadis + dtetaecdt
    551     endif
    552     teta = teta + dtetadis
    553     !------------------------------------------------------------------------
    554 
    555 
    556     !    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
    557     !   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
    558     !
    559 
    560     DO l = 1, llm
    561       DO ij = 1, iim
    562         tppn(ij) = aire(ij) * teta(ij, l)
    563         tpps(ij) = aire(ij + ip1jm) * teta(ij + ip1jm, l)
    564       ENDDO
    565       tpn = SSUM(iim, tppn, 1) / apoln
    566       tps = SSUM(iim, tpps, 1) / apols
    567 
    568       DO ij = 1, iip1
    569         teta(ij, l) = tpn
    570         teta(ij + ip1jm, l) = tps
    571       ENDDO
    572     ENDDO
    573 
    574     IF (1 == 0) THEN
    575       !!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
    576       !!!                     2) should probably not be here anyway
    577       !!! but are kept for those who would want to revert to previous behaviour
    578       DO ij = 1, iim
    579         tppn(ij) = aire(ij) * ps (ij)
    580         tpps(ij) = aire(ij + ip1jm) * ps (ij + ip1jm)
    581       ENDDO
    582       tpn = SSUM(iim, tppn, 1) / apoln
    583       tps = SSUM(iim, tpps, 1) / apols
    584 
    585       DO ij = 1, iip1
    586         ps(ij) = tpn
    587         ps(ij + ip1jm) = tps
    588       ENDDO
    589     endif ! of if (1 == 0)
    590 
    591   END IF ! of IF(apdiss)
    592 
    593   ! ajout debug
    594   ! IF( lafin ) THEN
    595   !   abort_message = 'Simulation finished'
    596   !   CALL abort_gcm(modname,abort_message,0)
    597   ! ENDIF
    598 
    599   !   ********************************************************************
    600   !   ********************************************************************
    601   !   .... fin de l'integration dynamique  et physique pour le pas itau ..
    602   !   ********************************************************************
    603   !   ********************************************************************
    604 
    605   !   preparation du pas d'integration suivant  ......
    606 
    607   CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1509')
    608 
    609   IF (.NOT.purmats) THEN
    610     ! ........................................................
    611     ! ..............  schema matsuno + leapfrog  ..............
    612     ! ........................................................
    613 
    614     IF(forward .OR. leapf) THEN
    615       itau = itau + 1
    616       ! iday= day_ini+itau/day_step
    617       ! time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
    618       !   IF(time.GT.1.) THEN
    619       !     time = time-1.
    620       !     iday = iday+1
    621       !   ENDIF
    622     ENDIF
    623 
    624     IF(itau == itaufinp1) THEN
    625       IF (flag_verif) THEN
    626         WRITE(79, *) 'ucov', ucov
    627         WRITE(80, *) 'vcov', vcov
    628         WRITE(81, *) 'teta', teta
    629         WRITE(82, *) 'ps', ps
    630         WRITE(83, *) 'q', q
    631         WRITE(85, *) 'q1 = ', q(:, :, 1)
    632         WRITE(86, *) 'q3 = ', q(:, :, 3)
    633       endif
    634 
    635       abort_message = 'Simulation finished'
    636 
    637       CALL abort_gcm(modname, abort_message, 0)
    638     ENDIF
    639     !-----------------------------------------------------------------------
    640     !   ecriture du fichier histoire moyenne:
    641     !   -------------------------------------
    642 
    643     IF(MOD(itau, iperiod)==0 .OR. itau==itaufin) THEN
    644       IF(itau==itaufin) THEN
    645         iav = 1
    646       ELSE
    647         iav = 0
    648       ENDIF
    649 
    650       ! Ehouarn: re-compute geopotential for outputs
    651       CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
    652 
    653       IF (ok_dynzon) THEN
    654         CALL bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &
    655                 ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
    656       END IF
    657       IF (ok_dyn_ave) THEN
    658         CALL writedynav(itau, vcov, &
    659                 ucov, teta, pk, phi, q, masse, ps, phis)
    660       ENDIF
    661 
    662     ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
    663 
    664     CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1584')
    665 
    666     !-----------------------------------------------------------------------
    667     !   ecriture de la bande histoire:
    668     !   ------------------------------
    669 
    670     IF(MOD(itau, iecri)==0) THEN
    671       ! ! Ehouarn: output only during LF or Backward Matsuno
    672       IF (leapf.OR.(.NOT.leapf.AND.(.NOT.forward))) THEN
    673         CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
    674         unat = 0.
    675         DO l = 1, llm
    676           unat(iip2:ip1jm, l) = ucov(iip2:ip1jm, l) / cu(iip2:ip1jm)
    677           vnat(:, l) = vcov(:, l) / cv(:)
    678         enddo
    679         IF (ok_dyn_ins) THEN
    680           ! WRITE(lunout,*) "leapfrog: CALL writehist, itau=",itau
    681           CALL writehist(itau, vcov, ucov, teta, phi, q, masse, ps, phis)
    682           ! CALL WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
    683           ! CALL WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
    684           ! CALL WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
    685           !  CALL WriteField('ps',reshape(ps,(/iip1,jmp1/)))
    686           !  CALL WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
    687         endif ! of if (ok_dyn_ins)
    688         ! For some Grads outputs of fields
    689         IF (output_grads_dyn) THEN
    690           INCLUDE "write_grads_dyn.h"
     631
     632      IF(itau == itaufinp1) THEN
     633        IF (flag_verif) THEN
     634          WRITE(79, *) 'ucov', ucov
     635          WRITE(80, *) 'vcov', vcov
     636          WRITE(81, *) 'teta', teta
     637          WRITE(82, *) 'ps', ps
     638          WRITE(83, *) 'q', q
     639          WRITE(85, *) 'q1 = ', q(:, :, 1)
     640          WRITE(86, *) 'q3 = ', q(:, :, 3)
    691641        endif
    692       endif ! of if (leapf.OR.(.NOT.leapf.AND.(.NOT.forward)))
    693     ENDIF ! of IF(MOD(itau,iecri).EQ.0)
    694 
    695     IF(itau==itaufin) THEN
    696 
    697 
    698       ! if (planet_type.EQ."earth") THEN
    699       ! Write an Earth-format restart file
    700       CALL dynredem1("restart.nc", start_time, &
    701               vcov, ucov, teta, q, masse, ps)
    702       ! END IF ! of if (planet_type.EQ."earth")
    703 
    704       CLOSE(99)
    705       IF (ok_guide) THEN
    706         ! ! set ok_guide to false to avoid extra output
    707         ! ! in following forward step
    708         ok_guide = .FALSE.
    709       endif
    710       ! !!! Ehouarn: Why not stop here and now?
    711     ENDIF ! of IF (itau.EQ.itaufin)
    712 
    713     !-----------------------------------------------------------------------
    714     !   gestion de l'integration temporelle:
    715     !   ------------------------------------
    716 
    717     IF(MOD(itau, iperiod)==0)    THEN
    718       GO TO 1
    719     ELSE IF (MOD(itau - 1, iperiod) == 0) THEN
    720 
    721       IF(forward)  THEN
    722         ! fin du pas forward et debut du pas backward
    723 
    724         forward = .FALSE.
    725         leapf = .FALSE.
    726         GO TO 2
    727 
    728       ELSE
    729         ! fin du pas backward et debut du premier pas leapfrog
    730 
    731         leapf = .TRUE.
    732         dt = 2. * dtvr
    733         GO TO 2
    734       END IF ! of IF (forward)
    735     ELSE
    736 
    737       ! ......   pas leapfrog  .....
    738 
    739       leapf = .TRUE.
    740       dt = 2. * dtvr
    741       GO TO 2
    742     END IF ! of IF (MOD(itau,iperiod).EQ.0)
    743     ! !    ELSEIF (MOD(itau-1,iperiod).EQ.0)
    744 
    745   ELSE ! of IF (.NOT.purmats)
    746 
    747     CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1664')
    748 
    749     ! ........................................................
    750     ! ..............       schema  matsuno        ...............
    751     ! ........................................................
    752     IF(forward)  THEN
    753 
    754       itau = itau + 1
    755       ! iday = day_ini+itau/day_step
    756       ! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
    757 
    758       !              IF(time.GT.1.) THEN
    759       !               time = time-1.
    760       !               iday = iday+1
    761       !              ENDIF
    762 
    763       forward = .FALSE.
    764       IF(itau == itaufinp1) THEN
     642
    765643        abort_message = 'Simulation finished'
     644
    766645        CALL abort_gcm(modname, abort_message, 0)
    767646      ENDIF
    768       GO TO 2
    769 
    770     ELSE ! of IF(forward) i.e. backward step
    771 
    772       CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1698')
     647      !-----------------------------------------------------------------------
     648      !   ecriture du fichier histoire moyenne:
     649      !   -------------------------------------
    773650
    774651      IF(MOD(itau, iperiod)==0 .OR. itau==itaufin) THEN
    … …  
    779656        ENDIF
    780657
    781         ! ! Ehouarn: re-compute geopotential for outputs
     658        ! Ehouarn: re-compute geopotential for outputs
    782659        CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
    783660
    … …  
    785662          CALL bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &
    786663                  ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
    787         ENDIF
     664        END IF
    788665        IF (ok_dyn_ave) THEN
    789666          CALL writedynav(itau, vcov, &
    … …  
    791668        ENDIF
    792669
    793       ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
     670      ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
     671
     672      CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1584')
     673
     674      !-----------------------------------------------------------------------
     675      !   ecriture de la bande histoire:
     676      !   ------------------------------
    794677
    795678      IF(MOD(itau, iecri)==0) THEN
    796         ! IF(MOD(itau,iecri*day_step).EQ.0) THEN
    797         CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
    798         unat = 0.
    799         DO l = 1, llm
    800           unat(iip2:ip1jm, l) = ucov(iip2:ip1jm, l) / cu(iip2:ip1jm)
    801           vnat(:, l) = vcov(:, l) / cv(:)
    802         enddo
    803         IF (ok_dyn_ins) THEN
    804           ! WRITE(lunout,*) "leapfrog: CALL writehist (b)",
    805           ! &                        itau,iecri
    806           CALL writehist(itau, vcov, ucov, teta, phi, q, masse, ps, phis)
    807         endif ! of if (ok_dyn_ins)
    808         ! For some Grads outputs
    809         IF (output_grads_dyn) THEN
    810           INCLUDE "write_grads_dyn.h"
    811         endif
    812 
    813       ENDIF ! of IF(MOD(itau,iecri         ).EQ.0)
     679        ! ! Ehouarn: output only during LF or Backward Matsuno
     680        IF (leapf.OR.(.NOT.leapf.AND.(.NOT.forward))) THEN
     681          CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
     682          unat = 0.
     683          DO l = 1, llm
     684            unat(iip2:ip1jm, l) = ucov(iip2:ip1jm, l) / cu(iip2:ip1jm)
     685            vnat(:, l) = vcov(:, l) / cv(:)
     686          enddo
     687          IF (ok_dyn_ins) THEN
     688            ! WRITE(lunout,*) "leapfrog: CALL writehist, itau=",itau
     689            CALL writehist(itau, vcov, ucov, teta, phi, q, masse, ps, phis)
     690            ! CALL WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
     691            ! CALL WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
     692            ! CALL WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
     693            !  CALL WriteField('ps',reshape(ps,(/iip1,jmp1/)))
     694            !  CALL WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
     695          endif ! of if (ok_dyn_ins)
     696          ! For some Grads outputs of fields
     697          IF (output_grads_dyn) THEN
     698            INCLUDE "write_grads_dyn.h"
     699          endif
     700        endif ! of if (leapf.OR.(.NOT.leapf.AND.(.NOT.forward)))
     701      ENDIF ! of IF(MOD(itau,iecri).EQ.0)
    814702
    815703      IF(itau==itaufin) THEN
     704
     705
    816706        ! if (planet_type.EQ."earth") THEN
     707        ! Write an Earth-format restart file
    817708        CALL dynredem1("restart.nc", start_time, &
    818709                vcov, ucov, teta, q, masse, ps)
    819710        ! END IF ! of if (planet_type.EQ."earth")
     711
     712        CLOSE(99)
    820713        IF (ok_guide) THEN
    821714          ! ! set ok_guide to false to avoid extra output
    … …  
    823716          ok_guide = .FALSE.
    824717        endif
    825       ENDIF ! of IF(itau.EQ.itaufin)
    826 
    827       forward = .TRUE.
    828       GO TO  1
    829 
    830     ENDIF ! of IF (forward)
    831 
    832   END IF ! of IF(.NOT.purmats)
    833 
    834 END SUBROUTINE leapfrog
     718        ! !!! Ehouarn: Why not stop here and now?
     719      ENDIF ! of IF (itau.EQ.itaufin)
     720
     721      !-----------------------------------------------------------------------
     722      !   gestion de l'integration temporelle:
     723      !   ------------------------------------
     724
     725      IF(MOD(itau, iperiod)==0)    THEN
     726        GO TO 1
     727      ELSE IF (MOD(itau - 1, iperiod) == 0) THEN
     728
     729        IF(forward)  THEN
     730          ! fin du pas forward et debut du pas backward
     731
     732          forward = .FALSE.
     733          leapf = .FALSE.
     734          GO TO 2
     735
     736        ELSE
     737          ! fin du pas backward et debut du premier pas leapfrog
     738
     739          leapf = .TRUE.
     740          dt = 2. * dtvr
     741          GO TO 2
     742        END IF ! of IF (forward)
     743      ELSE
     744
     745        ! ......   pas leapfrog  .....
     746
     747        leapf = .TRUE.
     748        dt = 2. * dtvr
     749        GO TO 2
     750      END IF ! of IF (MOD(itau,iperiod).EQ.0)
     751      ! !    ELSEIF (MOD(itau-1,iperiod).EQ.0)
     752
     753    ELSE ! of IF (.NOT.purmats)
     754
     755      CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1664')
     756
     757      ! ........................................................
     758      ! ..............       schema  matsuno        ...............
     759      ! ........................................................
     760      IF(forward)  THEN
     761
     762        itau = itau + 1
     763        ! iday = day_ini+itau/day_step
     764        ! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
     765
     766        !              IF(time.GT.1.) THEN
     767        !               time = time-1.
     768        !               iday = iday+1
     769        !              ENDIF
     770
     771        forward = .FALSE.
     772        IF(itau == itaufinp1) THEN
     773          abort_message = 'Simulation finished'
     774          CALL abort_gcm(modname, abort_message, 0)
     775        ENDIF
     776        GO TO 2
     777
     778      ELSE ! of IF(forward) i.e. backward step
     779
     780        CALL check_isotopes_seq(q, ip1jmp1, 'leapfrog 1698')
     781
     782        IF(MOD(itau, iperiod)==0 .OR. itau==itaufin) THEN
     783          IF(itau==itaufin) THEN
     784            iav = 1
     785          ELSE
     786            iav = 0
     787          ENDIF
     788
     789          ! ! Ehouarn: re-compute geopotential for outputs
     790          CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
     791
     792          IF (ok_dynzon) THEN
     793            CALL bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &
     794                    ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
     795          ENDIF
     796          IF (ok_dyn_ave) THEN
     797            CALL writedynav(itau, vcov, &
     798                    ucov, teta, pk, phi, q, masse, ps, phis)
     799          ENDIF
     800
     801        ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
     802
     803        IF(MOD(itau, iecri)==0) THEN
     804          ! IF(MOD(itau,iecri*day_step).EQ.0) THEN
     805          CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
     806          unat = 0.
     807          DO l = 1, llm
     808            unat(iip2:ip1jm, l) = ucov(iip2:ip1jm, l) / cu(iip2:ip1jm)
     809            vnat(:, l) = vcov(:, l) / cv(:)
     810          enddo
     811          IF (ok_dyn_ins) THEN
     812            ! WRITE(lunout,*) "leapfrog: CALL writehist (b)",
     813            ! &                        itau,iecri
     814            CALL writehist(itau, vcov, ucov, teta, phi, q, masse, ps, phis)
     815          endif ! of if (ok_dyn_ins)
     816          ! For some Grads outputs
     817          IF (output_grads_dyn) THEN
     818            INCLUDE "write_grads_dyn.h"
     819          endif
     820
     821        ENDIF ! of IF(MOD(itau,iecri         ).EQ.0)
     822
     823        IF(itau==itaufin) THEN
     824          ! if (planet_type.EQ."earth") THEN
     825          CALL dynredem1("restart.nc", start_time, &
     826                  vcov, ucov, teta, q, masse, ps)
     827          ! END IF ! of if (planet_type.EQ."earth")
     828          IF (ok_guide) THEN
     829            ! ! set ok_guide to false to avoid extra output
     830            ! ! in following forward step
     831            ok_guide = .FALSE.
     832          endif
     833        ENDIF ! of IF(itau.EQ.itaufin)
     834
     835        forward = .TRUE.
     836        GO TO  1
     837
     838      ENDIF ! of IF (forward)
     839
     840    END IF ! of IF(.NOT.purmats)
     841
     842  END SUBROUTINE leapfrog
     843
     844END MODULE lmdz_leapfrog
Note: See TracChangeset for help on using the changeset viewer.