Ignore:
Timestamp:
Jul 24, 2024, 4:39:59 PM (6 months ago)
Author:
abarral
Message:

Replace iniprint.h by lmdz_iniprint.f90
(lint) along the way

Location:
LMDZ6/branches/Amaury_dev/libf/dyn3d
Files:
16 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/abort_gcm.F90

    r5117 r5118  
    88  !! ug Pour les sorties XIOS
    99  USE lmdz_wxios
    10 
    11   include "iniprint.h"
     10  USE lmdz_iniprint, ONLY: lunout, prt_level
    1211
    1312  !

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

    r5117 r5118  
    1616  USE lmdz_description, ONLY: descript
    1717  USE lmdz_libmath, ONLY: minmax
     18  USE lmdz_iniprint, ONLY: lunout, prt_level
    1819
    1920  IMPLICIT NONE
    … …  
    2324  include "comdissip.h"
    2425  include "comgeom2.h"
    25   include "iniprint.h"
    2626
    2727  !---------------------------------------------------------------------------

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/bilan_dyn.F90

    r5117 r5118  
    1313  USE comvert_mod, ONLY: presnivs
    1414  USE temps_mod, ONLY: annee_ref, day_ref, itau_dyn
     15  USE lmdz_iniprint, ONLY: lunout, prt_level
    1516
    1617  IMPLICIT NONE
    … …  
    1920  include "paramet.h"
    2021  include "comgeom2.h"
    21   include "iniprint.h"
    2222
    2323  !====================================================================

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

    r5117 r5118  
    1616          alphax, alphay, taux, tauy
    1717  USE temps_mod, ONLY: calend, year_len
     18  USE lmdz_iniprint, ONLY: lunout, prt_level
    1819
    1920  IMPLICIT NONE
    … …  
    3536  include "paramet.h"
    3637  include "comdissnew.h"
    37   include "iniprint.h"
    3838
    3939  !   local:

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/dynetat0.F90

    r5117 r5118  
    2121  USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_INCA
    2222  USE lmdz_description, ONLY: descript
     23  USE lmdz_iniprint, ONLY: lunout, prt_level
    2324
    2425  IMPLICIT NONE
    … …  
    2627  include "paramet.h"
    2728  include "comgeom2.h"
    28   include "iniprint.h"
    2929!===============================================================================
    3030! Arguments:

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/dynredem.F90

    r5117 r5118  
    1 SUBROUTINE dynredem0(fichnom,iday_end,phis)
    2 
    3 !-------------------------------------------------------------------------------
    4 ! Write the NetCDF restart file (initialization).
    5 !-------------------------------------------------------------------------------
     1SUBROUTINE dynredem0(fichnom, iday_end, phis)
     2
     3  !-------------------------------------------------------------------------------
     4  ! Write the NetCDF restart file (initialization).
     5  !-------------------------------------------------------------------------------
    66  USE IOIPSL
    77  USE lmdz_strings, ONLY: maxlen
    88  USE infotrac, ONLY: nqtot, tracers
    9   USE netcdf, ONLY: nf90_create, nf90_def_dim, nf90_inq_varid, nf90_global,    &
    10                     nf90_close,  nf90_put_att, nf90_unlimited, nf90_clobber,  &
    11                     nf90_64bit_offset
     9  USE netcdf, ONLY: nf90_create, nf90_def_dim, nf90_inq_varid, nf90_global, &
     10          nf90_close, nf90_put_att, nf90_unlimited, nf90_clobber, &
     11          nf90_64bit_offset
    1212  USE dynredem_mod, ONLY: cre_var, put_var1, put_var2, err, modname, fil
    13   USE comvert_mod,  ONLY: ap, bp, presnivs, pa, preff, nivsig, nivsigs
     13  USE comvert_mod, ONLY: ap, bp, presnivs, pa, preff, nivsig, nivsigs
    1414  USE comconst_mod, ONLY: cpp, daysec, dtvr, g, kappa, omeg, rad
    1515  USE logic_mod, ONLY: fxyhypb, ysinus
    16   USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy, &
    17                               taux,tauy
     16  USE serre_mod, ONLY: clon, clat, grossismx, grossismy, dzoomx, dzoomy, &
     17          taux, tauy
    1818  USE temps_mod, ONLY: annee_ref, day_ref, itau_dyn, itaufin, start_time
    19   USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
     19  USE ener_mod, ONLY: etot0, ptot0, ztot0, stot0, ang0
    2020  USE lmdz_description, ONLY: descript
    21  
     21  USE lmdz_iniprint, ONLY: lunout, prt_level
     22
    2223  IMPLICIT NONE
    2324  include "dimensions.h"
    2425  include "paramet.h"
    2526  include "comgeom2.h"
    26   include "iniprint.h"
    27 !===============================================================================
    28 ! Arguments:
    29   CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
    30   INTEGER,          INTENT(IN) :: iday_end         !---
    31   REAL,             INTENT(IN) :: phis(iip1, jjp1) !--- GROUND GEOPOTENTIAL
    32 !===============================================================================
    33 ! Local variables:
     27  !===============================================================================
     28  ! Arguments:
     29  CHARACTER(LEN = *), INTENT(IN) :: fichnom          !--- FILE NAME
     30  INTEGER, INTENT(IN) :: iday_end         !---
     31  REAL, INTENT(IN) :: phis(iip1, jjp1) !--- GROUND GEOPOTENTIAL
     32  !===============================================================================
     33  ! Local variables:
    3434  INTEGER :: iq
    35   INTEGER, PARAMETER :: length=100
    36   REAL    :: tab_cntrl(length)                     !--- RUN PARAMETERS TABLE
    37 !   For NetCDF:
    38   CHARACTER(LEN=maxlen) :: unites
     35  INTEGER, PARAMETER :: length = 100
     36  REAL :: tab_cntrl(length)                     !--- RUN PARAMETERS TABLE
     37  !   For NetCDF:
     38  CHARACTER(LEN = maxlen) :: unites
    3939  INTEGER :: indexID
    4040  INTEGER :: rlonuID, rlonvID, rlatuID, rlatvID
    4141  INTEGER :: sID, sigID, nID, timID
    4242  INTEGER :: yyears0, jjour0, mmois0
    43   REAL    :: zjulian, hours
    44 !===============================================================================
    45   modname='dynredem0'; fil=fichnom
     43  REAL :: zjulian, hours
     44  !===============================================================================
     45  modname = 'dynredem0'; fil = fichnom
    4646  CALL ymds2ju(annee_ref, 1, iday_end, 0.0, zjulian)
    4747  CALL ju2ymds(zjulian, yyears0, mmois0, jjour0, hours)
    4848
    49   tab_cntrl(:)  = 0.
    50   tab_cntrl(1)  = REAL(iim)
    51   tab_cntrl(2)  = REAL(jjm)
    52   tab_cntrl(3)  = REAL(llm)
    53   tab_cntrl(4)  = REAL(day_ref)
    54   tab_cntrl(5)  = REAL(annee_ref)
    55   tab_cntrl(6)  = rad
    56   tab_cntrl(7)  = omeg
    57   tab_cntrl(8)  = g
    58   tab_cntrl(9)  = cpp
     49  tab_cntrl(:) = 0.
     50  tab_cntrl(1) = REAL(iim)
     51  tab_cntrl(2) = REAL(jjm)
     52  tab_cntrl(3) = REAL(llm)
     53  tab_cntrl(4) = REAL(day_ref)
     54  tab_cntrl(5) = REAL(annee_ref)
     55  tab_cntrl(6) = rad
     56  tab_cntrl(7) = omeg
     57  tab_cntrl(8) = g
     58  tab_cntrl(9) = cpp
    5959  tab_cntrl(10) = kappa
    6060  tab_cntrl(11) = daysec
    … …  
    6868  tab_cntrl(19) = preff
    6969
    70 !    .....    parameters for zoom    ......   
     70  !    .....    parameters for zoom    ......
    7171  tab_cntrl(20) = clon
    7272  tab_cntrl(21) = clat
    … …  
    7474  tab_cntrl(23) = grossismy
    7575
    76   IF ( fxyhypb )   THEN
     76  IF (fxyhypb)   THEN
    7777    tab_cntrl(24) = 1.
    7878    tab_cntrl(25) = dzoomx
    … …  
    8888    tab_cntrl(28) = 0.
    8989    tab_cntrl(29) = 0.
    90     IF( ysinus )  tab_cntrl(27) = 1.
     90    IF(ysinus)  tab_cntrl(27) = 1.
    9191  END IF
    9292  tab_cntrl(30) = REAL(iday_end)
    9393  tab_cntrl(31) = REAL(itau_dyn + itaufin)
    94 ! start_time: start_time of simulation (not necessarily 0.)
     94  ! start_time: start_time of simulation (not necessarily 0.)
    9595  tab_cntrl(32) = start_time
    9696
    97 !--- File creation
    98   CALL err(nf90_create(fichnom,IOR(nf90_clobber,nf90_64bit_offset),nid))
    99 
    100 !--- Some global attributes
    101   CALL err(nf90_put_att(nid,nf90_global,"title","Fichier demarrage dynamique"))
    102 
    103 !--- Dimensions
    104   CALL err(nf90_def_dim(nid,"index", length, indexID))
    105   CALL err(nf90_def_dim(nid,"rlonu", iip1,  rlonuID))
    106   CALL err(nf90_def_dim(nid,"rlatu", jjp1,  rlatuID))
    107   CALL err(nf90_def_dim(nid,"rlonv", iip1,  rlonvID))
    108   CALL err(nf90_def_dim(nid,"rlatv", jjm,    rlatvID))
    109   CALL err(nf90_def_dim(nid,"sigs",  llm,        sID))
    110   CALL err(nf90_def_dim(nid,"sig",   llmp1,    sigID))
    111   CALL err(nf90_def_dim(nid,"temps", nf90_unlimited, timID))
    112 
    113 !--- Define and save invariant fields
    114   CALL put_var1(nid,"controle","Parametres de controle" ,[indexID],tab_cntrl)
    115   CALL put_var1(nid,"rlonu"   ,"Longitudes des points U",[rlonuID],rlonu)
    116   CALL put_var1(nid,"rlatu"   ,"Latitudes des points U" ,[rlatuID],rlatu)
    117   CALL put_var1(nid,"rlonv"   ,"Longitudes des points V",[rlonvID],rlonv)
    118   CALL put_var1(nid,"rlatv"   ,"Latitudes des points V" ,[rlatvID],rlatv)
    119   CALL put_var1(nid,"nivsigs" ,"Numero naturel des couches s"    ,[sID]  ,nivsigs)
    120   CALL put_var1(nid,"nivsig"  ,"Numero naturel des couches sigma",[sigID],nivsig)
    121   CALL put_var1(nid,"ap"      ,"Coefficient A pour hybride"      ,[sigID],ap)
    122   CALL put_var1(nid,"bp"      ,"Coefficient B pour hybride"      ,[sigID],bp)
    123   CALL put_var1(nid,"presnivs",""                                ,[sID]  ,presnivs)
    124 ! covariant <-> contravariant <-> natural conversion coefficients
    125   CALL put_var2(nid,"cu","Coefficient de passage pour U",[rlonuID,rlatuID],cu)
    126   CALL put_var2(nid,"cv","Coefficient de passage pour V",[rlonvID,rlatvID],cv)
    127   CALL put_var2(nid,"aire","Aires de chaque maille"     ,[rlonvID,rlatuID],aire)
    128   CALL put_var2(nid,"phisinit","Geopotentiel au sol"    ,[rlonvID,rlatuID],phis)
    129 
    130 !--- Define fields saved later
    131   WRITE(unites,"('days since ',i4,'-',i2.2,'-',i2.2,' 00:00:00')") &
    132                yyears0,mmois0,jjour0
    133   CALL cre_var(nid,"temps","Temps de simulation",[timID],unites)
    134   CALL cre_var(nid,"ucov" ,"Vitesse U"  ,[rlonuID,rlatuID,sID,timID])
    135   CALL cre_var(nid,"vcov" ,"Vitesse V"  ,[rlonvID,rlatvID,sID,timID])
    136   CALL cre_var(nid,"teta" ,"Temperature",[rlonvID,rlatuID,sID,timID])
    137   DO iq=1,nqtot
    138     CALL cre_var(nid,tracers(iq)%name,tracers(iq)%longName,[rlonvID,rlatuID,sID,timID])
     97  !--- File creation
     98  CALL err(nf90_create(fichnom, IOR(nf90_clobber, nf90_64bit_offset), nid))
     99
     100  !--- Some global attributes
     101  CALL err(nf90_put_att(nid, nf90_global, "title", "Fichier demarrage dynamique"))
     102
     103  !--- Dimensions
     104  CALL err(nf90_def_dim(nid, "index", length, indexID))
     105  CALL err(nf90_def_dim(nid, "rlonu", iip1, rlonuID))
     106  CALL err(nf90_def_dim(nid, "rlatu", jjp1, rlatuID))
     107  CALL err(nf90_def_dim(nid, "rlonv", iip1, rlonvID))
     108  CALL err(nf90_def_dim(nid, "rlatv", jjm, rlatvID))
     109  CALL err(nf90_def_dim(nid, "sigs", llm, sID))
     110  CALL err(nf90_def_dim(nid, "sig", llmp1, sigID))
     111  CALL err(nf90_def_dim(nid, "temps", nf90_unlimited, timID))
     112
     113  !--- Define and save invariant fields
     114  CALL put_var1(nid, "controle", "Parametres de controle", [indexID], tab_cntrl)
     115  CALL put_var1(nid, "rlonu", "Longitudes des points U", [rlonuID], rlonu)
     116  CALL put_var1(nid, "rlatu", "Latitudes des points U", [rlatuID], rlatu)
     117  CALL put_var1(nid, "rlonv", "Longitudes des points V", [rlonvID], rlonv)
     118  CALL put_var1(nid, "rlatv", "Latitudes des points V", [rlatvID], rlatv)
     119  CALL put_var1(nid, "nivsigs", "Numero naturel des couches s", [sID], nivsigs)
     120  CALL put_var1(nid, "nivsig", "Numero naturel des couches sigma", [sigID], nivsig)
     121  CALL put_var1(nid, "ap", "Coefficient A pour hybride", [sigID], ap)
     122  CALL put_var1(nid, "bp", "Coefficient B pour hybride", [sigID], bp)
     123  CALL put_var1(nid, "presnivs", "", [sID], presnivs)
     124  ! covariant <-> contravariant <-> natural conversion coefficients
     125  CALL put_var2(nid, "cu", "Coefficient de passage pour U", [rlonuID, rlatuID], cu)
     126  CALL put_var2(nid, "cv", "Coefficient de passage pour V", [rlonvID, rlatvID], cv)
     127  CALL put_var2(nid, "aire", "Aires de chaque maille", [rlonvID, rlatuID], aire)
     128  CALL put_var2(nid, "phisinit", "Geopotentiel au sol", [rlonvID, rlatuID], phis)
     129
     130  !--- Define fields saved later
     131  WRITE(unites, "('days since ',i4,'-',i2.2,'-',i2.2,' 00:00:00')") &
     132          yyears0, mmois0, jjour0
     133  CALL cre_var(nid, "temps", "Temps de simulation", [timID], unites)
     134  CALL cre_var(nid, "ucov", "Vitesse U", [rlonuID, rlatuID, sID, timID])
     135  CALL cre_var(nid, "vcov", "Vitesse V", [rlonvID, rlatvID, sID, timID])
     136  CALL cre_var(nid, "teta", "Temperature", [rlonvID, rlatuID, sID, timID])
     137  DO iq = 1, nqtot
     138    CALL cre_var(nid, tracers(iq)%name, tracers(iq)%longName, [rlonvID, rlatuID, sID, timID])
    139139  END DO
    140   CALL cre_var(nid,"masse","Masse d air"    ,[rlonvID,rlatuID,sID,timID])
    141   CALL cre_var(nid,"ps"   ,"Pression au sol",[rlonvID,rlatuID    ,timID])
     140  CALL cre_var(nid, "masse", "Masse d air", [rlonvID, rlatuID, sID, timID])
     141  CALL cre_var(nid, "ps", "Pression au sol", [rlonvID, rlatuID, timID])
    142142  CALL err(nf90_close (nid))
    143143
    144   WRITE(lunout,*)TRIM(modname)//': iim,jjm,llm,iday_end',iim,jjm,llm,iday_end
    145   WRITE(lunout,*)TRIM(modname)//': rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
     144  WRITE(lunout, *)TRIM(modname) // ': iim,jjm,llm,iday_end', iim, jjm, llm, iday_end
     145  WRITE(lunout, *)TRIM(modname) // ': rad,omeg,g,cpp,kappa', rad, omeg, g, cpp, kappa
    146146
    147147END SUBROUTINE dynredem0
    … …  
    152152!-------------------------------------------------------------------------------
    153153
    154 SUBROUTINE dynredem1(fichnom,time,vcov,ucov,teta,q,masse,ps)
    155 
    156 !-------------------------------------------------------------------------------
    157 ! Purpose: Write the NetCDF restart file (append).
    158 !-------------------------------------------------------------------------------
     154SUBROUTINE dynredem1(fichnom, time, vcov, ucov, teta, q, masse, ps)
     155
     156  !-------------------------------------------------------------------------------
     157  ! Purpose: Write the NetCDF restart file (append).
     158  !-------------------------------------------------------------------------------
    159159  USE lmdz_strings, ONLY: maxlen
    160160  USE infotrac, ONLY: nqtot, tracers, type_trac
    161161  USE control_mod
    162   USE netcdf,   ONLY: nf90_open,  nf90_nowrite, nf90_get_var, nf90_inq_varid, &
    163                       nf90_close, nf90_write,  nf90_put_var, nf90_noerr
     162  USE netcdf, ONLY: nf90_open, nf90_nowrite, nf90_get_var, nf90_inq_varid, &
     163          nf90_close, nf90_write, nf90_put_var, nf90_noerr
    164164  USE dynredem_mod, ONLY: dynredem_write_u, dynredem_write_v, dynredem_read_u, &
    165                           err, modname, fil, msg
     165          err, modname, fil, msg
    166166  USE temps_mod, ONLY: itau_dyn, itaufin
    167167  USE lmdz_description, ONLY: descript
    168  
     168  USE lmdz_iniprint, ONLY: lunout, prt_level
     169
    169170  IMPLICIT NONE
    170171  include "dimensions.h"
    171172  include "paramet.h"
    172173  include "comgeom.h"
    173   include "iniprint.h"
    174 !===============================================================================
    175 ! Arguments:
    176   CHARACTER(LEN=*), INTENT(IN) :: fichnom              !-- FILE NAME
    177   REAL, INTENT(IN)    ::  time                         !-- TIME
    178   REAL, INTENT(IN)    ::  vcov(iip1,jjm, llm)          !-- V COVARIANT WIND
    179   REAL, INTENT(IN)    ::  ucov(iip1,jjp1,llm)          !-- U COVARIANT WIND
    180   REAL, INTENT(IN)    ::  teta(iip1,jjp1,llm)          !-- POTENTIAL TEMPERATURE
    181   REAL, INTENT(INOUT) ::     q(iip1,jjp1,llm,nqtot)    !-- TRACERS
    182   REAL, INTENT(IN)    :: masse(iip1,jjp1,llm)          !-- MASS PER CELL
    183   REAL, INTENT(IN)    ::    ps(iip1,jjp1)              !-- GROUND PRESSURE
    184 !===============================================================================
    185 ! Local variables:
     174  !===============================================================================
     175  ! Arguments:
     176  CHARACTER(LEN = *), INTENT(IN) :: fichnom              !-- FILE NAME
     177  REAL, INTENT(IN) :: time                         !-- TIME
     178  REAL, INTENT(IN) :: vcov(iip1, jjm, llm)          !-- V COVARIANT WIND
     179  REAL, INTENT(IN) :: ucov(iip1, jjp1, llm)          !-- U COVARIANT WIND
     180  REAL, INTENT(IN) :: teta(iip1, jjp1, llm)          !-- POTENTIAL TEMPERATURE
     181  REAL, INTENT(INOUT) :: q(iip1, jjp1, llm, nqtot)    !-- TRACERS
     182  REAL, INTENT(IN) :: masse(iip1, jjp1, llm)          !-- MASS PER CELL
     183  REAL, INTENT(IN) :: ps(iip1, jjp1)              !-- GROUND PRESSURE
     184  !===============================================================================
     185  ! Local variables:
    186186  INTEGER :: iq, nid, vID, ierr, nid_trac, vID_trac
    187   INTEGER, SAVE :: nb=0
    188   INTEGER, PARAMETER :: length=100
    189   REAL               :: tab_cntrl(length) ! tableau des parametres du run
    190   CHARACTER(LEN=maxlen) :: var, dum
    191   LOGICAL            :: lread_inca
    192 !===============================================================================
    193 
    194   modname='dynredem1'; fil=fichnom
    195   CALL err(nf90_open(fil,nf90_write,nid),"open",fil)
    196 
    197 !--- Write/extend time coordinate
     187  INTEGER, SAVE :: nb = 0
     188  INTEGER, PARAMETER :: length = 100
     189  REAL :: tab_cntrl(length) ! tableau des parametres du run
     190  CHARACTER(LEN = maxlen) :: var, dum
     191  LOGICAL :: lread_inca
     192  !===============================================================================
     193
     194  modname = 'dynredem1'; fil = fichnom
     195  CALL err(nf90_open(fil, nf90_write, nid), "open", fil)
     196
     197  !--- Write/extend time coordinate
    198198  nb = nb + 1
    199   var="temps"
    200   CALL err(nf90_inq_varid(nid,var,vID),"inq",var)
    201   CALL err(nf90_put_var(nid,vID,[time]),"put",var)
    202   WRITE(lunout,*)TRIM(modname)//": Saving for ", nb, time
    203 
    204 !--- Rewrite control table (itaufin undefined in dynredem0)
    205   var="controle"
    206   CALL err(nf90_inq_varid(nid,var,vID),"inq",var)
    207   CALL err(nf90_get_var(nid,vID,tab_cntrl),"get",var)
    208   tab_cntrl(31)=DBLE(itau_dyn + itaufin)
    209   CALL err(nf90_inq_varid(nid,var,vID),"inq",var)
    210   CALL err(nf90_put_var(nid,vID,tab_cntrl),"put",var)
    211 
    212 !--- Save fields
    213   CALL dynredem_write_u(nid,"ucov" ,ucov ,llm)
    214   CALL dynredem_write_v(nid,"vcov" ,vcov ,llm)
    215   CALL dynredem_write_u(nid,"teta" ,teta ,llm)
    216   CALL dynredem_write_u(nid,"masse",masse,llm)
    217   CALL dynredem_write_u(nid,"ps"   ,ps   ,1)
    218 
    219 !--- Tracers in file "start_trac.nc" (added by Anne)
    220   lread_inca=.FALSE.; fil="start_trac.nc"
    221   IF(ANY(type_trac == ['inca','inco'])) INQUIRE(FILE=fil,EXIST=lread_inca)
    222   IF(lread_inca) CALL err(nf90_open(fil,nf90_nowrite,nid_trac),"open")
    223 
    224 !--- Save tracers
    225   DO iq=1,nqtot; var=TRIM(tracers(iq)%name); ierr=-1
    226     IF(lread_inca) THEN                  !--- Possibly read from "start_trac.nc"
    227       fil="start_trac.nc"
    228       ierr=nf90_inq_varid(nid_trac,var,vID_trac)
    229       dum='inq'; IF(ierr==nf90_noerr) dum='fnd'
    230       WRITE(lunout,*)msg(dum,var)
    231 
    232 
    233       IF(ierr==nf90_noerr) CALL dynredem_read_u(nid_trac,var,q(:,:,:,iq),llm)
    234     END IF
    235     fil=fichnom
    236     CALL dynredem_write_u(nid,var,q(:,:,:,iq),llm)
     199  var = "temps"
     200  CALL err(nf90_inq_varid(nid, var, vID), "inq", var)
     201  CALL err(nf90_put_var(nid, vID, [time]), "put", var)
     202  WRITE(lunout, *)TRIM(modname) // ": Saving for ", nb, time
     203
     204  !--- Rewrite control table (itaufin undefined in dynredem0)
     205  var = "controle"
     206  CALL err(nf90_inq_varid(nid, var, vID), "inq", var)
     207  CALL err(nf90_get_var(nid, vID, tab_cntrl), "get", var)
     208  tab_cntrl(31) = DBLE(itau_dyn + itaufin)
     209  CALL err(nf90_inq_varid(nid, var, vID), "inq", var)
     210  CALL err(nf90_put_var(nid, vID, tab_cntrl), "put", var)
     211
     212  !--- Save fields
     213  CALL dynredem_write_u(nid, "ucov", ucov, llm)
     214  CALL dynredem_write_v(nid, "vcov", vcov, llm)
     215  CALL dynredem_write_u(nid, "teta", teta, llm)
     216  CALL dynredem_write_u(nid, "masse", masse, llm)
     217  CALL dynredem_write_u(nid, "ps", ps, 1)
     218
     219  !--- Tracers in file "start_trac.nc" (added by Anne)
     220  lread_inca = .FALSE.; fil = "start_trac.nc"
     221  IF(ANY(type_trac == ['inca', 'inco'])) INQUIRE(FILE = fil, EXIST = lread_inca)
     222  IF(lread_inca) CALL err(nf90_open(fil, nf90_nowrite, nid_trac), "open")
     223
     224  !--- Save tracers
     225  DO iq = 1, nqtot; var = TRIM(tracers(iq)%name); ierr = -1
     226  IF(lread_inca) THEN                  !--- Possibly read from "start_trac.nc"
     227    fil = "start_trac.nc"
     228    ierr = nf90_inq_varid(nid_trac, var, vID_trac)
     229    dum = 'inq'; IF(ierr==nf90_noerr) dum = 'fnd'
     230    WRITE(lunout, *)msg(dum, var)
     231
     232    IF(ierr==nf90_noerr) CALL dynredem_read_u(nid_trac, var, q(:, :, :, iq), llm)
     233  END IF
     234  fil = fichnom
     235  CALL dynredem_write_u(nid, var, q(:, :, :, iq), llm)
    237236  END DO
    238   CALL err(nf90_close(nid),"close")
    239   fil="start_trac.nc"
    240   IF(lread_inca) CALL err(nf90_close(nid_trac),"close")
     237  CALL err(nf90_close(nid), "close")
     238  fil = "start_trac.nc"
     239  IF(lread_inca) CALL err(nf90_close(nid_trac), "close")
    241240
    242241END SUBROUTINE dynredem1

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/fluxstokenc.F90

    r5117 r5118  
    66
    77  USE IOIPSL
     8  USE lmdz_iniprint, ONLY: lunout, prt_level
    89  !
    910  ! Auteur :  F. Hourdin
    … …  
    1819  include "comgeom.h"
    1920  include "tracstoke.h"
    20   include "iniprint.h"
    2121
    2222  REAL :: time_step, t_wrt, t_ops
    … …  
    157157  ENDIF ! if iadvtr.EQ.istdyn
    158158
    159 
    160159END SUBROUTINE fluxstokenc

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/friction.F90

    r5117 r5118  
    77  USE IOIPSL
    88  USE comconst_mod, ONLY: pi
     9  USE lmdz_iniprint, ONLY: lunout, prt_level
    910  IMPLICIT NONE
    1011
    … …  
    2425  include "paramet.h"
    2526  include "comgeom2.h"
    26   include "iniprint.h"
    2727  include "academic.h"
    2828
    … …  
    125125  ENDIF
    126126
    127 
    128127END SUBROUTINE friction
    129128

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/gcm.F90

    r5117 r5118  
    1 
    21! $Id$
    32
    … …  
    1312  USE control_mod
    1413  USE mod_const_mpi, ONLY: COMM_LMDZ
    15   USE temps_mod, ONLY: calend,start_time,annee_ref,day_ref, &
    16                      itau_dyn,itau_phy,day_ini,jD_ref,jH_ref,day_end
     14  USE temps_mod, ONLY: calend, start_time, annee_ref, day_ref, &
     15          itau_dyn, itau_phy, day_ini, jD_ref, jH_ref, day_end
    1716  USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, g, r, rad
    1817  USE logic_mod, ONLY: ecripar, iflag_phys, read_start
    … …  
    2120  USE lmdz_description, ONLY: descript
    2221
    23 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     22  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2423  ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
    2524  ! A nettoyer. On ne veut qu'une ou deux routines d'interface
    … …  
    2726  ! AB 2024/07/20: remplace CPP key by fortran logical, but ^ still relevant, see later use of iniphys later on
    2827  USE iniphysiq_mod, ONLY: iniphysiq
    29 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     28  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     29  USE lmdz_iniprint, ONLY: lunout, prt_level
    3030
    3131  IMPLICIT NONE
    … …  
    6565  include "comdissnew.h"
    6666  include "comgeom.h"
    67   include "iniprint.h"
    6867  include "tracstoke.h"
    6968
    … …  
    7170
    7271  !   variables dynamiques
    73   REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
    74   REAL teta(ip1jmp1,llm)                 ! temperature potentielle
    75   REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
     72  REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants
     73  REAL teta(ip1jmp1, llm)                 ! temperature potentielle
     74  REAL, ALLOCATABLE, DIMENSION(:, :, :) :: q! champs advectes
    7675  REAL ps(ip1jmp1)                       ! pression  au sol
    77 !  REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
    78   REAL masse(ip1jmp1,llm)                ! masse d'air
     76  !  REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
     77  REAL masse(ip1jmp1, llm)                ! masse d'air
    7978  REAL phis(ip1jmp1)                     ! geopotentiel au sol
    80 !  REAL phi(ip1jmp1,llm)                  ! geopotentiel
    81 !  REAL w(ip1jmp1,llm)                    ! vitesse verticale
     79  !  REAL phi(ip1jmp1,llm)                  ! geopotentiel
     80  !  REAL w(ip1jmp1,llm)                    ! vitesse verticale
    8281
    8382  ! variables dynamiques intermediaire pour le transport
    … …  
    8988
    9089  LOGICAL lafin
    91 
    9290
    9391  REAL time_step, t_wrt, t_ops
    … …  
    10199  !     tansformation d'energie cinetique en energie thermique
    102100  !     cree par la dissipation
    103 !  REAL dhecdt(ip1jmp1,llm)
     101  !  REAL dhecdt(ip1jmp1,llm)
    104102  !      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
    105103  !      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
    106 !  CHARACTER (len=15) :: ztit
     104  !  CHARACTER (len=15) :: ztit
    107105  !-jld
    108106
    109 
    110   CHARACTER (LEN=80) :: dynhist_file, dynhistave_file
    111   CHARACTER (LEN=20) :: modname
    112   CHARACTER (LEN=80) :: abort_message
     107  CHARACTER (LEN = 80) :: dynhist_file, dynhistave_file
     108  CHARACTER (LEN = 20) :: modname
     109  CHARACTER (LEN = 80) :: abort_message
    113110  ! locales pour gestion du temps
    114111  INTEGER :: an, mois, jour
    … …  
    123120  modname = 'gcm'
    124121  descript = 'Run GCM LMDZ'
    125   lafin    = .FALSE.
     122  lafin = .FALSE.
    126123  dynhist_file = 'dyn_hist.nc'
    127124  dynhistave_file = 'dyn_hist_ave.nc'
    … …  
    133130  !  ---------------------------------------
    134131
    135   CALL conf_gcm( 99, .TRUE.)
     132  CALL conf_gcm(99, .TRUE.)
    136133
    137134  IF (mod(iphysiq, iperiod) /= 0) CALL abort_gcm("conf_gcm", &
    138        "iphysiq must be a multiple of iperiod", 1)
    139 
    140   use_filtre_fft=.FALSE.
    141   CALL getin('use_filtre_fft',use_filtre_fft)
     135          "iphysiq must be a multiple of iperiod", 1)
     136
     137  use_filtre_fft = .FALSE.
     138  CALL getin('use_filtre_fft', use_filtre_fft)
    142139  IF (use_filtre_fft) CALL abort_gcm("gcm", 'FFT filter is not available in ' &
    143140          // 'the sequential version of the dynamics.', 1)
    144141
    145 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     142  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    146143  ! Initialisation de XIOS
    147 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     144  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    148145
    149146  IF (using_xios) THEN
    … …  
    152149
    153150
    154 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     151  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    155152  ! FH 2008/05/02
    156153  ! A nettoyer. On ne veut qu'une ou deux routines d'interface
    157154  ! dynamique -> physique pour l'initialisation
    158 !#ifdef CPP_PHYS
    159 !  CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
    160 !  !      CALL InitComgeomphy ! now done in iniphysiq
    161 !#endif
    162 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     155  !#ifdef CPP_PHYS
     156  !  CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
     157  !  !      CALL InitComgeomphy ! now done in iniphysiq
     158  !#endif
     159  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    163160  !-----------------------------------------------------------------------
    164161  !   Choix du calendrier
    … …  
    168165
    169166  IF (calend == 'earth_360d') THEN
    170      CALL ioconf_calendar('360_day')
    171      WRITE(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
     167    CALL ioconf_calendar('360_day')
     168    WRITE(lunout, *)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
    172169  ELSE IF (calend == 'earth_365d') THEN
    173      CALL ioconf_calendar('noleap')
    174      WRITE(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
     170    CALL ioconf_calendar('noleap')
     171    WRITE(lunout, *)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
    175172  ELSE IF (calend == 'gregorian') THEN
    176      CALL ioconf_calendar('gregorian')
    177      WRITE(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
     173    CALL ioconf_calendar('gregorian')
     174    WRITE(lunout, *)'CALENDRIER CHOISI: Terrestre bissextile'
    178175  else
    179      abort_message = 'Mauvais choix de calendrier'
    180      CALL abort_gcm(modname,abort_message,1)
     176    abort_message = 'Mauvais choix de calendrier'
     177    CALL abort_gcm(modname, abort_message, 1)
    181178  ENDIF
    182179  !-----------------------------------------------------------------------
    … …  
    196193
    197194  ! Allocation de la tableau q : champs advectes   
    198   allocate(q(ip1jmp1,llm,nqtot))
     195  allocate(q(ip1jmp1, llm, nqtot))
    199196
    200197  !-----------------------------------------------------------------------
    … …  
    204201  !  lecture du fichier start.nc
    205202  IF (read_start) THEN
    206      ! we still need to run iniacademic to initialize some
    207      ! constants & fields, if we run the 'newtonian' or 'SW' cases:
    208      IF (iflag_phys/=1) THEN
    209         CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
    210      endif
    211 
    212      !        if (planet_type.EQ."earth") THEN
    213      ! Load an Earth-format start file
    214      CALL dynetat0("start.nc",vcov,ucov, &
    215           teta,q,masse,ps,phis, time_0)
    216      !        endif ! of if (planet_type.EQ."earth")
    217 
    218      !       WRITE(73,*) 'ucov',ucov
    219      !       WRITE(74,*) 'vcov',vcov
    220      !       WRITE(75,*) 'teta',teta
    221      !       WRITE(76,*) 'ps',ps
    222      !       WRITE(77,*) 'q',q
     203    ! we still need to run iniacademic to initialize some
     204    ! constants & fields, if we run the 'newtonian' or 'SW' cases:
     205    IF (iflag_phys/=1) THEN
     206      CALL iniacademic(vcov, ucov, teta, q, masse, ps, phis, time_0)
     207    endif
     208
     209    !        if (planet_type.EQ."earth") THEN
     210    ! Load an Earth-format start file
     211    CALL dynetat0("start.nc", vcov, ucov, &
     212            teta, q, masse, ps, phis, time_0)
     213    !        endif ! of if (planet_type.EQ."earth")
     214
     215    !       WRITE(73,*) 'ucov',ucov
     216    !       WRITE(74,*) 'vcov',vcov
     217    !       WRITE(75,*) 'teta',teta
     218    !       WRITE(76,*) 'ps',ps
     219    !       WRITE(77,*) 'q',q
    223220
    224221  ENDIF ! of if (read_start)
    … …  
    226223
    227224  ! le cas echeant, creation d un etat initial
    228   IF (prt_level > 9) WRITE(lunout,*) &
    229        'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
     225  IF (prt_level > 9) WRITE(lunout, *) &
     226          'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
    230227  IF (.NOT.read_start) THEN
    231      start_time=0.
    232      annee_ref=anneeref
    233      CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
     228    start_time = 0.
     229    annee_ref = anneeref
     230    CALL iniacademic(vcov, ucov, teta, q, masse, ps, phis, time_0)
    234231  ENDIF
    235232
    … …  
    240237  !  on recalcule eventuellement le pas de temps
    241238
    242   IF(MOD(day_step,iperiod)/=0) THEN
    243      abort_message = &
    244           'Il faut choisir un nb de pas par jour multiple de iperiod'
    245      CALL abort_gcm(modname,abort_message,1)
    246   ENDIF
    247 
    248   IF(MOD(day_step,iphysiq)/=0) THEN
    249      abort_message = &
    250           'Il faut choisir un nb de pas par jour multiple de iphysiq'
    251      CALL abort_gcm(modname,abort_message,1)
    252   ENDIF
    253 
    254   zdtvr    = daysec/REAL(day_step)
     239  IF(MOD(day_step, iperiod)/=0) THEN
     240    abort_message = &
     241            'Il faut choisir un nb de pas par jour multiple de iperiod'
     242    CALL abort_gcm(modname, abort_message, 1)
     243  ENDIF
     244
     245  IF(MOD(day_step, iphysiq)/=0) THEN
     246    abort_message = &
     247            'Il faut choisir un nb de pas par jour multiple de iphysiq'
     248    CALL abort_gcm(modname, abort_message, 1)
     249  ENDIF
     250
     251  zdtvr = daysec / REAL(day_step)
    255252  IF(dtvr/=zdtvr) THEN
    256      WRITE(lunout,*) &
    257           'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
     253    WRITE(lunout, *) &
     254            'WARNING!!! changement de pas de temps', dtvr, '>', zdtvr
    258255  ENDIF
    259256
    … …  
    261258
    262259  IF (start_time /= starttime) THEN
    263      WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le' &
    264           ,' fichier restart ne correspond pas a celle lue dans le run.def'
    265      IF (raz_date == 1) THEN
    266         WRITE(lunout,*)'Je prends l''heure lue dans run.def'
    267         start_time = starttime
    268      ELSE
    269         CALL abort_gcm("gcm", "'Je m''arrete'", 1)
    270      ENDIF
     260    WRITE(lunout, *)' GCM: Attention l''heure de depart lue dans le' &
     261            , ' fichier restart ne correspond pas a celle lue dans le run.def'
     262    IF (raz_date == 1) THEN
     263      WRITE(lunout, *)'Je prends l''heure lue dans run.def'
     264      start_time = starttime
     265    ELSE
     266      CALL abort_gcm("gcm", "'Je m''arrete'", 1)
     267    ENDIF
    271268  ENDIF
    272269  IF (raz_date == 1) THEN
    273      annee_ref = anneeref
    274      day_ref = dayref
    275      day_ini = dayref
    276      itau_dyn = 0
    277      itau_phy = 0
    278      time_0 = 0.
    279      WRITE(lunout,*) &
    280           'GCM: On reinitialise a la date lue dans gcm.def'
     270    annee_ref = anneeref
     271    day_ref = dayref
     272    day_ini = dayref
     273    itau_dyn = 0
     274    itau_phy = 0
     275    time_0 = 0.
     276    WRITE(lunout, *) &
     277            'GCM: On reinitialise a la date lue dans gcm.def'
    281278  ELSE IF (annee_ref /= anneeref .OR. day_ref /= dayref) THEN
    282      WRITE(lunout,*) &
    283           'GCM: Attention les dates initiales lues dans le fichier'
    284      WRITE(lunout,*) &
    285           ' restart ne correspondent pas a celles lues dans '
    286      WRITE(lunout,*)' gcm.def'
    287      WRITE(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
    288      WRITE(lunout,*)' day_ref=',day_ref," dayref=",dayref
    289      WRITE(lunout,*)' Pas de remise a zero'
     279    WRITE(lunout, *) &
     280            'GCM: Attention les dates initiales lues dans le fichier'
     281    WRITE(lunout, *) &
     282            ' restart ne correspondent pas a celles lues dans '
     283    WRITE(lunout, *)' gcm.def'
     284    WRITE(lunout, *)' annee_ref=', annee_ref, " anneeref=", anneeref
     285    WRITE(lunout, *)' day_ref=', day_ref, " dayref=", dayref
     286    WRITE(lunout, *)' Pas de remise a zero'
    290287  ENDIF
    291288
    … …  
    323320  CALL ioconf_startdate(INT(jD_ref), jH_ref)
    324321
    325   WRITE(lunout,*)'DEBUG'
    326   WRITE(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
    327   WRITE(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
    328   CALL ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
    329   WRITE(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
    330   WRITE(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
    331 
     322  WRITE(lunout, *)'DEBUG'
     323  WRITE(lunout, *)'annee_ref, mois, day_ref, heure, jD_ref'
     324  WRITE(lunout, *)annee_ref, mois, day_ref, heure, jD_ref
     325  CALL ju2ymds(jD_ref + jH_ref, an, mois, jour, heure)
     326  WRITE(lunout, *)'jD_ref+jH_ref,an, mois, jour, heure'
     327  WRITE(lunout, *)jD_ref + jH_ref, an, mois, jour, heure
    332328
    333329  IF (iflag_phys==1) THEN
    334      ! these initialisations have already been done (via iniacademic)
    335      ! if running in SW or Newtonian mode
    336      !-----------------------------------------------------------------------
    337      !   Initialisation des constantes dynamiques :
    338      !   ------------------------------------------
    339      dtvr = zdtvr
    340      CALL iniconst
    341 
    342      !-----------------------------------------------------------------------
    343      !   Initialisation de la geometrie :
    344      !   --------------------------------
    345      CALL inigeom
    346 
    347      !-----------------------------------------------------------------------
    348      !   Initialisation du filtre :
    349      !   --------------------------
    350      CALL inifilr
     330    ! these initialisations have already been done (via iniacademic)
     331    ! if running in SW or Newtonian mode
     332    !-----------------------------------------------------------------------
     333    !   Initialisation des constantes dynamiques :
     334    !   ------------------------------------------
     335    dtvr = zdtvr
     336    CALL iniconst
     337
     338    !-----------------------------------------------------------------------
     339    !   Initialisation de la geometrie :
     340    !   --------------------------------
     341    CALL inigeom
     342
     343    !-----------------------------------------------------------------------
     344    !   Initialisation du filtre :
     345    !   --------------------------
     346    CALL inifilr
    351347  ENDIF ! of if (iflag_phys.EQ.1)
    352348
    … …  
    355351  !   ----------------------------------
    356352
    357   CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   , &
    358        tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
     353  CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, &
     354          tetagdiv, tetagrot, tetatemp, vert_prof_dissip)
    359355
    360356  !  numero de stockage pour les fichiers de redemarrage:
    … …  
    364360  !   ------------------------
    365361
    366 
    367362  IF (nday>=0) THEN
    368      day_end = day_ini + nday
     363    day_end = day_ini + nday
    369364  else
    370      day_end = day_ini - nday/day_step
    371   ENDIF
    372   WRITE(lunout,300)day_ini,day_end
    373 300 FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
     365    day_end = day_ini - nday / day_step
     366  ENDIF
     367  WRITE(lunout, 300)day_ini, day_end
     368  300 FORMAT('1'/, 15x, 'run du jour', i7, 2x, 'au jour', i7//)
    374369
    375370  CALL ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
    376   write (lunout,301)jour, mois, an
     371  write (lunout, 301)jour, mois, an
    377372  CALL ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
    378   write (lunout,302)jour, mois, an
    379 301 FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
    380 302 FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
     373  write (lunout, 302)jour, mois, an
     374  301 FORMAT('1'/, 15x, 'run du ', i2, '/', i2, '/', i4)
     375  302 FORMAT('1'/, 15x, '    au ', i2, '/', i2, '/', i4)
    381376
    382377  !-----------------------------------------------------------------------
    … …  
    387382    ! Physics:
    388383    IF (CPPKEY_PHYS) THEN
    389       CALL iniphysiq(iim,jjm,llm, &
    390           (jjm-1)*iim+2,comm_lmdz, &
    391           daysec,day_ini,dtphys/nsplit_phys, &
    392           rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, &
    393           iflag_phys)
     384      CALL iniphysiq(iim, jjm, llm, &
     385              (jjm - 1) * iim + 2, comm_lmdz, &
     386              daysec, day_ini, dtphys / nsplit_phys, &
     387              rlatu, rlatv, rlonu, rlonv, aire, cu, cv, rad, g, r, cpp, &
     388              iflag_phys)
    394389    END IF
    395390  ENDIF ! of IF ((iflag_phys==1).OR.(iflag_phys>=100))
    … …  
    405400  time_step = zdtvr
    406401  IF (ok_dyn_ins) THEN
    407      ! initialize output file for instantaneous outputs
    408      ! t_ops = iecri * daysec ! do operations every t_ops
    409      t_ops =((1.0*iecri)/day_step) * daysec 
    410      t_wrt = daysec ! iecri * daysec ! write output every t_wrt
    411      CALL inithist(day_ref,annee_ref,time_step, &
    412           t_ops,t_wrt)
    413   ENDIF
    414 
    415   IF (ok_dyn_ave) THEN 
    416      ! initialize output file for averaged outputs
    417      t_ops = iperiod * time_step ! do operations every t_ops
    418      t_wrt = periodav * daysec   ! write output every t_wrt
    419      CALL initdynav(day_ref,annee_ref,time_step, &
    420           t_ops,t_wrt)
     402    ! initialize output file for instantaneous outputs
     403    ! t_ops = iecri * daysec ! do operations every t_ops
     404    t_ops = ((1.0 * iecri) / day_step) * daysec
     405    t_wrt = daysec ! iecri * daysec ! write output every t_wrt
     406    CALL inithist(day_ref, annee_ref, time_step, &
     407            t_ops, t_wrt)
     408  ENDIF
     409
     410  IF (ok_dyn_ave) THEN
     411    ! initialize output file for averaged outputs
     412    t_ops = iperiod * time_step ! do operations every t_ops
     413    t_wrt = periodav * daysec   ! write output every t_wrt
     414    CALL initdynav(day_ref, annee_ref, time_step, &
     415            t_ops, t_wrt)
    421416  END IF
    422   dtav = iperiod*dtvr/daysec
     417  dtav = iperiod * dtvr / daysec
    423418
    424419  !  Choix des frequences de stokage pour le offline
    425420  !      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
    426421  !      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
    427   istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
    428   istphy=istdyn/iphysiq     
     422  istdyn = day_step / 4     ! stockage toutes les 6h=1jour/12
     423  istphy = istdyn / iphysiq
    429424
    430425  !-----------------------------------------------------------------------
    … …  
    438433  !       WRITE(78,*) 'q',q
    439434
    440 
    441   CALL leapfrog(ucov,vcov,teta,ps,masse,phis,q,time_0)
     435  CALL leapfrog(ucov, vcov, teta, ps, masse, phis, q, time_0)
    442436
    443437END PROGRAM gcm

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/guide_mod.F90

    r5117 r5118  
    1 
    21! $Id$
    32
    43MODULE guide_mod
    54
    6 !=======================================================================
    7 !   Auteur:  F.Hourdin
    8 !            F. Codron 01/09
    9 !=======================================================================
     5  !=======================================================================
     6  !   Auteur:  F.Hourdin
     7  !            F. Codron 01/09
     8  !=======================================================================
    109
    1110  USE getparam, ONLY: ini_getparam, fin_getparam, getpar
    … …  
    2019  IMPLICIT NONE
    2120
    22 ! ---------------------------------------------
    23 ! Declarations des cles logiques et parametres
    24 ! ---------------------------------------------
    25   INTEGER, PRIVATE, SAVE  :: iguide_read,iguide_int,iguide_sav
    26   INTEGER, PRIVATE, SAVE  :: nlevnc, guide_plevs
    27   LOGICAL, PRIVATE, SAVE  :: guide_u,guide_v,guide_T,guide_Q,guide_P
    28   LOGICAL, PRIVATE, SAVE  :: guide_hr,guide_teta 
    29   LOGICAL, PRIVATE, SAVE  :: guide_BL,guide_reg,guide_add,gamma4,guide_zon
    30   LOGICAL, PRIVATE, SAVE  :: invert_p,invert_y,ini_anal
    31   LOGICAL, PRIVATE, SAVE  :: guide_2D,guide_sav,guide_modele
    32 !FC
    33   LOGICAL, PRIVATE, SAVE  :: convert_Pa
    34  
    35   REAL, PRIVATE, SAVE     :: tau_min_u,tau_max_u
    36   REAL, PRIVATE, SAVE     :: tau_min_v,tau_max_v
    37   REAL, PRIVATE, SAVE     :: tau_min_T,tau_max_T
    38   REAL, PRIVATE, SAVE     :: tau_min_Q,tau_max_Q
    39   REAL, PRIVATE, SAVE     :: tau_min_P,tau_max_P
    40 
    41   REAL, PRIVATE, SAVE     :: lat_min_g,lat_max_g
    42   REAL, PRIVATE, SAVE     :: lon_min_g,lon_max_g
    43   REAL, PRIVATE, SAVE     :: tau_lon,tau_lat
    44 
    45   REAL, PRIVATE, SAVE     :: plim_guide_BL
    46 
    47   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: alpha_u,alpha_v
    48   REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE, SAVE     :: alpha_T,alpha_Q
    49   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: alpha_P,alpha_pcor
    50  
    51 ! ---------------------------------------------
    52 ! Variables de guidage
    53 ! ---------------------------------------------
    54 ! Variables des fichiers de guidage
    55   REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE   :: unat1,unat2
    56   REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE   :: vnat1,vnat2
    57   REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE   :: tnat1,tnat2
    58   REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE   :: qnat1,qnat2
    59   REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE   :: pnat1,pnat2
    60   REAL, ALLOCATABLE, DIMENSION(:,:),   PRIVATE, SAVE   :: psnat1,psnat2
    61   REAL, ALLOCATABLE, DIMENSION(:),     PRIVATE, SAVE   :: apnc,bpnc
    62 ! Variables aux dimensions du modele
    63   REAL, ALLOCATABLE, DIMENSION(:,:),   PRIVATE, SAVE   :: ugui1,ugui2
    64   REAL, ALLOCATABLE, DIMENSION(:,:),   PRIVATE, SAVE   :: vgui1,vgui2
    65   REAL, ALLOCATABLE, DIMENSION(:,:),   PRIVATE, SAVE   :: tgui1,tgui2
    66   REAL, ALLOCATABLE, DIMENSION(:,:),   PRIVATE, SAVE   :: qgui1,qgui2
    67   REAL, ALLOCATABLE, DIMENSION(:),   PRIVATE, SAVE   :: psgui1,psgui2
     21  ! ---------------------------------------------
     22  ! Declarations des cles logiques et parametres
     23  ! ---------------------------------------------
     24  INTEGER, PRIVATE, SAVE :: iguide_read, iguide_int, iguide_sav
     25  INTEGER, PRIVATE, SAVE :: nlevnc, guide_plevs
     26  LOGICAL, PRIVATE, SAVE :: guide_u, guide_v, guide_T, guide_Q, guide_P
     27  LOGICAL, PRIVATE, SAVE :: guide_hr, guide_teta
     28  LOGICAL, PRIVATE, SAVE :: guide_BL, guide_reg, guide_add, gamma4, guide_zon
     29  LOGICAL, PRIVATE, SAVE :: invert_p, invert_y, ini_anal
     30  LOGICAL, PRIVATE, SAVE :: guide_2D, guide_sav, guide_modele
     31  !FC
     32  LOGICAL, PRIVATE, SAVE :: convert_Pa
     33
     34  REAL, PRIVATE, SAVE :: tau_min_u, tau_max_u
     35  REAL, PRIVATE, SAVE :: tau_min_v, tau_max_v
     36  REAL, PRIVATE, SAVE :: tau_min_T, tau_max_T
     37  REAL, PRIVATE, SAVE :: tau_min_Q, tau_max_Q
     38  REAL, PRIVATE, SAVE :: tau_min_P, tau_max_P
     39
     40  REAL, PRIVATE, SAVE :: lat_min_g, lat_max_g
     41  REAL, PRIVATE, SAVE :: lon_min_g, lon_max_g
     42  REAL, PRIVATE, SAVE :: tau_lon, tau_lat
     43
     44  REAL, PRIVATE, SAVE :: plim_guide_BL
     45
     46  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_u, alpha_v
     47  REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE, SAVE :: alpha_T, alpha_Q
     48  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_P, alpha_pcor
     49
     50  ! ---------------------------------------------
     51  ! Variables de guidage
     52  ! ---------------------------------------------
     53  ! Variables des fichiers de guidage
     54  REAL, ALLOCATABLE, DIMENSION(:, :, :), PRIVATE, SAVE :: unat1, unat2
     55  REAL, ALLOCATABLE, DIMENSION(:, :, :), PRIVATE, SAVE :: vnat1, vnat2
     56  REAL, ALLOCATABLE, DIMENSION(:, :, :), PRIVATE, SAVE :: tnat1, tnat2
     57  REAL, ALLOCATABLE, DIMENSION(:, :, :), PRIVATE, SAVE :: qnat1, qnat2
     58  REAL, ALLOCATABLE, DIMENSION(:, :, :), PRIVATE, SAVE :: pnat1, pnat2
     59  REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE, SAVE :: psnat1, psnat2
     60  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: apnc, bpnc
     61  ! Variables aux dimensions du modele
     62  REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE, SAVE :: ugui1, ugui2
     63  REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE, SAVE :: vgui1, vgui2
     64  REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE, SAVE :: tgui1, tgui2
     65  REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE, SAVE :: qgui1, qgui2
     66  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: psgui1, psgui2
    6867
    6968CONTAINS
    70 !=======================================================================
     69  !=======================================================================
    7170
    7271  SUBROUTINE guide_init
    … …  
    7776
    7877    IMPLICIT NONE
    79  
     78
    8079    INCLUDE "dimensions.h"
    8180    INCLUDE "paramet.h"
    8281
    83     INTEGER                :: error,ncidpl,rid,rcod
    84     CHARACTER (len = 80)   :: abort_message
    85     CHARACTER (len = 20)   :: modname = 'guide_init'
    86     CHARACTER (len = 20)   :: namedim
    87 
    88 ! ---------------------------------------------
    89 ! Lecture des parametres: 
    90 ! ---------------------------------------------
     82    INTEGER :: error, ncidpl, rid, rcod
     83    CHARACTER (len = 80) :: abort_message
     84    CHARACTER (len = 20) :: modname = 'guide_init'
     85    CHARACTER (len = 20) :: namedim
     86
     87    ! ---------------------------------------------
     88    ! Lecture des parametres:
     89    ! ---------------------------------------------
    9190    CALL ini_getparam("nudging_parameters_out.txt")
    92 ! Variables guidees
    93     CALL getpar('guide_u',.TRUE.,guide_u,'guidage de u')
    94     CALL getpar('guide_v',.TRUE.,guide_v,'guidage de v')
    95     CALL getpar('guide_T',.TRUE.,guide_T,'guidage de T')
    96     CALL getpar('guide_P',.TRUE.,guide_P,'guidage de P')
    97     CALL getpar('guide_Q',.TRUE.,guide_Q,'guidage de Q')
    98     CALL getpar('guide_hr',.TRUE.,guide_hr,'guidage de Q par H.R')
    99     CALL getpar('guide_teta',.FALSE.,guide_teta,'guidage de T par Teta')
    100 
    101     CALL getpar('guide_add',.FALSE.,guide_add,'forçage constant?')
    102     CALL getpar('guide_zon',.FALSE.,guide_zon,'guidage moy zonale')
     91    ! Variables guidees
     92    CALL getpar('guide_u', .TRUE., guide_u, 'guidage de u')
     93    CALL getpar('guide_v', .TRUE., guide_v, 'guidage de v')
     94    CALL getpar('guide_T', .TRUE., guide_T, 'guidage de T')
     95    CALL getpar('guide_P', .TRUE., guide_P, 'guidage de P')
     96    CALL getpar('guide_Q', .TRUE., guide_Q, 'guidage de Q')
     97    CALL getpar('guide_hr', .TRUE., guide_hr, 'guidage de Q par H.R')
     98    CALL getpar('guide_teta', .FALSE., guide_teta, 'guidage de T par Teta')
     99
     100    CALL getpar('guide_add', .FALSE., guide_add, 'forçage constant?')
     101    CALL getpar('guide_zon', .FALSE., guide_zon, 'guidage moy zonale')
    103102    IF (guide_zon .AND. abs(grossismx - 1.) > 0.01) &
    104          CALL abort_gcm("guide_init", &
    105          "zonal nudging requires grid regular in longitude", 1)
    106 
    107 !   Constantes de rappel. Unite : fraction de jour
    108     CALL getpar('tau_min_u',0.02,tau_min_u,'Cste de rappel min, u')
    109     CALL getpar('tau_max_u', 10.,tau_max_u,'Cste de rappel max, u')
    110     CALL getpar('tau_min_v',0.02,tau_min_v,'Cste de rappel min, v')
    111     CALL getpar('tau_max_v', 10.,tau_max_v,'Cste de rappel max, v')
    112     CALL getpar('tau_min_T',0.02,tau_min_T,'Cste de rappel min, T')
    113     CALL getpar('tau_max_T', 10.,tau_max_T,'Cste de rappel max, T')
    114     CALL getpar('tau_min_Q',0.02,tau_min_Q,'Cste de rappel min, Q')
    115     CALL getpar('tau_max_Q', 10.,tau_max_Q,'Cste de rappel max, Q')
    116     CALL getpar('tau_min_P',0.02,tau_min_P,'Cste de rappel min, P')
    117     CALL getpar('tau_max_P', 10.,tau_max_P,'Cste de rappel max, P')
    118     CALL getpar('gamma4',.FALSE.,gamma4,'Zone sans rappel elargie')
    119     CALL getpar('guide_BL',.TRUE.,guide_BL,'guidage dans C.Lim')
    120     CALL getpar('plim_guide_BL',85000.,plim_guide_BL,'BL top presnivs value')
    121 
    122 
    123 ! Sauvegarde du forçage
    124     CALL getpar('guide_sav',.FALSE.,guide_sav,'sauvegarde guidage')
    125     CALL getpar('iguide_sav',4,iguide_sav,'freq. sauvegarde guidage')
     103            CALL abort_gcm("guide_init", &
     104                    "zonal nudging requires grid regular in longitude", 1)
     105
     106    !   Constantes de rappel. Unite : fraction de jour
     107    CALL getpar('tau_min_u', 0.02, tau_min_u, 'Cste de rappel min, u')
     108    CALL getpar('tau_max_u', 10., tau_max_u, 'Cste de rappel max, u')
     109    CALL getpar('tau_min_v', 0.02, tau_min_v, 'Cste de rappel min, v')
     110    CALL getpar('tau_max_v', 10., tau_max_v, 'Cste de rappel max, v')
     111    CALL getpar('tau_min_T', 0.02, tau_min_T, 'Cste de rappel min, T')
     112    CALL getpar('tau_max_T', 10., tau_max_T, 'Cste de rappel max, T')
     113    CALL getpar('tau_min_Q', 0.02, tau_min_Q, 'Cste de rappel min, Q')
     114    CALL getpar('tau_max_Q', 10., tau_max_Q, 'Cste de rappel max, Q')
     115    CALL getpar('tau_min_P', 0.02, tau_min_P, 'Cste de rappel min, P')
     116    CALL getpar('tau_max_P', 10., tau_max_P, 'Cste de rappel max, P')
     117    CALL getpar('gamma4', .FALSE., gamma4, 'Zone sans rappel elargie')
     118    CALL getpar('guide_BL', .TRUE., guide_BL, 'guidage dans C.Lim')
     119    CALL getpar('plim_guide_BL', 85000., plim_guide_BL, 'BL top presnivs value')
     120
     121
     122    ! Sauvegarde du forçage
     123    CALL getpar('guide_sav', .FALSE., guide_sav, 'sauvegarde guidage')
     124    CALL getpar('iguide_sav', 4, iguide_sav, 'freq. sauvegarde guidage')
    126125    ! frequences f>0: fx/jour; f<0: tous les f jours; f=0: 1 seule fois.
    127126    IF (iguide_sav>0) THEN
    128        iguide_sav=day_step/iguide_sav
     127      iguide_sav = day_step / iguide_sav
    129128    ELSE if (iguide_sav == 0) THEN
    130        iguide_sav = huge(0)
     129      iguide_sav = huge(0)
    131130    ELSE
    132        iguide_sav=day_step*iguide_sav
    133     ENDIF
    134 
    135 ! Guidage regional seulement (sinon constant ou suivant le zoom)
    136     CALL getpar('guide_reg',.FALSE.,guide_reg,'guidage regional')
    137     CALL getpar('lat_min_g',-90.,lat_min_g,'Latitude mini guidage ')
    138     CALL getpar('lat_max_g', 90.,lat_max_g,'Latitude maxi guidage ')
    139     CALL getpar('lon_min_g',-180.,lon_min_g,'longitude mini guidage ')
    140     CALL getpar('lon_max_g', 180.,lon_max_g,'longitude maxi guidage ')
    141     CALL getpar('tau_lat', 5.,tau_lat,'raideur lat guide regional ')
    142     CALL getpar('tau_lon', 5.,tau_lon,'raideur lon guide regional ')
    143 
    144 ! Parametres pour lecture des fichiers
    145     CALL getpar('iguide_read',4,iguide_read,'freq. lecture guidage')
    146     CALL getpar('iguide_int',4,iguide_int,'freq. interpolation vert')
     131      iguide_sav = day_step * iguide_sav
     132    ENDIF
     133
     134    ! Guidage regional seulement (sinon constant ou suivant le zoom)
     135    CALL getpar('guide_reg', .FALSE., guide_reg, 'guidage regional')
     136    CALL getpar('lat_min_g', -90., lat_min_g, 'Latitude mini guidage ')
     137    CALL getpar('lat_max_g', 90., lat_max_g, 'Latitude maxi guidage ')
     138    CALL getpar('lon_min_g', -180., lon_min_g, 'longitude mini guidage ')
     139    CALL getpar('lon_max_g', 180., lon_max_g, 'longitude maxi guidage ')
     140    CALL getpar('tau_lat', 5., tau_lat, 'raideur lat guide regional ')
     141    CALL getpar('tau_lon', 5., tau_lon, 'raideur lon guide regional ')
     142
     143    ! Parametres pour lecture des fichiers
     144    CALL getpar('iguide_read', 4, iguide_read, 'freq. lecture guidage')
     145    CALL getpar('iguide_int', 4, iguide_int, 'freq. interpolation vert')
    147146    IF (iguide_int==0) THEN
    148         iguide_int=1
     147      iguide_int = 1
    149148    ELSEIF (iguide_int>0) THEN
    150         iguide_int=day_step/iguide_int
     149      iguide_int = day_step / iguide_int
    151150    ELSE
    152         iguide_int=day_step*iguide_int
    153     ENDIF
    154     CALL getpar('guide_plevs',0,guide_plevs,'niveaux pression fichiers guidage')
     151      iguide_int = day_step * iguide_int
     152    ENDIF
     153    CALL getpar('guide_plevs', 0, guide_plevs, 'niveaux pression fichiers guidage')
    155154    ! Pour compatibilite avec ancienne version avec guide_modele
    156     CALL getpar('guide_modele',.FALSE.,guide_modele,'niveaux pression ap+bp*psol')
     155    CALL getpar('guide_modele', .FALSE., guide_modele, 'niveaux pression ap+bp*psol')
    157156    IF (guide_modele) THEN
    158         guide_plevs=1
    159     ENDIF
    160 !FC
    161     CALL getpar('convert_Pa',.TRUE.,convert_Pa,'Convert Pressure levels in Pa')
     157      guide_plevs = 1
     158    ENDIF
     159    !FC
     160    CALL getpar('convert_Pa', .TRUE., convert_Pa, 'Convert Pressure levels in Pa')
    162161    ! Fin raccord
    163     CALL getpar('ini_anal',.FALSE.,ini_anal,'Etat initial = analyse')
    164     CALL getpar('guide_invertp',.TRUE.,invert_p,'niveaux p inverses')
    165     CALL getpar('guide_inverty',.TRUE.,invert_y,'inversion N-S')
    166     CALL getpar('guide_2D',.FALSE.,guide_2D,'fichier guidage lat-P')
     162    CALL getpar('ini_anal', .FALSE., ini_anal, 'Etat initial = analyse')
     163    CALL getpar('guide_invertp', .TRUE., invert_p, 'niveaux p inverses')
     164    CALL getpar('guide_inverty', .TRUE., invert_y, 'inversion N-S')
     165    CALL getpar('guide_2D', .FALSE., guide_2D, 'fichier guidage lat-P')
    167166
    168167    CALL fin_getparam
    169    
    170 ! ---------------------------------------------
    171 ! Determination du nombre de niveaux verticaux
    172 ! des fichiers guidage
    173 ! ---------------------------------------------
    174     ncidpl=-99
     168
     169    ! ---------------------------------------------
     170    ! Determination du nombre de niveaux verticaux
     171    ! des fichiers guidage
     172    ! ---------------------------------------------
     173    ncidpl = -99
    175174    IF (guide_plevs==1) THEN
    176        IF (ncidpl==-99) THEN
    177           rcod=nf90_open('apbp.nc',nf90_nowrite, ncidpl)
    178           IF (rcod/=nf90_noerr) THEN
    179              abort_message=' Nudging error -> no file apbp.nc'
    180              CALL abort_gcm(modname,abort_message,1)
    181           endif
    182        endif
     175      IF (ncidpl==-99) THEN
     176        rcod = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
     177        IF (rcod/=nf90_noerr) THEN
     178          abort_message = ' Nudging error -> no file apbp.nc'
     179          CALL abort_gcm(modname, abort_message, 1)
     180        endif
     181      endif
    183182    elseif (guide_plevs==2) THEN
    184        IF (ncidpl==-99) THEN
    185           rcod=nf90_open('P.nc',nf90_nowrite,ncidpl)
    186           IF (rcod/=nf90_noerr) THEN
    187              abort_message=' Nudging error -> no file P.nc'
    188              CALL abort_gcm(modname,abort_message,1)
    189           endif
    190        endif
     183      IF (ncidpl==-99) THEN
     184        rcod = nf90_open('P.nc', nf90_nowrite, ncidpl)
     185        IF (rcod/=nf90_noerr) THEN
     186          abort_message = ' Nudging error -> no file P.nc'
     187          CALL abort_gcm(modname, abort_message, 1)
     188        endif
     189      endif
    191190
    192191    elseif (guide_u) THEN
    193            IF (ncidpl==-99) THEN
    194                rcod=nf90_open('u.nc',nf90_nowrite,ncidpl)
    195                IF (rcod/=nf90_noerr) THEN
    196                   CALL abort_gcm(modname, &
    197                        ' Nudging error -> no file u.nc',1)
    198                endif
    199            endif
     192      IF (ncidpl==-99) THEN
     193        rcod = nf90_open('u.nc', nf90_nowrite, ncidpl)
     194        IF (rcod/=nf90_noerr) THEN
     195          CALL abort_gcm(modname, &
     196                  ' Nudging error -> no file u.nc', 1)
     197        endif
     198      endif
    200199
    201200    elseif (guide_v) THEN
    202            IF (ncidpl==-99) THEN
    203                rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
    204                IF (rcod/=nf90_noerr) THEN
    205                   CALL abort_gcm(modname, &
    206                        ' Nudging error -> no file v.nc',1)
    207                endif
    208            endif
     201      IF (ncidpl==-99) THEN
     202        rcod = nf90_open('v.nc', nf90_nowrite, ncidpl)
     203        IF (rcod/=nf90_noerr) THEN
     204          CALL abort_gcm(modname, &
     205                  ' Nudging error -> no file v.nc', 1)
     206        endif
     207      endif
    209208    elseif (guide_T) THEN
    210            IF (ncidpl==-99) THEN
    211                rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
    212                IF (rcod/=nf90_noerr) THEN
    213                   CALL abort_gcm(modname, &
    214                        ' Nudging error -> no file T.nc',1)
    215                endif
    216            endif
     209      IF (ncidpl==-99) THEN
     210        rcod = nf90_open('T.nc', nf90_nowrite, ncidpl)
     211        IF (rcod/=nf90_noerr) THEN
     212          CALL abort_gcm(modname, &
     213                  ' Nudging error -> no file T.nc', 1)
     214        endif
     215      endif
    217216    elseif (guide_Q) THEN
    218            IF (ncidpl==-99) THEN
    219                rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
    220                IF (rcod/=nf90_noerr) THEN
    221                   CALL abort_gcm(modname, &
    222                        ' Nudging error -> no file hur.nc',1)
    223                endif
    224            endif
    225 
    226 
    227     endif
    228     error=nf90_inq_dimid(ncidpl,'LEVEL',rid)
    229     IF (error/=nf90_noerr) error=nf90_inq_dimid(ncidpl,'PRESSURE',rid)
     217      IF (ncidpl==-99) THEN
     218        rcod = nf90_open('hur.nc', nf90_nowrite, ncidpl)
     219        IF (rcod/=nf90_noerr) THEN
     220          CALL abort_gcm(modname, &
     221                  ' Nudging error -> no file hur.nc', 1)
     222        endif
     223      endif
     224
     225    endif
     226    error = nf90_inq_dimid(ncidpl, 'LEVEL', rid)
     227    IF (error/=nf90_noerr) error = nf90_inq_dimid(ncidpl, 'PRESSURE', rid)
    230228    IF (error/=nf90_noerr) THEN
    231         CALL abort_gcm(modname,'Nudging: error reading pressure levels',1)
    232     ENDIF
    233     error=nf90_inquire_dimension(ncidpl,rid,len=nlevnc)
    234     WRITE(*,*)trim(modname)//' : number of vertical levels nlevnc', nlevnc
     229      CALL abort_gcm(modname, 'Nudging: error reading pressure levels', 1)
     230    ENDIF
     231    error = nf90_inquire_dimension(ncidpl, rid, len = nlevnc)
     232    WRITE(*, *)trim(modname) // ' : number of vertical levels nlevnc', nlevnc
    235233    rcod = nf90_close(ncidpl)
    236234
    237 ! ---------------------------------------------
    238 ! Allocation des variables
    239 ! ---------------------------------------------
    240     abort_message='nudging allocation error'
     235    ! ---------------------------------------------
     236    ! Allocation des variables
     237    ! ---------------------------------------------
     238    abort_message = 'nudging allocation error'
    241239
    242240    ALLOCATE(apnc(nlevnc), stat = error)
    243     IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
     241    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
    244242    ALLOCATE(bpnc(nlevnc), stat = error)
    245     IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    246     apnc=0.;bpnc=0.
     243    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     244    apnc = 0.;bpnc = 0.
    247245
    248246    ALLOCATE(alpha_pcor(llm), stat = error)
    249     IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
     247    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
    250248    ALLOCATE(alpha_u(ip1jmp1), stat = error)
    251     IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
     249    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
    252250    ALLOCATE(alpha_v(ip1jm), stat = error)
    253     IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
     251    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
    254252    ALLOCATE(alpha_T(iip1, jjp1), stat = error)
    255     IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
     253    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
    256254    ALLOCATE(alpha_Q(iip1, jjp1), stat = error)
    257     IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
     255    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
    258256    ALLOCATE(alpha_P(ip1jmp1), stat = error)
    259     IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    260     alpha_u=0.;alpha_v=0;alpha_T=0;alpha_Q=0;alpha_P=0
    261    
     257    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     258    alpha_u = 0.;alpha_v = 0;alpha_T = 0;alpha_Q = 0;alpha_P = 0
     259
    262260    IF (guide_u) THEN
    263         ALLOCATE(unat1(iip1,jjp1,nlevnc), stat = error)
    264         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    265         ALLOCATE(ugui1(ip1jmp1,llm), stat = error)
    266         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    267         ALLOCATE(unat2(iip1,jjp1,nlevnc), stat = error)
    268         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    269         ALLOCATE(ugui2(ip1jmp1,llm), stat = error)
    270         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    271         unat1=0.;unat2=0.;ugui1=0.;ugui2=0.
     261      ALLOCATE(unat1(iip1, jjp1, nlevnc), stat = error)
     262      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     263      ALLOCATE(ugui1(ip1jmp1, llm), stat = error)
     264      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     265      ALLOCATE(unat2(iip1, jjp1, nlevnc), stat = error)
     266      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     267      ALLOCATE(ugui2(ip1jmp1, llm), stat = error)
     268      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     269      unat1 = 0.;unat2 = 0.;ugui1 = 0.;ugui2 = 0.
    272270    ENDIF
    273271
    274272    IF (guide_T) THEN
    275         ALLOCATE(tnat1(iip1,jjp1,nlevnc), stat = error)
    276         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    277         ALLOCATE(tgui1(ip1jmp1,llm), stat = error)
    278         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    279         ALLOCATE(tnat2(iip1,jjp1,nlevnc), stat = error)
    280         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    281         ALLOCATE(tgui2(ip1jmp1,llm), stat = error)
    282         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    283         tnat1=0.;tnat2=0.;tgui1=0.;tgui2=0.
    284     ENDIF
    285      
     273      ALLOCATE(tnat1(iip1, jjp1, nlevnc), stat = error)
     274      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     275      ALLOCATE(tgui1(ip1jmp1, llm), stat = error)
     276      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     277      ALLOCATE(tnat2(iip1, jjp1, nlevnc), stat = error)
     278      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     279      ALLOCATE(tgui2(ip1jmp1, llm), stat = error)
     280      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     281      tnat1 = 0.;tnat2 = 0.;tgui1 = 0.;tgui2 = 0.
     282    ENDIF
     283
    286284    IF (guide_Q) THEN
    287         ALLOCATE(qnat1(iip1,jjp1,nlevnc), stat = error)
    288         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    289         ALLOCATE(qgui1(ip1jmp1,llm), stat = error)
    290         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    291         ALLOCATE(qnat2(iip1,jjp1,nlevnc), stat = error)
    292         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    293         ALLOCATE(qgui2(ip1jmp1,llm), stat = error)
    294         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    295         qnat1=0.;qnat2=0.;qgui1=0.;qgui2=0.
     285      ALLOCATE(qnat1(iip1, jjp1, nlevnc), stat = error)
     286      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     287      ALLOCATE(qgui1(ip1jmp1, llm), stat = error)
     288      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     289      ALLOCATE(qnat2(iip1, jjp1, nlevnc), stat = error)
     290      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     291      ALLOCATE(qgui2(ip1jmp1, llm), stat = error)
     292      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     293      qnat1 = 0.;qnat2 = 0.;qgui1 = 0.;qgui2 = 0.
    296294    ENDIF
    297295
    298296    IF (guide_v) THEN
    299         ALLOCATE(vnat1(iip1,jjm,nlevnc), stat = error)
    300         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    301         ALLOCATE(vgui1(ip1jm,llm), stat = error)
    302         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    303         ALLOCATE(vnat2(iip1,jjm,nlevnc), stat = error)
    304         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    305         ALLOCATE(vgui2(ip1jm,llm), stat = error)
    306         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    307         vnat1=0.;vnat2=0.;vgui1=0.;vgui2=0.
     297      ALLOCATE(vnat1(iip1, jjm, nlevnc), stat = error)
     298      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     299      ALLOCATE(vgui1(ip1jm, llm), stat = error)
     300      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     301      ALLOCATE(vnat2(iip1, jjm, nlevnc), stat = error)
     302      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     303      ALLOCATE(vgui2(ip1jm, llm), stat = error)
     304      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     305      vnat1 = 0.;vnat2 = 0.;vgui1 = 0.;vgui2 = 0.
    308306    ENDIF
    309307
    310308    IF (guide_plevs==2) THEN
    311         ALLOCATE(pnat1(iip1,jjp1,nlevnc), stat = error)
    312         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    313         ALLOCATE(pnat2(iip1,jjp1,nlevnc), stat = error)
    314         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    315         pnat1=0.;pnat2=0.;
     309      ALLOCATE(pnat1(iip1, jjp1, nlevnc), stat = error)
     310      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     311      ALLOCATE(pnat2(iip1, jjp1, nlevnc), stat = error)
     312      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     313      pnat1 = 0.;pnat2 = 0.;
    316314    ENDIF
    317315
    318316    IF (guide_P.OR.guide_plevs==1) THEN
    319         ALLOCATE(psnat1(iip1,jjp1), stat = error)
    320         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    321         ALLOCATE(psnat2(iip1,jjp1), stat = error)
    322         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    323         psnat1=0.;psnat2=0.;
     317      ALLOCATE(psnat1(iip1, jjp1), stat = error)
     318      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     319      ALLOCATE(psnat2(iip1, jjp1), stat = error)
     320      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     321      psnat1 = 0.;psnat2 = 0.;
    324322    ENDIF
    325323    IF (guide_P) THEN
    326         ALLOCATE(psgui2(ip1jmp1), stat = error)
    327         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    328         ALLOCATE(psgui1(ip1jmp1), stat = error)
    329         IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
    330         psgui1=0.;psgui2=0.
    331     ENDIF
    332 
    333 ! ---------------------------------------------
    334 !   Lecture du premier etat de guidage.
    335 ! ---------------------------------------------
     324      ALLOCATE(psgui2(ip1jmp1), stat = error)
     325      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     326      ALLOCATE(psgui1(ip1jmp1), stat = error)
     327      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
     328      psgui1 = 0.;psgui2 = 0.
     329    ENDIF
     330
     331    ! ---------------------------------------------
     332    !   Lecture du premier etat de guidage.
     333    ! ---------------------------------------------
    336334    IF (guide_2D) THEN
    337         CALL guide_read2D(1)
     335      CALL guide_read2D(1)
    338336    ELSE
    339         CALL guide_read(1)
    340     ENDIF
    341     IF (guide_v) vnat1=vnat2
    342     IF (guide_u) unat1=unat2
    343     IF (guide_T) tnat1=tnat2
    344     IF (guide_Q) qnat1=qnat2
    345     IF (guide_plevs==2) pnat1=pnat2
    346     IF (guide_P.OR.guide_plevs==1) psnat1=psnat2
     337      CALL guide_read(1)
     338    ENDIF
     339    IF (guide_v) vnat1 = vnat2
     340    IF (guide_u) unat1 = unat2
     341    IF (guide_T) tnat1 = tnat2
     342    IF (guide_Q) qnat1 = qnat2
     343    IF (guide_plevs==2) pnat1 = pnat2
     344    IF (guide_P.OR.guide_plevs==1) psnat1 = psnat2
    347345
    348346  END SUBROUTINE guide_init
    349347
    350 !=======================================================================
    351   SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps)
     348  !=======================================================================
     349  SUBROUTINE guide_main(itau, ucov, vcov, teta, q, masse, ps)
    352350
    353351    USE exner_hyb_m, ONLY: exner_hyb
    354352    USE exner_milieu_m, ONLY: exner_milieu
    355353    USE control_mod, ONLY: day_step, iperiod
    356     USE comconst_mod, ONLY: cpp, dtvr, daysec,kappa
     354    USE comconst_mod, ONLY: cpp, dtvr, daysec, kappa
    357355    USE comvert_mod, ONLY: ap, bp, preff, presnivs, pressure_exner
    358  
     356    USE lmdz_iniprint, ONLY: lunout, prt_level
     357
    359358    IMPLICIT NONE
    360  
     359
    361360    INCLUDE "dimensions.h"
    362361    INCLUDE "paramet.h"
    363     INCLUDE "iniprint.h"
    364362
    365363
    366364    ! Variables entree
    367     INTEGER,                       INTENT(IN)    :: itau !pas de temps
    368     REAL, DIMENSION (ip1jmp1,llm), INTENT(INOUT) :: ucov,teta,q,masse
    369     REAL, DIMENSION (ip1jm,llm),  INTENT(INOUT) :: vcov
    370     REAL, DIMENSION (ip1jmp1),     INTENT(INOUT) :: ps
     365    INTEGER, INTENT(IN) :: itau !pas de temps
     366    REAL, DIMENSION (ip1jmp1, llm), INTENT(INOUT) :: ucov, teta, q, masse
     367    REAL, DIMENSION (ip1jm, llm), INTENT(INOUT) :: vcov
     368    REAL, DIMENSION (ip1jmp1), INTENT(INOUT) :: ps
    371369
    372370    ! Variables locales
    373     LOGICAL, SAVE :: first=.TRUE.
    374     LOGICAL       :: f_out ! sortie guidage
    375     REAL, DIMENSION (ip1jmp1,llm) :: f_add ! var aux: champ de guidage
    376     REAL :: pk(ip1jmp1,llm) ! Exner at mid-layers
     371    LOGICAL, SAVE :: first = .TRUE.
     372    LOGICAL :: f_out ! sortie guidage
     373    REAL, DIMENSION (ip1jmp1, llm) :: f_add ! var aux: champ de guidage
     374    REAL :: pk(ip1jmp1, llm) ! Exner at mid-layers
    377375    REAL :: pks(ip1jmp1) ! Exner at the surface
    378376    REAL :: unskap ! 1./kappa
    379     REAL, DIMENSION (ip1jmp1,llmp1) :: p ! Pressure at inter-layers
     377    REAL, DIMENSION (ip1jmp1, llmp1) :: p ! Pressure at inter-layers
    380378    ! Compteurs temps:
    381     INTEGER, SAVE :: step_rea,count_no_rea,itau_test ! lecture guidage
    382     REAL          :: ditau, dday_step
    383     REAL          :: tau,reste ! position entre 2 etats de guidage
    384     REAL, SAVE    :: factt ! pas de temps en fraction de jour
    385    
    386     INTEGER       :: l
    387     CHARACTER(LEN=20) :: modname="guide_main"
    388     CHARACTER (len = 80)   :: abort_message
    389 
    390 
    391 !-----------------------------------------------------------------------
    392 ! Initialisations au premier passage
    393 !-----------------------------------------------------------------------
     379    INTEGER, SAVE :: step_rea, count_no_rea, itau_test ! lecture guidage
     380    REAL :: ditau, dday_step
     381    REAL :: tau, reste ! position entre 2 etats de guidage
     382    REAL, SAVE :: factt ! pas de temps en fraction de jour
     383
     384    INTEGER :: l
     385    CHARACTER(LEN = 20) :: modname = "guide_main"
     386    CHARACTER (len = 80) :: abort_message
     387
     388
     389    !-----------------------------------------------------------------------
     390    ! Initialisations au premier passage
     391    !-----------------------------------------------------------------------
    394392    IF (first) THEN
    395         first=.FALSE.
    396         CALL guide_init
    397         itau_test=1001
    398         step_rea=1
    399         count_no_rea=0
    400 ! Calcul des constantes de rappel
    401         factt=dtvr*iperiod/daysec
    402         CALL tau2alpha(3,iip1,jjm ,factt,tau_min_v,tau_max_v,alpha_v)
    403         CALL tau2alpha(2,iip1,jjp1,factt,tau_min_u,tau_max_u,alpha_u)
    404         CALL tau2alpha(1,iip1,jjp1,factt,tau_min_T,tau_max_T,alpha_T)
    405         CALL tau2alpha(1,iip1,jjp1,factt,tau_min_P,tau_max_P,alpha_P)
    406         CALL tau2alpha(1,iip1,jjp1,factt,tau_min_Q,tau_max_Q,alpha_Q)
    407 ! correction de rappel dans couche limite
    408         IF (guide_BL) THEN
    409              alpha_pcor(:)=1.
    410         else
    411             do l=1,llm
    412                alpha_pcor(l)=(1.+tanh(((plim_guide_BL-presnivs(l))/preff)/0.05))/2.
    413             enddo
    414         endif
    415 ! ini_anal: etat initial egal au guidage       
    416         IF (ini_anal) THEN
    417             CALL guide_interp(ps,teta)
    418             IF (guide_u) ucov=ugui2
    419             IF (guide_v) vcov=ugui2
    420             IF (guide_T) teta=tgui2
    421             IF (guide_Q) q=qgui2
    422             IF (guide_P) THEN
    423                 ps=psgui2
    424                 CALL pression(ip1jmp1,ap,bp,ps,p)
    425                 CALL massdair(p,masse)
    426             ENDIF
    427 
    428         ENDIF
    429 ! Verification structure guidage
    430 !        IF (guide_u) THEN
    431 !            CALL writefield('unat',unat1)
    432 !            CALL writefield('ucov',RESHAPE(ucov,(/iip1,jjp1,llm/)))
    433 !        ENDIF
    434 !        IF (guide_T) THEN
    435 !            CALL writefield('tnat',tnat1)
    436 !            CALL writefield('teta',RESHAPE(teta,(/iip1,jjp1,llm/)))
    437 !        ENDIF
     393      first = .FALSE.
     394      CALL guide_init
     395      itau_test = 1001
     396      step_rea = 1
     397      count_no_rea = 0
     398      ! Calcul des constantes de rappel
     399      factt = dtvr * iperiod / daysec
     400      CALL tau2alpha(3, iip1, jjm, factt, tau_min_v, tau_max_v, alpha_v)
     401      CALL tau2alpha(2, iip1, jjp1, factt, tau_min_u, tau_max_u, alpha_u)
     402      CALL tau2alpha(1, iip1, jjp1, factt, tau_min_T, tau_max_T, alpha_T)
     403      CALL tau2alpha(1, iip1, jjp1, factt, tau_min_P, tau_max_P, alpha_P)
     404      CALL tau2alpha(1, iip1, jjp1, factt, tau_min_Q, tau_max_Q, alpha_Q)
     405      ! correction de rappel dans couche limite
     406      IF (guide_BL) THEN
     407        alpha_pcor(:) = 1.
     408      else
     409        do l = 1, llm
     410          alpha_pcor(l) = (1. + tanh(((plim_guide_BL - presnivs(l)) / preff) / 0.05)) / 2.
     411        enddo
     412      endif
     413      ! ini_anal: etat initial egal au guidage
     414      IF (ini_anal) THEN
     415        CALL guide_interp(ps, teta)
     416        IF (guide_u) ucov = ugui2
     417        IF (guide_v) vcov = ugui2
     418        IF (guide_T) teta = tgui2
     419        IF (guide_Q) q = qgui2
     420        IF (guide_P) THEN
     421          ps = psgui2
     422          CALL pression(ip1jmp1, ap, bp, ps, p)
     423          CALL massdair(p, masse)
     424        ENDIF
     425
     426      ENDIF
     427      ! Verification structure guidage
     428      !        IF (guide_u) THEN
     429      !            CALL writefield('unat',unat1)
     430      !            CALL writefield('ucov',RESHAPE(ucov,(/iip1,jjp1,llm/)))
     431      !        ENDIF
     432      !        IF (guide_T) THEN
     433      !            CALL writefield('tnat',tnat1)
     434      !            CALL writefield('teta',RESHAPE(teta,(/iip1,jjp1,llm/)))
     435      !        ENDIF
    438436
    439437    ENDIF !first
    440438
    441 !-----------------------------------------------------------------------
    442 ! Lecture des fichiers de guidage ?
    443 !-----------------------------------------------------------------------
     439    !-----------------------------------------------------------------------
     440    ! Lecture des fichiers de guidage ?
     441    !-----------------------------------------------------------------------
    444442    IF (iguide_read/=0) THEN
    445       ditau=real(itau)
    446       dday_step=real(day_step)
     443      ditau = real(itau)
     444      dday_step = real(day_step)
    447445      IF (iguide_read<0) THEN
    448           tau=ditau/dday_step/REAL(iguide_read)
     446        tau = ditau / dday_step / REAL(iguide_read)
    449447      ELSE
    450           tau=REAL(iguide_read)*ditau/dday_step
    451       ENDIF
    452       reste=tau-AINT(tau)
     448        tau = REAL(iguide_read) * ditau / dday_step
     449      ENDIF
     450      reste = tau - AINT(tau)
    453451      IF (reste==0.) THEN
    454           IF (itau_test==itau) THEN
    455             WRITE(lunout,*)trim(modname)//' second pass in advreel at itau=',&
    456             itau
    457               abort_message='stopped'
    458               CALL abort_gcm(modname,abort_message,1)
     452        IF (itau_test==itau) THEN
     453          WRITE(lunout, *)trim(modname) // ' second pass in advreel at itau=', &
     454                  itau
     455          abort_message = 'stopped'
     456          CALL abort_gcm(modname, abort_message, 1)
     457        ELSE
     458          IF (guide_v) vnat1 = vnat2
     459          IF (guide_u) unat1 = unat2
     460          IF (guide_T) tnat1 = tnat2
     461          IF (guide_Q) qnat1 = qnat2
     462          IF (guide_plevs==2) pnat1 = pnat2
     463          IF (guide_P.OR.guide_plevs==1) psnat1 = psnat2
     464          step_rea = step_rea + 1
     465          itau_test = itau
     466          WRITE(*, *)trim(modname) // ' Reading nudging files, step ', &
     467                  step_rea, 'after ', count_no_rea, ' skips'
     468          IF (guide_2D) THEN
     469            CALL guide_read2D(step_rea)
    459470          ELSE
    460               IF (guide_v) vnat1=vnat2
    461               IF (guide_u) unat1=unat2
    462               IF (guide_T) tnat1=tnat2
    463               IF (guide_Q) qnat1=qnat2
    464               IF (guide_plevs==2) pnat1=pnat2
    465               IF (guide_P.OR.guide_plevs==1) psnat1=psnat2
    466               step_rea=step_rea+1
    467               itau_test=itau
    468               WRITE(*,*)trim(modname)//' Reading nudging files, step ',&
    469                      step_rea,'after ',count_no_rea,' skips'
    470               IF (guide_2D) THEN
    471                   CALL guide_read2D(step_rea)
    472               ELSE
    473                   CALL guide_read(step_rea)
    474               ENDIF
    475               count_no_rea=0
     471            CALL guide_read(step_rea)
    476472          ENDIF
     473          count_no_rea = 0
     474        ENDIF
    477475      ELSE
    478         count_no_rea=count_no_rea+1
     476        count_no_rea = count_no_rea + 1
    479477
    480478      ENDIF
    481479    ENDIF !iguide_read=0
    482480
    483 !-----------------------------------------------------------------------
    484 ! Interpolation et conversion des champs de guidage
    485 !-----------------------------------------------------------------------
    486     IF (MOD(itau,iguide_int)==0) THEN
    487         CALL guide_interp(ps,teta)
    488     ENDIF
    489 ! Repartition entre 2 etats de guidage
     481    !-----------------------------------------------------------------------
     482    ! Interpolation et conversion des champs de guidage
     483    !-----------------------------------------------------------------------
     484    IF (MOD(itau, iguide_int)==0) THEN
     485      CALL guide_interp(ps, teta)
     486    ENDIF
     487    ! Repartition entre 2 etats de guidage
    490488    IF (iguide_read/=0) THEN
    491         tau=reste
     489      tau = reste
    492490    ELSE
    493         tau=1.
    494     ENDIF
    495 
    496 !-----------------------------------------------------------------------
    497 !   Ajout des champs de guidage
    498 !-----------------------------------------------------------------------
    499 ! Sauvegarde du guidage?
    500     f_out=((MOD(itau,iguide_sav)==0).AND.guide_sav)
     491      tau = 1.
     492    ENDIF
     493
     494    !-----------------------------------------------------------------------
     495    !   Ajout des champs de guidage
     496    !-----------------------------------------------------------------------
     497    ! Sauvegarde du guidage?
     498    f_out = ((MOD(itau, iguide_sav)==0).AND.guide_sav)
    501499    IF (f_out) THEN
    502500      ! compute pressures at layer interfaces
    503       CALL pression(ip1jmp1,ap,bp,ps,p)
     501      CALL pression(ip1jmp1, ap, bp, ps, p)
    504502      IF (pressure_exner) THEN
    505         CALL exner_hyb(ip1jmp1,ps,p,pks,pk)
     503        CALL exner_hyb(ip1jmp1, ps, p, pks, pk)
    506504      else
    507         CALL exner_milieu(ip1jmp1,ps,p,pks,pk)
    508       endif
    509       unskap=1./kappa
     505        CALL exner_milieu(ip1jmp1, ps, p, pks, pk)
     506      endif
     507      unskap = 1. / kappa
    510508      ! Now compute pressures at mid-layer
    511       do l=1,llm
    512         p(:,l)=preff*(pk(:,l)/cpp)**unskap
     509      do l = 1, llm
     510        p(:, l) = preff * (pk(:, l) / cpp)**unskap
    513511      enddo
    514       CALL guide_out("SP",jjp1,llm,p(:,1:llm))
    515     ENDIF
    516    
     512      CALL guide_out("SP", jjp1, llm, p(:, 1:llm))
     513    ENDIF
     514
    517515    IF (guide_u) THEN
    518         IF (guide_add) THEN
    519            f_add=(1.-tau)*ugui1+tau*ugui2
    520         else
    521            f_add=(1.-tau)*ugui1+tau*ugui2-ucov
    522         endif
    523         IF (guide_zon) CALL guide_zonave(1,jjp1,llm,f_add)
    524         CALL guide_addfield(ip1jmp1,llm,f_add,alpha_u)
    525         IF (f_out) CALL guide_out("ua",jjp1,llm,(1.-tau)*ugui1+tau*ugui2)
    526         IF (f_out) CALL guide_out("u",jjp1,llm,ucov)
    527         IF (f_out) CALL guide_out("ucov",jjp1,llm,f_add/factt)
    528         ucov=ucov+f_add
     516      IF (guide_add) THEN
     517        f_add = (1. - tau) * ugui1 + tau * ugui2
     518      else
     519        f_add = (1. - tau) * ugui1 + tau * ugui2 - ucov
     520      endif
     521      IF (guide_zon) CALL guide_zonave(1, jjp1, llm, f_add)
     522      CALL guide_addfield(ip1jmp1, llm, f_add, alpha_u)
     523      IF (f_out) CALL guide_out("ua", jjp1, llm, (1. - tau) * ugui1 + tau * ugui2)
     524      IF (f_out) CALL guide_out("u", jjp1, llm, ucov)
     525      IF (f_out) CALL guide_out("ucov", jjp1, llm, f_add / factt)
     526      ucov = ucov + f_add
    529527    endif
    530528
    531529    IF (guide_T) THEN
    532         IF (guide_add) THEN
    533            f_add=(1.-tau)*tgui1+tau*tgui2
    534         else
    535            f_add=(1.-tau)*tgui1+tau*tgui2-teta
    536         endif
    537         IF (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
    538         CALL guide_addfield(ip1jmp1,llm,f_add,alpha_T)
    539         IF (f_out) CALL guide_out("teta",jjp1,llm,f_add/factt)
    540         teta=teta+f_add
     530      IF (guide_add) THEN
     531        f_add = (1. - tau) * tgui1 + tau * tgui2
     532      else
     533        f_add = (1. - tau) * tgui1 + tau * tgui2 - teta
     534      endif
     535      IF (guide_zon) CALL guide_zonave(2, jjp1, llm, f_add)
     536      CALL guide_addfield(ip1jmp1, llm, f_add, alpha_T)
     537      IF (f_out) CALL guide_out("teta", jjp1, llm, f_add / factt)
     538      teta = teta + f_add
    541539    endif
    542540
    543541    IF (guide_P) THEN
    544         IF (guide_add) THEN
    545            f_add(1:ip1jmp1,1)=(1.-tau)*psgui1+tau*psgui2
    546         else
    547            f_add(1:ip1jmp1,1)=(1.-tau)*psgui1+tau*psgui2-ps
    548         endif
    549         IF (guide_zon) CALL guide_zonave(2,jjp1,1,f_add(1:ip1jmp1,1))
    550         CALL guide_addfield(ip1jmp1,1,f_add(1:ip1jmp1,1),alpha_P)
    551 !        IF (f_out) CALL guide_out("ps",jjp1,1,f_add(1:ip1jmp1,1)/factt)
    552         ps=ps+f_add(1:ip1jmp1,1)
    553         CALL pression(ip1jmp1,ap,bp,ps,p)
    554         CALL massdair(p,masse)
     542      IF (guide_add) THEN
     543        f_add(1:ip1jmp1, 1) = (1. - tau) * psgui1 + tau * psgui2
     544      else
     545        f_add(1:ip1jmp1, 1) = (1. - tau) * psgui1 + tau * psgui2 - ps
     546      endif
     547      IF (guide_zon) CALL guide_zonave(2, jjp1, 1, f_add(1:ip1jmp1, 1))
     548      CALL guide_addfield(ip1jmp1, 1, f_add(1:ip1jmp1, 1), alpha_P)
     549      !        IF (f_out) CALL guide_out("ps",jjp1,1,f_add(1:ip1jmp1,1)/factt)
     550      ps = ps + f_add(1:ip1jmp1, 1)
     551      CALL pression(ip1jmp1, ap, bp, ps, p)
     552      CALL massdair(p, masse)
    555553    endif
    556554
    557555    IF (guide_Q) THEN
    558         IF (guide_add) THEN
    559            f_add=(1.-tau)*qgui1+tau*qgui2
    560         else
    561            f_add=(1.-tau)*qgui1+tau*qgui2-q
    562         endif
    563         IF (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
    564         CALL guide_addfield(ip1jmp1,llm,f_add,alpha_Q)
    565         IF (f_out) CALL guide_out("q",jjp1,llm,f_add/factt)
    566         q=q+f_add
     556      IF (guide_add) THEN
     557        f_add = (1. - tau) * qgui1 + tau * qgui2
     558      else
     559        f_add = (1. - tau) * qgui1 + tau * qgui2 - q
     560      endif
     561      IF (guide_zon) CALL guide_zonave(2, jjp1, llm, f_add)
     562      CALL guide_addfield(ip1jmp1, llm, f_add, alpha_Q)
     563      IF (f_out) CALL guide_out("q", jjp1, llm, f_add / factt)
     564      q = q + f_add
    567565    endif
    568566
    569567    IF (guide_v) THEN
    570         IF (guide_add) THEN
    571            f_add(1:ip1jm,:)=(1.-tau)*vgui1+tau*vgui2
    572         else
    573            f_add(1:ip1jm,:)=(1.-tau)*vgui1+tau*vgui2-vcov
    574         endif
    575         IF (guide_zon) CALL guide_zonave(2,jjm,llm,f_add(1:ip1jm,:))
    576         CALL guide_addfield(ip1jm,llm,f_add(1:ip1jm,:),alpha_v)
    577         IF (f_out) CALL guide_out("v",jjm,llm,vcov)
    578         IF (f_out) CALL guide_out("va",jjm,llm,(1.-tau)*vgui1+tau*vgui2)
    579         IF (f_out) CALL guide_out("vcov",jjm,llm,f_add(1:ip1jm,:)/factt)
    580         vcov=vcov+f_add(1:ip1jm,:)
     568      IF (guide_add) THEN
     569        f_add(1:ip1jm, :) = (1. - tau) * vgui1 + tau * vgui2
     570      else
     571        f_add(1:ip1jm, :) = (1. - tau) * vgui1 + tau * vgui2 - vcov
     572      endif
     573      IF (guide_zon) CALL guide_zonave(2, jjm, llm, f_add(1:ip1jm, :))
     574      CALL guide_addfield(ip1jm, llm, f_add(1:ip1jm, :), alpha_v)
     575      IF (f_out) CALL guide_out("v", jjm, llm, vcov)
     576      IF (f_out) CALL guide_out("va", jjm, llm, (1. - tau) * vgui1 + tau * vgui2)
     577      IF (f_out) CALL guide_out("vcov", jjm, llm, f_add(1:ip1jm, :) / factt)
     578      vcov = vcov + f_add(1:ip1jm, :)
    581579    endif
    582580  END SUBROUTINE guide_main
    583581
    584 !=======================================================================
    585   SUBROUTINE guide_addfield(hsize,vsize,field,alpha)
    586 ! field1=a*field1+alpha*field2
     582  !=======================================================================
     583  SUBROUTINE guide_addfield(hsize, vsize, field, alpha)
     584    ! field1=a*field1+alpha*field2
    587585
    588586    IMPLICIT NONE
    589587
    590588    ! input variables
    591     INTEGER,                      INTENT(IN)    :: hsize
    592     INTEGER,                      INTENT(IN)    :: vsize
    593     REAL, DIMENSION(hsize),       INTENT(IN)    :: alpha
    594     REAL, DIMENSION(hsize,vsize), INTENT(INOUT) :: field
     589    INTEGER, INTENT(IN) :: hsize
     590    INTEGER, INTENT(IN) :: vsize
     591    REAL, DIMENSION(hsize), INTENT(IN) :: alpha
     592    REAL, DIMENSION(hsize, vsize), INTENT(INOUT) :: field
    595593
    596594    ! Local variables
    597595    INTEGER :: l
    598596
    599     do l=1,vsize
    600         field(:,l)=alpha*field(:,l)*alpha_pcor(l)
     597    do l = 1, vsize
     598      field(:, l) = alpha * field(:, l) * alpha_pcor(l)
    601599    enddo
    602600
    603601  END SUBROUTINE guide_addfield
    604602
    605 !=======================================================================
    606   SUBROUTINE guide_zonave(typ,hsize,vsize,field)
     603  !=======================================================================
     604  SUBROUTINE guide_zonave(typ, hsize, vsize, field)
    607605
    608606    USE comconst_mod, ONLY: pi
    609    
     607
    610608    IMPLICIT NONE
    611609
    … …  
    613611    INCLUDE "paramet.h"
    614612    INCLUDE "comgeom.h"
    615    
     613
    616614    ! input/output variables
    617     INTEGER,                           INTENT(IN)    :: typ
    618     INTEGER,                           INTENT(IN)    :: vsize
    619     INTEGER,                           INTENT(IN)    :: hsize
    620     REAL, DIMENSION(hsize*iip1,vsize), INTENT(INOUT) :: field
     615    INTEGER, INTENT(IN) :: typ
     616    INTEGER, INTENT(IN) :: vsize
     617    INTEGER, INTENT(IN) :: hsize
     618    REAL, DIMENSION(hsize * iip1, vsize), INTENT(INOUT) :: field
    621619
    622620    ! Local variables
    623     LOGICAL, SAVE                :: first=.TRUE.
     621    LOGICAL, SAVE :: first = .TRUE.
    624622    INTEGER, DIMENSION (2), SAVE :: imin, imax ! averaging domain
    625     INTEGER                      :: i,j,l,ij
    626     REAL, DIMENSION (iip1)       :: lond       ! longitude in Deg.
    627     REAL, DIMENSION (hsize,vsize):: fieldm     ! zon-averaged field
     623    INTEGER :: i, j, l, ij
     624    REAL, DIMENSION (iip1) :: lond       ! longitude in Deg.
     625    REAL, DIMENSION (hsize, vsize) :: fieldm     ! zon-averaged field
    628626
    629627    IF (first) THEN
    630         first=.FALSE.
    631 !Compute domain for averaging
    632         lond=rlonu*180./pi
    633         imin(1)=1;imax(1)=iip1;
    634         imin(2)=1;imax(2)=iip1;
    635         IF (guide_reg) THEN
    636             DO i=1,iim
    637                 IF (lond(i)<lon_min_g) imin(1)=i
    638                 IF (lond(i)<=lon_max_g) imax(1)=i
    639             ENDDO
    640             lond=rlonv*180./pi
    641             DO i=1,iim
    642                 IF (lond(i)<lon_min_g) imin(2)=i
    643                 IF (lond(i)<=lon_max_g) imax(2)=i
    644             ENDDO
    645         ENDIF
    646     ENDIF
    647 
    648     fieldm=0.
    649     DO l=1,vsize
    650     ! Compute zonal average
    651         DO j=1,hsize
    652             DO i=imin(typ),imax(typ)
    653                 ij=(j-1)*iip1+i
    654                 fieldm(j,l)=fieldm(j,l)+field(ij,l)
    655             ENDDO
    656         ENDDO
    657         fieldm(:,l)=fieldm(:,l)/REAL(imax(typ)-imin(typ)+1)
    658     ! Compute forcing
    659         DO j=1,hsize
    660             DO i=1,iip1
    661                 ij=(j-1)*iip1+i
    662                 field(ij,l)=fieldm(j,l)
    663             ENDDO
     628      first = .FALSE.
     629      !Compute domain for averaging
     630      lond = rlonu * 180. / pi
     631      imin(1) = 1;imax(1) = iip1;
     632      imin(2) = 1;imax(2) = iip1;
     633      IF (guide_reg) THEN
     634        DO i = 1, iim
     635          IF (lond(i)<lon_min_g) imin(1) = i
     636          IF (lond(i)<=lon_max_g) imax(1) = i
    664637        ENDDO
     638        lond = rlonv * 180. / pi
     639        DO i = 1, iim
     640          IF (lond(i)<lon_min_g) imin(2) = i
     641          IF (lond(i)<=lon_max_g) imax(2) = i
     642        ENDDO
     643      ENDIF
     644    ENDIF
     645
     646    fieldm = 0.
     647    DO l = 1, vsize
     648      ! Compute zonal average
     649      DO j = 1, hsize
     650        DO i = imin(typ), imax(typ)
     651          ij = (j - 1) * iip1 + i
     652          fieldm(j, l) = fieldm(j, l) + field(ij, l)
     653        ENDDO
     654      ENDDO
     655      fieldm(:, l) = fieldm(:, l) / REAL(imax(typ) - imin(typ) + 1)
     656      ! Compute forcing
     657      DO j = 1, hsize
     658        DO i = 1, iip1
     659          ij = (j - 1) * iip1 + i
     660          field(ij, l) = fieldm(j, l)
     661        ENDDO
     662      ENDDO
    665663    ENDDO
    666664
    667665  END SUBROUTINE guide_zonave
    668666
    669 !=======================================================================
    670   SUBROUTINE guide_interp(psi,teta)
    671  
    672   USE exner_hyb_m, ONLY: exner_hyb
    673   USE exner_milieu_m, ONLY: exner_milieu
    674   USE comconst_mod, ONLY: kappa, cpp
    675   USE comvert_mod, ONLY: preff, pressure_exner, bp, ap
    676   USE lmdz_q_sat, ONLY: q_sat
    677   IMPLICIT NONE
    678 
    679   include "dimensions.h"
    680   include "paramet.h"
    681   include "comgeom2.h"
    682 
    683   REAL, DIMENSION (iip1,jjp1),     INTENT(IN) :: psi ! Psol gcm
    684   REAL, DIMENSION (iip1,jjp1,llm), INTENT(IN) :: teta ! Temp. Pot. gcm
    685 
    686   LOGICAL, SAVE                      :: first=.TRUE.
    687   ! Variables pour niveaux pression:
    688   REAL, DIMENSION (iip1,jjp1,nlevnc) :: plnc1,plnc2 !niveaux pression guidage
    689   REAL, DIMENSION (iip1,jjp1,llm)    :: plunc,plsnc !niveaux pression modele
    690   REAL, DIMENSION (iip1,jjm,llm)     :: plvnc       !niveaux pression modele
    691   REAL, DIMENSION (iip1,jjp1,llmp1)  :: p           ! pression intercouches
    692   REAL, DIMENSION (iip1,jjp1,llm)    :: pls, pext   ! var intermediaire
    693   REAL, DIMENSION (iip1,jjp1,llm)    :: pbarx
    694   REAL, DIMENSION (iip1,jjm,llm)     :: pbary
    695   ! Variables pour fonction Exner (P milieu couche)
    696   REAL, DIMENSION (iip1,jjp1,llm)    :: pk
    697   REAL, DIMENSION (iip1,jjp1)        :: pks   
    698   REAL                               :: prefkap,unskap
    699   ! Pression de vapeur saturante
    700   REAL, DIMENSION (ip1jmp1,llm)      :: qsat
    701   !Variables intermediaires interpolation
    702   REAL, DIMENSION (iip1,jjp1,llm)    :: zu1,zu2
    703   REAL, DIMENSION (iip1,jjm,llm)     :: zv1,zv2
    704  
    705   INTEGER                            :: i,j,l,ij
    706   CHARACTER(LEN=20),PARAMETER :: modname="guide_interp"
    707  
    708     WRITE(*,*)trim(modname)//': interpolate nudging variables'
    709 ! -----------------------------------------------------------------
    710 ! Calcul des niveaux de pression champs guidage
    711 ! -----------------------------------------------------------------
    712 IF (guide_modele) THEN
    713     do i=1,iip1
    714         do j=1,jjp1
    715             do l=1,nlevnc
    716                 plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j)
    717                 plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j)
    718             enddo
    719         enddo
    720     enddo
    721 else
    722     do i=1,iip1
    723         do j=1,jjp1
    724             do l=1,nlevnc
    725                 plnc2(i,j,l)=apnc(l)
    726                 plnc1(i,j,l)=apnc(l)
    727            enddo
    728         enddo
    729     enddo
    730 
    731 END IF
    732     IF (first) THEN
    733         first=.FALSE.
    734         WRITE(*,*)trim(modname)//' : check vertical level order'
    735         WRITE(*,*)trim(modname)//' LMDZ :'
    736         do l=1,llm
    737           WRITE(*,*)trim(modname)//' PL(',l,')=',(ap(l)+ap(l+1))/2. &
    738                   +psi(1,jjp1)*(bp(l)+bp(l+1))/2.
    739         enddo
    740         WRITE(*,*)trim(modname)//' nudging file :'
    741         do l=1,nlevnc
    742           WRITE(*,*)trim(modname)//' PL(',l,')=',plnc2(1,1,l)
    743         enddo
    744         WRITE(*,*)trim(modname)//' invert ordering: invert_p=',invert_p
    745         IF (guide_u) THEN
    746             do l=1,nlevnc
    747               WRITE(*,*)trim(modname)//' U(',l,')=',unat2(1,1,l)
    748             enddo
    749         endif
    750         IF (guide_T) THEN
    751             do l=1,nlevnc
    752               WRITE(*,*)trim(modname)//' T(',l,')=',tnat2(1,1,l)
    753             enddo
    754         endif
    755     endif
    756    
    757 ! -----------------------------------------------------------------
    758 ! Calcul niveaux pression modele
    759 ! -----------------------------------------------------------------
    760     CALL pression( ip1jmp1, ap, bp, psi, p )
    761     IF (pressure_exner) THEN
    762       CALL exner_hyb(ip1jmp1,psi,p,pks,pk)
    763     else
    764       CALL exner_milieu(ip1jmp1,psi,p,pks,pk)
    765     endif
    766 !    ....  Calcul de pls , pression au milieu des couches ,en Pascals
    767     unskap=1./kappa
    768     prefkap =  preff  ** kappa
    769     DO l = 1, llm
    770         DO j=1,jjp1
    771             DO i =1, iip1
    772                 pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
    773             ENDDO
    774         ENDDO
    775     ENDDO
    776 
    777 !   calcul des pressions pour les grilles u et v
    778     do l=1,llm
    779         do j=1,jjp1
    780             do i=1,iip1
    781                 pext(i,j,l)=pls(i,j,l)*aire(i,j)
    782             enddo
    783         enddo
    784     enddo
    785     CALL massbar(pext, pbarx, pbary )
    786     do l=1,llm
    787         do j=1,jjp1
    788             do i=1,iip1
    789                 plunc(i,j,l)=pbarx(i,j,l)/aireu(i,j)
    790                 plsnc(i,j,l)=pls(i,j,l)
    791             enddo
    792         enddo
    793     enddo
    794     do l=1,llm
    795         do j=1,jjm
    796             do i=1,iip1
    797                 plvnc(i,j,l)=pbary(i,j,l)/airev(i,j)
    798             enddo
    799         enddo
    800     enddo
    801 
    802 ! -----------------------------------------------------------------
    803 ! Interpolation champs guidage sur niveaux modele (+inversion N/S)
    804 ! Conversion en variables gcm (ucov, vcov...)
    805 ! -----------------------------------------------------------------
    806     IF (guide_P) THEN
    807         do j=1,jjp1
    808             do i=1,iim
    809                 ij=(j-1)*iip1+i
    810                 psgui1(ij)=psnat1(i,j)
    811                 psgui2(ij)=psnat2(i,j)
    812             enddo
    813             psgui1(iip1*j)=psnat1(1,j)
    814             psgui2(iip1*j)=psnat2(1,j)
    815         enddo
    816     endif
    817 
    818     IF (guide_u) THEN
    819         CALL pres2lev(unat1,zu1,nlevnc,llm,plnc1,plunc,iip1,jjp1,invert_p)
    820         CALL pres2lev(unat2,zu2,nlevnc,llm,plnc2,plunc,iip1,jjp1,invert_p)
    821         do l=1,llm
    822             do j=1,jjp1
    823                 do i=1,iim
    824                     ij=(j-1)*iip1+i
    825                     ugui1(ij,l)=zu1(i,j,l)*cu(i,j)
    826                     ugui2(ij,l)=zu2(i,j,l)*cu(i,j)
    827                 enddo
    828                 ugui1(j*iip1,l)=ugui1((j-1)*iip1+1,l)   
    829                 ugui2(j*iip1,l)=ugui2((j-1)*iip1+1,l)   
    830             enddo
    831             do i=1,iip1
    832                 ugui1(i,l)=0.
    833                 ugui1(ip1jm+i,l)=0.
    834                 ugui2(i,l)=0.
    835                 ugui2(ip1jm+i,l)=0.
    836             enddo
    837         enddo
    838     ENDIF
    839    
    840     IF (guide_T) THEN
    841         CALL pres2lev(tnat1,zu1,nlevnc,llm,plnc1,plsnc,iip1,jjp1,invert_p)
    842         CALL pres2lev(tnat2,zu2,nlevnc,llm,plnc2,plsnc,iip1,jjp1,invert_p)
    843         do l=1,llm
    844             do j=1,jjp1
    845                 IF (guide_teta) THEN
    846                     do i=1,iim
    847                         ij=(j-1)*iip1+i
    848                         tgui1(ij,l)=zu1(i,j,l)
    849                         tgui2(ij,l)=zu2(i,j,l)
    850                     enddo
    851                 ELSE
    852                     do i=1,iim
    853                         ij=(j-1)*iip1+i
    854                         tgui1(ij,l)=zu1(i,j,l)*cpp/pk(i,j,l)
    855                         tgui2(ij,l)=zu2(i,j,l)*cpp/pk(i,j,l)
    856                     enddo
    857                 ENDIF
    858                 tgui1(j*iip1,l)=tgui1((j-1)*iip1+1,l)   
    859                 tgui2(j*iip1,l)=tgui2((j-1)*iip1+1,l)   
    860             enddo
    861             do i=1,iip1
    862                 tgui1(i,l)=tgui1(1,l)
    863                 tgui1(ip1jm+i,l)=tgui1(ip1jm+1,l)
    864                 tgui2(i,l)=tgui2(1,l)
    865                 tgui2(ip1jm+i,l)=tgui2(ip1jm+1,l)
    866             enddo
    867         enddo
    868     ENDIF
    869 
    870     IF (guide_v) THEN
    871 
    872         CALL pres2lev(vnat1,zv1,nlevnc,llm,plnc1(:,:jjm,:),plvnc,iip1,jjm,invert_p)
    873         CALL pres2lev(vnat2,zv2,nlevnc,llm,plnc2(:,:jjm,:),plvnc,iip1,jjm,invert_p)
    874 
    875         do l=1,llm
    876             do j=1,jjm
    877                 do i=1,iim
    878                     ij=(j-1)*iip1+i
    879                     vgui1(ij,l)=zv1(i,j,l)*cv(i,j)
    880                     vgui2(ij,l)=zv2(i,j,l)*cv(i,j)
    881                 enddo
    882                 vgui1(j*iip1,l)=vgui1((j-1)*iip1+1,l)   
    883                 vgui2(j*iip1,l)=vgui2((j-1)*iip1+1,l)   
    884             enddo
    885         enddo
    886     ENDIF
    887    
    888     IF (guide_Q) THEN
    889         ! On suppose qu'on a la bonne variable dans le fichier de guidage:
    890         ! Hum.Rel si guide_hr, Hum.Spec. sinon.
    891         CALL pres2lev(qnat1,zu1,nlevnc,llm,plnc1,plsnc,iip1,jjp1,invert_p)
    892         CALL pres2lev(qnat2,zu2,nlevnc,llm,plnc2,plsnc,iip1,jjp1,invert_p)
    893         do l=1,llm
    894             do j=1,jjp1
    895                 do i=1,iim
    896                     ij=(j-1)*iip1+i
    897                     qgui1(ij,l)=zu1(i,j,l)
    898                     qgui2(ij,l)=zu2(i,j,l)
    899                 enddo
    900                 qgui1(j*iip1,l)=qgui1((j-1)*iip1+1,l)   
    901                 qgui2(j*iip1,l)=qgui2((j-1)*iip1+1,l)   
    902             enddo
    903             do i=1,iip1
    904                 qgui1(i,l)=qgui1(1,l)
    905                 qgui1(ip1jm+i,l)=qgui1(ip1jm+1,l)
    906                 qgui2(i,l)=qgui2(1,l)
    907                 qgui2(ip1jm+i,l)=qgui2(ip1jm+1,l)
    908             enddo
    909         enddo
    910         IF (guide_hr) THEN
    911             CALL q_sat(iip1*jjp1*llm,teta*pk/cpp,plsnc,qsat)
    912             qgui1=qgui1*qsat*0.01 !hum. rel. en %
    913             qgui2=qgui2*qsat*0.01
    914         ENDIF
    915     ENDIF
    916 
    917   END SUBROUTINE guide_interp
    918 
    919 !=======================================================================
    920   SUBROUTINE tau2alpha(typ,pim,pjm,factt,taumin,taumax,alpha)
    921 
    922 ! Calcul des constantes de rappel alpha (=1/tau)
    923 
    924     USE comconst_mod, ONLY: pi
    925     USE serre_mod, ONLY: clon, clat, grossismx, grossismy
    926    
     667  !=======================================================================
     668  SUBROUTINE guide_interp(psi, teta)
     669
     670    USE exner_hyb_m, ONLY: exner_hyb
     671    USE exner_milieu_m, ONLY: exner_milieu
     672    USE comconst_mod, ONLY: kappa, cpp
     673    USE comvert_mod, ONLY: preff, pressure_exner, bp, ap
     674    USE lmdz_q_sat, ONLY: q_sat
    927675    IMPLICIT NONE
    928676
    … …  
    931679    include "comgeom2.h"
    932680
    933 ! input arguments :
     681    REAL, DIMENSION (iip1, jjp1), INTENT(IN) :: psi ! Psol gcm
     682    REAL, DIMENSION (iip1, jjp1, llm), INTENT(IN) :: teta ! Temp. Pot. gcm
     683
     684    LOGICAL, SAVE :: first = .TRUE.
     685    ! Variables pour niveaux pression:
     686    REAL, DIMENSION (iip1, jjp1, nlevnc) :: plnc1, plnc2 !niveaux pression guidage
     687    REAL, DIMENSION (iip1, jjp1, llm) :: plunc, plsnc !niveaux pression modele
     688    REAL, DIMENSION (iip1, jjm, llm) :: plvnc       !niveaux pression modele
     689    REAL, DIMENSION (iip1, jjp1, llmp1) :: p           ! pression intercouches
     690    REAL, DIMENSION (iip1, jjp1, llm) :: pls, pext   ! var intermediaire
     691    REAL, DIMENSION (iip1, jjp1, llm) :: pbarx
     692    REAL, DIMENSION (iip1, jjm, llm) :: pbary
     693    ! Variables pour fonction Exner (P milieu couche)
     694    REAL, DIMENSION (iip1, jjp1, llm) :: pk
     695    REAL, DIMENSION (iip1, jjp1) :: pks
     696    REAL :: prefkap, unskap
     697    ! Pression de vapeur saturante
     698    REAL, DIMENSION (ip1jmp1, llm) :: qsat
     699    !Variables intermediaires interpolation
     700    REAL, DIMENSION (iip1, jjp1, llm) :: zu1, zu2
     701    REAL, DIMENSION (iip1, jjm, llm) :: zv1, zv2
     702
     703    INTEGER :: i, j, l, ij
     704    CHARACTER(LEN = 20), PARAMETER :: modname = "guide_interp"
     705
     706    WRITE(*, *)trim(modname) // ': interpolate nudging variables'
     707    ! -----------------------------------------------------------------
     708    ! Calcul des niveaux de pression champs guidage
     709    ! -----------------------------------------------------------------
     710    IF (guide_modele) THEN
     711      do i = 1, iip1
     712        do j = 1, jjp1
     713          do l = 1, nlevnc
     714            plnc2(i, j, l) = apnc(l) + bpnc(l) * psnat2(i, j)
     715            plnc1(i, j, l) = apnc(l) + bpnc(l) * psnat1(i, j)
     716          enddo
     717        enddo
     718      enddo
     719    else
     720      do i = 1, iip1
     721        do j = 1, jjp1
     722          do l = 1, nlevnc
     723            plnc2(i, j, l) = apnc(l)
     724            plnc1(i, j, l) = apnc(l)
     725          enddo
     726        enddo
     727      enddo
     728
     729    END IF
     730    IF (first) THEN
     731      first = .FALSE.
     732      WRITE(*, *)trim(modname) // ' : check vertical level order'
     733      WRITE(*, *)trim(modname) // ' LMDZ :'
     734      do l = 1, llm
     735        WRITE(*, *)trim(modname) // ' PL(', l, ')=', (ap(l) + ap(l + 1)) / 2. &
     736                + psi(1, jjp1) * (bp(l) + bp(l + 1)) / 2.
     737      enddo
     738      WRITE(*, *)trim(modname) // ' nudging file :'
     739      do l = 1, nlevnc
     740        WRITE(*, *)trim(modname) // ' PL(', l, ')=', plnc2(1, 1, l)
     741      enddo
     742      WRITE(*, *)trim(modname) // ' invert ordering: invert_p=', invert_p
     743      IF (guide_u) THEN
     744        do l = 1, nlevnc
     745          WRITE(*, *)trim(modname) // ' U(', l, ')=', unat2(1, 1, l)
     746        enddo
     747      endif
     748      IF (guide_T) THEN
     749        do l = 1, nlevnc
     750          WRITE(*, *)trim(modname) // ' T(', l, ')=', tnat2(1, 1, l)
     751        enddo
     752      endif
     753    endif
     754
     755    ! -----------------------------------------------------------------
     756    ! Calcul niveaux pression modele
     757    ! -----------------------------------------------------------------
     758    CALL pression(ip1jmp1, ap, bp, psi, p)
     759    IF (pressure_exner) THEN
     760      CALL exner_hyb(ip1jmp1, psi, p, pks, pk)
     761    else
     762      CALL exner_milieu(ip1jmp1, psi, p, pks, pk)
     763    endif
     764    !    ....  Calcul de pls , pression au milieu des couches ,en Pascals
     765    unskap = 1. / kappa
     766    prefkap = preff  ** kappa
     767    DO l = 1, llm
     768      DO j = 1, jjp1
     769        DO i = 1, iip1
     770          pls(i, j, l) = preff * (pk(i, j, l) / cpp) ** unskap
     771        ENDDO
     772      ENDDO
     773    ENDDO
     774
     775    !   calcul des pressions pour les grilles u et v
     776    do l = 1, llm
     777      do j = 1, jjp1
     778        do i = 1, iip1
     779          pext(i, j, l) = pls(i, j, l) * aire(i, j)
     780        enddo
     781      enddo
     782    enddo
     783    CALL massbar(pext, pbarx, pbary)
     784    do l = 1, llm
     785      do j = 1, jjp1
     786        do i = 1, iip1
     787          plunc(i, j, l) = pbarx(i, j, l) / aireu(i, j)
     788          plsnc(i, j, l) = pls(i, j, l)
     789        enddo
     790      enddo
     791    enddo
     792    do l = 1, llm
     793      do j = 1, jjm
     794        do i = 1, iip1
     795          plvnc(i, j, l) = pbary(i, j, l) / airev(i, j)
     796        enddo
     797      enddo
     798    enddo
     799
     800    ! -----------------------------------------------------------------
     801    ! Interpolation champs guidage sur niveaux modele (+inversion N/S)
     802    ! Conversion en variables gcm (ucov, vcov...)
     803    ! -----------------------------------------------------------------
     804    IF (guide_P) THEN
     805      do j = 1, jjp1
     806        do i = 1, iim
     807          ij = (j - 1) * iip1 + i
     808          psgui1(ij) = psnat1(i, j)
     809          psgui2(ij) = psnat2(i, j)
     810        enddo
     811        psgui1(iip1 * j) = psnat1(1, j)
     812        psgui2(iip1 * j) = psnat2(1, j)
     813      enddo
     814    endif
     815
     816    IF (guide_u) THEN
     817      CALL pres2lev(unat1, zu1, nlevnc, llm, plnc1, plunc, iip1, jjp1, invert_p)
     818      CALL pres2lev(unat2, zu2, nlevnc, llm, plnc2, plunc, iip1, jjp1, invert_p)
     819      do l = 1, llm
     820        do j = 1, jjp1
     821          do i = 1, iim
     822            ij = (j - 1) * iip1 + i
     823            ugui1(ij, l) = zu1(i, j, l) * cu(i, j)
     824            ugui2(ij, l) = zu2(i, j, l) * cu(i, j)
     825          enddo
     826          ugui1(j * iip1, l) = ugui1((j - 1) * iip1 + 1, l)
     827          ugui2(j * iip1, l) = ugui2((j - 1) * iip1 + 1, l)
     828        enddo
     829        do i = 1, iip1
     830          ugui1(i, l) = 0.
     831          ugui1(ip1jm + i, l) = 0.
     832          ugui2(i, l) = 0.
     833          ugui2(ip1jm + i, l) = 0.
     834        enddo
     835      enddo
     836    ENDIF
     837
     838    IF (guide_T) THEN
     839      CALL pres2lev(tnat1, zu1, nlevnc, llm, plnc1, plsnc, iip1, jjp1, invert_p)
     840      CALL pres2lev(tnat2, zu2, nlevnc, llm, plnc2, plsnc, iip1, jjp1, invert_p)
     841      do l = 1, llm
     842        do j = 1, jjp1
     843          IF (guide_teta) THEN
     844            do i = 1, iim
     845              ij = (j - 1) * iip1 + i
     846              tgui1(ij, l) = zu1(i, j, l)
     847              tgui2(ij, l) = zu2(i, j, l)
     848            enddo
     849          ELSE
     850            do i = 1, iim
     851              ij = (j - 1) * iip1 + i
     852              tgui1(ij, l) = zu1(i, j, l) * cpp / pk(i, j, l)
     853              tgui2(ij, l) = zu2(i, j, l) * cpp / pk(i, j, l)
     854            enddo
     855          ENDIF
     856          tgui1(j * iip1, l) = tgui1((j - 1) * iip1 + 1, l)
     857          tgui2(j * iip1, l) = tgui2((j - 1) * iip1 + 1, l)
     858        enddo
     859        do i = 1, iip1
     860          tgui1(i, l) = tgui1(1, l)
     861          tgui1(ip1jm + i, l) = tgui1(ip1jm + 1, l)
     862          tgui2(i, l) = tgui2(1, l)
     863          tgui2(ip1jm + i, l) = tgui2(ip1jm + 1, l)
     864        enddo
     865      enddo
     866    ENDIF
     867
     868    IF (guide_v) THEN
     869
     870      CALL pres2lev(vnat1, zv1, nlevnc, llm, plnc1(:, :jjm, :), plvnc, iip1, jjm, invert_p)
     871      CALL pres2lev(vnat2, zv2, nlevnc, llm, plnc2(:, :jjm, :), plvnc, iip1, jjm, invert_p)
     872
     873      do l = 1, llm
     874        do j = 1, jjm
     875          do i = 1, iim
     876            ij = (j - 1) * iip1 + i
     877            vgui1(ij, l) = zv1(i, j, l) * cv(i, j)
     878            vgui2(ij, l) = zv2(i, j, l) * cv(i, j)
     879          enddo
     880          vgui1(j * iip1, l) = vgui1((j - 1) * iip1 + 1, l)
     881          vgui2(j * iip1, l) = vgui2((j - 1) * iip1 + 1, l)
     882        enddo
     883      enddo
     884    ENDIF
     885
     886    IF (guide_Q) THEN
     887      ! On suppose qu'on a la bonne variable dans le fichier de guidage:
     888      ! Hum.Rel si guide_hr, Hum.Spec. sinon.
     889      CALL pres2lev(qnat1, zu1, nlevnc, llm, plnc1, plsnc, iip1, jjp1, invert_p)
     890      CALL pres2lev(qnat2, zu2, nlevnc, llm, plnc2, plsnc, iip1, jjp1, invert_p)
     891      do l = 1, llm
     892        do j = 1, jjp1
     893          do i = 1, iim
     894            ij = (j - 1) * iip1 + i
     895            qgui1(ij, l) = zu1(i, j, l)
     896            qgui2(ij, l) = zu2(i, j, l)
     897          enddo
     898          qgui1(j * iip1, l) = qgui1((j - 1) * iip1 + 1, l)
     899          qgui2(j * iip1, l) = qgui2((j - 1) * iip1 + 1, l)
     900        enddo
     901        do i = 1, iip1
     902          qgui1(i, l) = qgui1(1, l)
     903          qgui1(ip1jm + i, l) = qgui1(ip1jm + 1, l)
     904          qgui2(i, l) = qgui2(1, l)
     905          qgui2(ip1jm + i, l) = qgui2(ip1jm + 1, l)
     906        enddo
     907      enddo
     908      IF (guide_hr) THEN
     909        CALL q_sat(iip1 * jjp1 * llm, teta * pk / cpp, plsnc, qsat)
     910        qgui1 = qgui1 * qsat * 0.01 !hum. rel. en %
     911        qgui2 = qgui2 * qsat * 0.01
     912      ENDIF
     913    ENDIF
     914
     915  END SUBROUTINE guide_interp
     916
     917  !=======================================================================
     918  SUBROUTINE tau2alpha(typ, pim, pjm, factt, taumin, taumax, alpha)
     919
     920    ! Calcul des constantes de rappel alpha (=1/tau)
     921
     922    USE comconst_mod, ONLY: pi
     923    USE serre_mod, ONLY: clon, clat, grossismx, grossismy
     924
     925    IMPLICIT NONE
     926
     927    include "dimensions.h"
     928    include "paramet.h"
     929    include "comgeom2.h"
     930
     931    ! input arguments :
    934932    INTEGER, INTENT(IN) :: typ    ! u(2),v(3), ou scalaire(1)
    935     INTEGER, INTENT(IN) :: pim,pjm ! dimensions en lat, lon
    936     REAL, INTENT(IN)    :: factt   ! pas de temps en fraction de jour
    937     REAL, INTENT(IN)    :: taumin,taumax
    938 ! output arguments:
    939     REAL, DIMENSION(pim,pjm), INTENT(OUT) :: alpha
    940  
    941 !  local variables:
    942     LOGICAL, SAVE               :: first=.TRUE.
    943     REAL, SAVE                  :: gamma,dxdy_min,dxdy_max
    944     REAL, DIMENSION (iip1,jjp1) :: zdx,zdy
    945     REAL, DIMENSION (iip1,jjp1) :: dxdys,dxdyu
    946     REAL, DIMENSION (iip1,jjm) :: dxdyv
     933    INTEGER, INTENT(IN) :: pim, pjm ! dimensions en lat, lon
     934    REAL, INTENT(IN) :: factt   ! pas de temps en fraction de jour
     935    REAL, INTENT(IN) :: taumin, taumax
     936    ! output arguments:
     937    REAL, DIMENSION(pim, pjm), INTENT(OUT) :: alpha
     938
     939    !  local variables:
     940    LOGICAL, SAVE :: first = .TRUE.
     941    REAL, SAVE :: gamma, dxdy_min, dxdy_max
     942    REAL, DIMENSION (iip1, jjp1) :: zdx, zdy
     943    REAL, DIMENSION (iip1, jjp1) :: dxdys, dxdyu
     944    REAL, DIMENSION (iip1, jjm) :: dxdyv
    947945    REAL dxdy_
    948     REAL zlat,zlon
    949     REAL alphamin,alphamax,xi
    950     INTEGER i,j,ilon,ilat
    951     CHARACTER(LEN=20),parameter :: modname="tau2alpha"
    952     CHARACTER (len = 80)   :: abort_message
    953 
    954 
    955     alphamin=factt/taumax
    956     alphamax=factt/taumin
     946    REAL zlat, zlon
     947    REAL alphamin, alphamax, xi
     948    INTEGER i, j, ilon, ilat
     949    CHARACTER(LEN = 20), parameter :: modname = "tau2alpha"
     950    CHARACTER (len = 80) :: abort_message
     951
     952    alphamin = factt / taumax
     953    alphamax = factt / taumin
    957954    IF (guide_reg.OR.guide_add) THEN
    958         alpha=alphamax
    959 !-----------------------------------------------------------------------
    960 ! guide_reg: alpha=alpha_min dans region, 0. sinon.
    961 !-----------------------------------------------------------------------
    962         IF (guide_reg) THEN
    963             do j=1,pjm
    964                 do i=1,pim
    965                     IF (typ==2) THEN
    966                        zlat=rlatu(j)*180./pi
    967                        zlon=rlonu(i)*180./pi
    968                     elseif (typ==1) THEN
    969                        zlat=rlatu(j)*180./pi
    970                        zlon=rlonv(i)*180./pi
    971                     elseif (typ==3) THEN
    972                        zlat=rlatv(j)*180./pi
    973                        zlon=rlonv(i)*180./pi
    974                     endif
    975                     alpha(i,j)=alphamax/16.* &
    976                               (1.+tanh((zlat-lat_min_g)/tau_lat))* &
    977                               (1.+tanh((lat_max_g-zlat)/tau_lat))* &
    978                               (1.+tanh((zlon-lon_min_g)/tau_lon))* &
    979                               (1.+tanh((lon_max_g-zlon)/tau_lon))
    980                 enddo
    981             enddo
    982         ENDIF
     955      alpha = alphamax
     956      !-----------------------------------------------------------------------
     957      ! guide_reg: alpha=alpha_min dans region, 0. sinon.
     958      !-----------------------------------------------------------------------
     959      IF (guide_reg) THEN
     960        do j = 1, pjm
     961          do i = 1, pim
     962            IF (typ==2) THEN
     963              zlat = rlatu(j) * 180. / pi
     964              zlon = rlonu(i) * 180. / pi
     965            elseif (typ==1) THEN
     966              zlat = rlatu(j) * 180. / pi
     967              zlon = rlonv(i) * 180. / pi
     968            elseif (typ==3) THEN
     969              zlat = rlatv(j) * 180. / pi
     970              zlon = rlonv(i) * 180. / pi
     971            endif
     972            alpha(i, j) = alphamax / 16. * &
     973                    (1. + tanh((zlat - lat_min_g) / tau_lat)) * &
     974                    (1. + tanh((lat_max_g - zlat) / tau_lat)) * &
     975                    (1. + tanh((zlon - lon_min_g) / tau_lon)) * &
     976                    (1. + tanh((lon_max_g - zlon) / tau_lon))
     977          enddo
     978        enddo
     979      ENDIF
    983980    ELSE
    984 !-----------------------------------------------------------------------
    985 ! Sinon, alpha varie entre alpha_min et alpha_max suivant le zoom.
    986 !-----------------------------------------------------------------------
    987 !Calcul de l'aire des mailles
    988         do j=2,jjm
    989             do i=2,iip1
    990                zdx(i,j)=0.5*(cu(i-1,j)+cu(i,j))/cos(rlatu(j))
    991             enddo
    992             zdx(1,j)=zdx(iip1,j)
    993         enddo
    994         do j=2,jjm
    995             do i=1,iip1
    996                zdy(i,j)=0.5*(cv(i,j-1)+cv(i,j))
    997             enddo
    998         enddo
    999         do i=1,iip1
    1000             zdx(i,1)=zdx(i,2)
    1001             zdx(i,jjp1)=zdx(i,jjm)
    1002             zdy(i,1)=zdy(i,2)
    1003             zdy(i,jjp1)=zdy(i,jjm)
    1004         enddo
    1005         do j=1,jjp1
    1006             do i=1,iip1
    1007                dxdys(i,j)=sqrt(zdx(i,j)*zdx(i,j)+zdy(i,j)*zdy(i,j))
    1008             enddo
    1009         enddo
    1010         IF (typ==2) THEN
    1011             do j=1,jjp1
    1012                 do i=1,iim
    1013                    dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
    1014                 enddo
    1015                 dxdyu(iip1,j)=dxdyu(1,j)
    1016             enddo
    1017         ENDIF
    1018         IF (typ==3) THEN
    1019             do j=1,jjm
    1020                 do i=1,iip1
    1021                    dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i,j+1))
    1022                 enddo
    1023             enddo
    1024         ENDIF
    1025 ! Premier appel: calcul des aires min et max et de gamma.
    1026         IF (first) THEN
    1027             first=.FALSE.
    1028             ! coordonnees du centre du zoom
    1029             CALL coordij(clon,clat,ilon,ilat)
    1030             ! aire de la maille au centre du zoom
    1031             dxdy_min=dxdys(ilon,ilat)
    1032             ! dxdy maximale de la maille
    1033             dxdy_max=0.
    1034             do j=1,jjp1
    1035                 do i=1,iip1
    1036                      dxdy_max=max(dxdy_max,dxdys(i,j))
    1037                 enddo
    1038             enddo
    1039             ! Calcul de gamma
    1040             IF (abs(grossismx-1.)<0.1.OR.abs(grossismy-1.)<0.1) THEN
    1041               WRITE(*,*)trim(modname)//' ATTENTION modele peu zoome'
    1042               WRITE(*,*)trim(modname)//' ATTENTION on prend une constante de guidage cste'
    1043               gamma=0.
     981      !-----------------------------------------------------------------------
     982      ! Sinon, alpha varie entre alpha_min et alpha_max suivant le zoom.
     983      !-----------------------------------------------------------------------
     984      !Calcul de l'aire des mailles
     985      do j = 2, jjm
     986        do i = 2, iip1
     987          zdx(i, j) = 0.5 * (cu(i - 1, j) + cu(i, j)) / cos(rlatu(j))
     988        enddo
     989        zdx(1, j) = zdx(iip1, j)
     990      enddo
     991      do j = 2, jjm
     992        do i = 1, iip1
     993          zdy(i, j) = 0.5 * (cv(i, j - 1) + cv(i, j))
     994        enddo
     995      enddo
     996      do i = 1, iip1
     997        zdx(i, 1) = zdx(i, 2)
     998        zdx(i, jjp1) = zdx(i, jjm)
     999        zdy(i, 1) = zdy(i, 2)
     1000        zdy(i, jjp1) = zdy(i, jjm)
     1001      enddo
     1002      do j = 1, jjp1
     1003        do i = 1, iip1
     1004          dxdys(i, j) = sqrt(zdx(i, j) * zdx(i, j) + zdy(i, j) * zdy(i, j))
     1005        enddo
     1006      enddo
     1007      IF (typ==2) THEN
     1008        do j = 1, jjp1
     1009          do i = 1, iim
     1010            dxdyu(i, j) = 0.5 * (dxdys(i, j) + dxdys(i + 1, j))
     1011          enddo
     1012          dxdyu(iip1, j) = dxdyu(1, j)
     1013        enddo
     1014      ENDIF
     1015      IF (typ==3) THEN
     1016        do j = 1, jjm
     1017          do i = 1, iip1
     1018            dxdyv(i, j) = 0.5 * (dxdys(i, j) + dxdys(i, j + 1))
     1019          enddo
     1020        enddo
     1021      ENDIF
     1022      ! Premier appel: calcul des aires min et max et de gamma.
     1023      IF (first) THEN
     1024        first = .FALSE.
     1025        ! coordonnees du centre du zoom
     1026        CALL coordij(clon, clat, ilon, ilat)
     1027        ! aire de la maille au centre du zoom
     1028        dxdy_min = dxdys(ilon, ilat)
     1029        ! dxdy maximale de la maille
     1030        dxdy_max = 0.
     1031        do j = 1, jjp1
     1032          do i = 1, iip1
     1033            dxdy_max = max(dxdy_max, dxdys(i, j))
     1034          enddo
     1035        enddo
     1036        ! Calcul de gamma
     1037        IF (abs(grossismx - 1.)<0.1.OR.abs(grossismy - 1.)<0.1) THEN
     1038          WRITE(*, *)trim(modname) // ' ATTENTION modele peu zoome'
     1039          WRITE(*, *)trim(modname) // ' ATTENTION on prend une constante de guidage cste'
     1040          gamma = 0.
     1041        else
     1042          gamma = (dxdy_max - 2. * dxdy_min) / (dxdy_max - dxdy_min)
     1043          WRITE(*, *)trim(modname) // ' gamma=', gamma
     1044          IF (gamma<1.e-5) THEN
     1045            WRITE(*, *)trim(modname) // ' gamma =', gamma, '<1e-5'
     1046            abort_message = 'stopped'
     1047            CALL abort_gcm(modname, abort_message, 1)
     1048          endif
     1049          gamma = log(0.5) / log(gamma)
     1050          IF (gamma4) THEN
     1051            gamma = min(gamma, 4.)
     1052          endif
     1053          WRITE(*, *)trim(modname) // ' gamma=', gamma
     1054        endif
     1055      ENDIF !first
     1056
     1057      do j = 1, pjm
     1058        do i = 1, pim
     1059          IF (typ==1) THEN
     1060            dxdy_ = dxdys(i, j)
     1061            zlat = rlatu(j) * 180. / pi
     1062          elseif (typ==2) THEN
     1063            dxdy_ = dxdyu(i, j)
     1064            zlat = rlatu(j) * 180. / pi
     1065          elseif (typ==3) THEN
     1066            dxdy_ = dxdyv(i, j)
     1067            zlat = rlatv(j) * 180. / pi
     1068          endif
     1069          IF (abs(grossismx - 1.)<0.1.OR.abs(grossismy - 1.)<0.1) THEN
     1070            ! pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
     1071            alpha(i, j) = alphamin
     1072          else
     1073            xi = ((dxdy_max - dxdy_) / (dxdy_max - dxdy_min))**gamma
     1074            xi = min(xi, 1.)
     1075            IF(lat_min_g<=zlat .AND. zlat<=lat_max_g) THEN
     1076              alpha(i, j) = xi * alphamin + (1. - xi) * alphamax
    10441077            else
    1045               gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
    1046               WRITE(*,*)trim(modname)//' gamma=',gamma
    1047               IF (gamma<1.e-5) THEN
    1048                 WRITE(*,*)trim(modname)//' gamma =',gamma,'<1e-5'
    1049                 abort_message='stopped'
    1050                 CALL abort_gcm(modname,abort_message,1)
    1051               endif
    1052               gamma=log(0.5)/log(gamma)
    1053               IF (gamma4) THEN
    1054                 gamma=min(gamma,4.)
    1055               endif
    1056               WRITE(*,*)trim(modname)//' gamma=',gamma
     1078              alpha(i, j) = 0.
    10571079            endif
    1058         ENDIF !first
    1059 
    1060         do j=1,pjm
    1061             do i=1,pim
    1062                 IF (typ==1) THEN
    1063                    dxdy_=dxdys(i,j)
    1064                    zlat=rlatu(j)*180./pi
    1065                 elseif (typ==2) THEN
    1066                    dxdy_=dxdyu(i,j)
    1067                    zlat=rlatu(j)*180./pi
    1068                 elseif (typ==3) THEN
    1069                    dxdy_=dxdyv(i,j)
    1070                    zlat=rlatv(j)*180./pi
    1071                 endif
    1072                 IF (abs(grossismx-1.)<0.1.OR.abs(grossismy-1.)<0.1) THEN
    1073                 ! pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
    1074                     alpha(i,j)=alphamin
    1075                 else
    1076                     xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
    1077                     xi=min(xi,1.)
    1078                     IF(lat_min_g<=zlat .AND. zlat<=lat_max_g) THEN
    1079                         alpha(i,j)=xi*alphamin+(1.-xi)*alphamax
    1080                     else
    1081                         alpha(i,j)=0.
    1082                     endif
    1083                 endif
    1084             enddo
    1085         enddo
     1080          endif
     1081        enddo
     1082      enddo
    10861083    ENDIF ! guide_reg
    10871084
    … …  
    10901087  END SUBROUTINE tau2alpha
    10911088
    1092 !=======================================================================
     1089  !=======================================================================
    10931090  SUBROUTINE guide_read(timestep)
    10941091    IMPLICIT NONE
    … …  
    10971094    include "paramet.h"
    10981095
    1099     INTEGER, INTENT(IN)   :: timestep
    1100 
    1101     LOGICAL, SAVE         :: first=.TRUE.
    1102 ! Identification fichiers et variables NetCDF:
    1103     INTEGER, SAVE         :: ncidu,varidu,ncidv,varidv,ncidp,varidp
    1104     INTEGER, SAVE         :: ncidQ,varidQ,ncidt,varidt,ncidps,varidps
    1105     INTEGER               :: ncidpl,varidpl,varidap,varidbp,dimid,lendim
    1106 ! Variables auxiliaires NetCDF:
    1107     INTEGER, DIMENSION(4) :: start,count
    1108     INTEGER               :: status,rcode
    1109     CHARACTER (len = 80)   :: abort_message
    1110     CHARACTER (len = 20)   :: modname = 'guide_read'
    1111     CHARACTER (len = 20)   :: namedim
    1112 
    1113 ! -----------------------------------------------------------------
    1114 ! Premier appel: initialisation de la lecture des fichiers
    1115 ! -----------------------------------------------------------------
     1096    INTEGER, INTENT(IN) :: timestep
     1097
     1098    LOGICAL, SAVE :: first = .TRUE.
     1099    ! Identification fichiers et variables NetCDF:
     1100    INTEGER, SAVE :: ncidu, varidu, ncidv, varidv, ncidp, varidp
     1101    INTEGER, SAVE :: ncidQ, varidQ, ncidt, varidt, ncidps, varidps
     1102    INTEGER :: ncidpl, varidpl, varidap, varidbp, dimid, lendim
     1103    ! Variables auxiliaires NetCDF:
     1104    INTEGER, DIMENSION(4) :: start, count
     1105    INTEGER :: status, rcode
     1106    CHARACTER (len = 80) :: abort_message
     1107    CHARACTER (len = 20) :: modname = 'guide_read'
     1108    CHARACTER (len = 20) :: namedim
     1109
     1110    ! -----------------------------------------------------------------
     1111    ! Premier appel: initialisation de la lecture des fichiers
     1112    ! -----------------------------------------------------------------
    11161113    IF (first) THEN
    1117          ncidpl=-99
    1118          WRITE(*,*) trim(modname)//': opening nudging files '
    1119 ! Niveaux de pression si non constants
    1120          IF (guide_plevs==1) THEN
    1121              WRITE(*,*) trim(modname)//' Reading nudging on model levels'
    1122              rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
    1123              IF (rcode/=nf90_noerr) THEN
    1124               abort_message='Nudging: error -> no file apbp.nc'
    1125               CALL abort_gcm(modname,abort_message,1)
    1126              ENDIF
    1127              rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
    1128              IF (rcode/=nf90_noerr) THEN
    1129               abort_message='Nudging: error -> no AP variable in file apbp.nc'
    1130               CALL abort_gcm(modname,abort_message,1)
    1131              ENDIF
    1132              rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
    1133              IF (rcode/=nf90_noerr) THEN
    1134               abort_message='Nudging: error -> no BP variable in file apbp.nc'
    1135               CALL abort_gcm(modname,abort_message,1)
    1136              ENDIF
    1137              WRITE(*,*) trim(modname)//' ncidpl,varidap',ncidpl,varidap
    1138          endif
    1139 
    1140 ! Pression si guidage sur niveaux P variables
    1141          IF (guide_plevs==2) THEN
    1142              rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
    1143              IF (rcode/=nf90_noerr) THEN
    1144               abort_message='Nudging: error -> no file P.nc'
    1145               CALL abort_gcm(modname,abort_message,1)
    1146              ENDIF
    1147              rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
    1148              IF (rcode/=nf90_noerr) THEN
    1149               abort_message='Nudging: error -> no PRES variable in file P.nc'
    1150               CALL abort_gcm(modname,abort_message,1)
    1151              ENDIF
    1152              WRITE(*,*) trim(modname)//' ncidp,varidp',ncidp,varidp
    1153              IF (ncidpl==-99) ncidpl=ncidp
    1154          endif
    1155 
    1156 ! Vent zonal
    1157          IF (guide_u) THEN
    1158              rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
    1159              IF (rcode/=nf90_noerr) THEN
    1160               abort_message='Nudging: error -> no file u.nc'
    1161               CALL abort_gcm(modname,abort_message,1)
    1162              ENDIF
    1163              rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
    1164              IF (rcode/=nf90_noerr) THEN
    1165               abort_message='Nudging: error -> no UWND variable in file u.nc'
    1166               CALL abort_gcm(modname,abort_message,1)
    1167              ENDIF
    1168              WRITE(*,*) trim(modname)//' ncidu,varidu',ncidu,varidu
    1169              IF (ncidpl==-99) ncidpl=ncidu
    1170 
    1171              status=nf90_inq_dimid(ncidu, "LONU", dimid)
    1172              status=nf90_inquire_dimension(ncidu,dimid,namedim,lendim)
    1173              IF (lendim /= iip1) THEN
    1174                 abort_message='dimension LONU different from iip1 in u.nc'
    1175                 CALL abort_gcm(modname,abort_message,1)
    1176              ENDIF
    1177 
    1178              status=nf90_inq_dimid(ncidu, "LATU", dimid)
    1179              status=nf90_inquire_dimension(ncidu,dimid,namedim,lendim)
    1180              IF (lendim /= jjp1) THEN
    1181                 abort_message='dimension LATU different from jjp1 in u.nc'
    1182                 CALL abort_gcm(modname,abort_message,1)
    1183              ENDIF
    1184 
    1185          endif
    1186 
    1187 ! Vent meridien
    1188          IF (guide_v) THEN
    1189              rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
    1190              IF (rcode/=nf90_noerr) THEN
    1191               abort_message='Nudging: error -> no file v.nc'
    1192               CALL abort_gcm(modname,abort_message,1)
    1193              ENDIF
    1194              rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
    1195              IF (rcode/=nf90_noerr) THEN
    1196               abort_message='Nudging: error -> no VWND variable in file v.nc'
    1197               CALL abort_gcm(modname,abort_message,1)
    1198              ENDIF
    1199              WRITE(*,*) trim(modname)//' ncidv,varidv',ncidv,varidv
    1200              IF (ncidpl==-99) ncidpl=ncidv
    1201              
    1202              status=nf90_inq_dimid(ncidv, "LONV", dimid)
    1203              status=nf90_inquire_dimension(ncidv,dimid,namedim,lendim)
    1204              
    1205                 IF (lendim /= iip1) THEN
    1206                 abort_message='dimension LONV different from iip1 in v.nc'
    1207                 CALL abort_gcm(modname,abort_message,1)
    1208              ENDIF
    1209 
    1210 
    1211              status=nf90_inq_dimid(ncidv, "LATV", dimid)
    1212              status=nf90_inquire_dimension(ncidv,dimid,namedim,lendim)
    1213              IF (lendim /= jjm) THEN
    1214                 abort_message='dimension LATV different from jjm in v.nc'
    1215                 CALL abort_gcm(modname,abort_message,1)
    1216              ENDIF
    1217        
    1218          endif
    1219 
    1220 ! Temperature
    1221          IF (guide_T) THEN
    1222              rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
    1223              IF (rcode/=nf90_noerr) THEN
    1224               abort_message='Nudging: error -> no file T.nc'
    1225               CALL abort_gcm(modname,abort_message,1)
    1226              ENDIF
    1227              rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
    1228              IF (rcode/=nf90_noerr) THEN
    1229               abort_message='Nudging: error -> no AIR variable in file T.nc'
    1230               CALL abort_gcm(modname,abort_message,1)
    1231              ENDIF
    1232              WRITE(*,*) trim(modname)//' ncidT,varidT',ncidt,varidt
    1233              IF (ncidpl==-99) ncidpl=ncidt
    1234 
    1235              status=nf90_inq_dimid(ncidt, "LONV", dimid)
    1236              status=nf90_inquire_dimension(ncidt,dimid,namedim,lendim)
    1237              IF (lendim /= iip1) THEN
    1238                 abort_message='dimension LONV different from iip1 in T.nc'
    1239                 CALL abort_gcm(modname,abort_message,1)
    1240              ENDIF
    1241 
    1242              status=nf90_inq_dimid(ncidt, "LATU", dimid)
    1243              status=nf90_inquire_dimension(ncidt,dimid,namedim,lendim)
    1244              IF (lendim /= jjp1) THEN
    1245                 abort_message='dimension LATU different from jjp1 in T.nc'
    1246                 CALL abort_gcm(modname,abort_message,1)
    1247              ENDIF
    1248 
    1249          endif
    1250 
    1251 ! Humidite
    1252          IF (guide_Q) THEN
    1253              rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
    1254              IF (rcode/=nf90_noerr) THEN
    1255               abort_message='Nudging: error -> no file hur.nc'
    1256               CALL abort_gcm(modname,abort_message,1)
    1257              ENDIF
    1258              rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
    1259              IF (rcode/=nf90_noerr) THEN
    1260               abort_message='Nudging: error -> no RH variable in file hur.nc'
    1261               CALL abort_gcm(modname,abort_message,1)
    1262              ENDIF
    1263              WRITE(*,*) trim(modname)//' ncidQ,varidQ',ncidQ,varidQ
    1264              IF (ncidpl==-99) ncidpl=ncidQ
    1265 
    1266              status=nf90_inq_dimid(ncidQ, "LONV", dimid)
    1267              status=nf90_inquire_dimension(ncidQ,dimid,namedim,lendim)
    1268              IF (lendim /= iip1) THEN
    1269                 abort_message='dimension LONV different from iip1 in hur.nc'
    1270                 CALL abort_gcm(modname,abort_message,1)
    1271              ENDIF
    1272 
    1273              status=nf90_inq_dimid(ncidQ, "LATU", dimid)
    1274              status=nf90_inquire_dimension(ncidQ,dimid,namedim,lendim)
    1275              IF (lendim /= jjp1) THEN
    1276                 abort_message='dimension LATU different from jjp1 in hur.nc'
    1277                 CALL abort_gcm(modname,abort_message,1)
    1278              ENDIF
    1279 
    1280          endif
    1281 
    1282 ! Pression de surface
    1283          IF ((guide_P).OR.(guide_modele)) THEN
    1284              rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
    1285              IF (rcode/=nf90_noerr) THEN
    1286               abort_message='Nudging: error -> no file ps.nc'
    1287               CALL abort_gcm(modname,abort_message,1)
    1288              ENDIF
    1289              rcode = nf90_inq_varid(ncidps, 'SP', varidps)
    1290              IF (rcode/=nf90_noerr) THEN
    1291               abort_message='Nudging: error -> no SP variable in file ps.nc'
    1292               CALL abort_gcm(modname,abort_message,1)
    1293              ENDIF
    1294              WRITE(*,*) trim(modname)//' ncidps,varidps',ncidps,varidps
    1295          endif
    1296 ! Coordonnee verticale
    1297          IF (guide_plevs==0) THEN
    1298               rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
    1299               IF (rcode/=0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
    1300               WRITE(*,*) trim(modname)//' ncidpl,varidpl',ncidpl,varidpl
    1301          endif
    1302 ! Coefs ap, bp pour calcul de la pression aux differents niveaux
    1303          IF (guide_plevs==1) THEN
    1304              status=nf90_get_var(ncidpl,varidap,apnc,[1],[nlevnc])
    1305              status=nf90_get_var(ncidpl,varidbp,bpnc,[1],[nlevnc])
    1306          ELSEIF (guide_plevs==0) THEN
    1307              status=nf90_get_var(ncidpl,varidpl,apnc,[1],[nlevnc])
    1308 !FC Pour les corrections la pression est deja en Pascals on commente la ligne ci-dessous
    1309              IF(convert_Pa) apnc=apnc*100.! conversion en Pascals
    1310              bpnc(:)=0.
    1311          endif
    1312          first=.FALSE.
    1313      endif ! (first)
    1314 
    1315 ! -----------------------------------------------------------------
    1316 !   lecture des champs u, v, T, Q, ps
    1317 ! -----------------------------------------------------------------
    1318 
    1319 !  dimensions pour les champs scalaires et le vent zonal
    1320      start(1)=1
    1321      start(2)=1
    1322      start(3)=1
    1323      start(4)=timestep
    1324 
    1325      count(1)=iip1
    1326      count(2)=jjp1
    1327      count(3)=nlevnc
    1328      count(4)=1
    1329 
    1330 ! Pression
    1331      IF (guide_plevs==2) THEN
    1332          status=nf90_get_var(ncidp,varidp,pnat2,start,count)
    1333          IF (invert_y) THEN
    1334 !           PRINT*,"Invertion impossible actuellement"
    1335 !           CALL abort_gcm(modname,abort_message,1)
    1336            CALL invert_lat(iip1,jjp1,nlevnc,pnat2)
    1337          ENDIF
    1338      endif
    1339 
    1340 !  Vent zonal
    1341      IF (guide_u) THEN
    1342          status=nf90_get_var(ncidu,varidu,unat2,start,count)
    1343          IF (invert_y) THEN
    1344            CALL invert_lat(iip1,jjp1,nlevnc,unat2)
    1345          ENDIF
    1346      endif
    1347 
    1348 !  Temperature
    1349      IF (guide_T) THEN
    1350          status=nf90_get_var(ncidt,varidt,tnat2,start,count)
    1351          IF (invert_y) THEN
    1352            CALL invert_lat(iip1,jjp1,nlevnc,tnat2)
    1353          ENDIF
    1354      endif
    1355 
    1356 !  Humidite
    1357      IF (guide_Q) THEN
    1358          status=nf90_get_var(ncidQ,varidQ,qnat2,start,count)
    1359          IF (invert_y) THEN
    1360            CALL invert_lat(iip1,jjp1,nlevnc,qnat2)
    1361          ENDIF
    1362          
    1363      endif
    1364 
    1365 !  Vent meridien
    1366      IF (guide_v) THEN
    1367          count(2)=jjm
    1368          status=nf90_get_var(ncidv,varidv,vnat2,start,count)
    1369          IF (invert_y) THEN
    1370            CALL invert_lat(iip1,jjm,nlevnc,vnat2)
    1371          ENDIF
    1372      endif
    1373 
    1374 !  Pression de surface
    1375      IF ((guide_P).OR.(guide_modele))  THEN
    1376          start(3)=timestep
    1377          start(4)=0
    1378          count(2)=jjp1
    1379          count(3)=1
    1380          count(4)=0
    1381          status=nf90_get_var(ncidps,varidps,psnat2,start,count)
    1382          IF (invert_y) THEN
    1383            CALL invert_lat(iip1,jjp1,1,psnat2)
    1384          ENDIF
    1385      endif
     1114      ncidpl = -99
     1115      WRITE(*, *) trim(modname) // ': opening nudging files '
     1116      ! Niveaux de pression si non constants
     1117      IF (guide_plevs==1) THEN
     1118        WRITE(*, *) trim(modname) // ' Reading nudging on model levels'
     1119        rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
     1120        IF (rcode/=nf90_noerr) THEN
     1121          abort_message = 'Nudging: error -> no file apbp.nc'
     1122          CALL abort_gcm(modname, abort_message, 1)
     1123        ENDIF
     1124        rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
     1125        IF (rcode/=nf90_noerr) THEN
     1126          abort_message = 'Nudging: error -> no AP variable in file apbp.nc'
     1127          CALL abort_gcm(modname, abort_message, 1)
     1128        ENDIF
     1129        rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
     1130        IF (rcode/=nf90_noerr) THEN
     1131          abort_message = 'Nudging: error -> no BP variable in file apbp.nc'
     1132          CALL abort_gcm(modname, abort_message, 1)
     1133        ENDIF
     1134        WRITE(*, *) trim(modname) // ' ncidpl,varidap', ncidpl, varidap
     1135      endif
     1136
     1137      ! Pression si guidage sur niveaux P variables
     1138      IF (guide_plevs==2) THEN
     1139        rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
     1140        IF (rcode/=nf90_noerr) THEN
     1141          abort_message = 'Nudging: error -> no file P.nc'
     1142          CALL abort_gcm(modname, abort_message, 1)
     1143        ENDIF
     1144        rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
     1145        IF (rcode/=nf90_noerr) THEN
     1146          abort_message = 'Nudging: error -> no PRES variable in file P.nc'
     1147          CALL abort_gcm(modname, abort_message, 1)
     1148        ENDIF
     1149        WRITE(*, *) trim(modname) // ' ncidp,varidp', ncidp, varidp
     1150        IF (ncidpl==-99) ncidpl = ncidp
     1151      endif
     1152
     1153      ! Vent zonal
     1154      IF (guide_u) THEN
     1155        rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
     1156        IF (rcode/=nf90_noerr) THEN
     1157          abort_message = 'Nudging: error -> no file u.nc'
     1158          CALL abort_gcm(modname, abort_message, 1)
     1159        ENDIF
     1160        rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
     1161        IF (rcode/=nf90_noerr) THEN
     1162          abort_message = 'Nudging: error -> no UWND variable in file u.nc'
     1163          CALL abort_gcm(modname, abort_message, 1)
     1164        ENDIF
     1165        WRITE(*, *) trim(modname) // ' ncidu,varidu', ncidu, varidu
     1166        IF (ncidpl==-99) ncidpl = ncidu
     1167
     1168        status = nf90_inq_dimid(ncidu, "LONU", dimid)
     1169        status = nf90_inquire_dimension(ncidu, dimid, namedim, lendim)
     1170        IF (lendim /= iip1) THEN
     1171          abort_message = 'dimension LONU different from iip1 in u.nc'
     1172          CALL abort_gcm(modname, abort_message, 1)
     1173        ENDIF
     1174
     1175        status = nf90_inq_dimid(ncidu, "LATU", dimid)
     1176        status = nf90_inquire_dimension(ncidu, dimid, namedim, lendim)
     1177        IF (lendim /= jjp1) THEN
     1178          abort_message = 'dimension LATU different from jjp1 in u.nc'
     1179          CALL abort_gcm(modname, abort_message, 1)
     1180        ENDIF
     1181
     1182      endif
     1183
     1184      ! Vent meridien
     1185      IF (guide_v) THEN
     1186        rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
     1187        IF (rcode/=nf90_noerr) THEN
     1188          abort_message = 'Nudging: error -> no file v.nc'
     1189          CALL abort_gcm(modname, abort_message, 1)
     1190        ENDIF
     1191        rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
     1192        IF (rcode/=nf90_noerr) THEN
     1193          abort_message = 'Nudging: error -> no VWND variable in file v.nc'
     1194          CALL abort_gcm(modname, abort_message, 1)
     1195        ENDIF
     1196        WRITE(*, *) trim(modname) // ' ncidv,varidv', ncidv, varidv
     1197        IF (ncidpl==-99) ncidpl = ncidv
     1198
     1199        status = nf90_inq_dimid(ncidv, "LONV", dimid)
     1200        status = nf90_inquire_dimension(ncidv, dimid, namedim, lendim)
     1201
     1202        IF (lendim /= iip1) THEN
     1203          abort_message = 'dimension LONV different from iip1 in v.nc'
     1204          CALL abort_gcm(modname, abort_message, 1)
     1205        ENDIF
     1206
     1207        status = nf90_inq_dimid(ncidv, "LATV", dimid)
     1208        status = nf90_inquire_dimension(ncidv, dimid, namedim, lendim)
     1209        IF (lendim /= jjm) THEN
     1210          abort_message = 'dimension LATV different from jjm in v.nc'
     1211          CALL abort_gcm(modname, abort_message, 1)
     1212        ENDIF
     1213
     1214      endif
     1215
     1216      ! Temperature
     1217      IF (guide_T) THEN
     1218        rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
     1219        IF (rcode/=nf90_noerr) THEN
     1220          abort_message = 'Nudging: error -> no file T.nc'
     1221          CALL abort_gcm(modname, abort_message, 1)
     1222        ENDIF
     1223        rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
     1224        IF (rcode/=nf90_noerr) THEN
     1225          abort_message = 'Nudging: error -> no AIR variable in file T.nc'
     1226          CALL abort_gcm(modname, abort_message, 1)
     1227        ENDIF
     1228        WRITE(*, *) trim(modname) // ' ncidT,varidT', ncidt, varidt
     1229        IF (ncidpl==-99) ncidpl = ncidt
     1230
     1231        status = nf90_inq_dimid(ncidt, "LONV", dimid)
     1232        status = nf90_inquire_dimension(ncidt, dimid, namedim, lendim)
     1233        IF (lendim /= iip1) THEN
     1234          abort_message = 'dimension LONV different from iip1 in T.nc'
     1235          CALL abort_gcm(modname, abort_message, 1)
     1236        ENDIF
     1237
     1238        status = nf90_inq_dimid(ncidt, "LATU", dimid)
     1239        status = nf90_inquire_dimension(ncidt, dimid, namedim, lendim)
     1240        IF (lendim /= jjp1) THEN
     1241          abort_message = 'dimension LATU different from jjp1 in T.nc'
     1242          CALL abort_gcm(modname, abort_message, 1)
     1243        ENDIF
     1244
     1245      endif
     1246
     1247      ! Humidite
     1248      IF (guide_Q) THEN
     1249        rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
     1250        IF (rcode/=nf90_noerr) THEN
     1251          abort_message = 'Nudging: error -> no file hur.nc'
     1252          CALL abort_gcm(modname, abort_message, 1)
     1253        ENDIF
     1254        rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
     1255        IF (rcode/=nf90_noerr) THEN
     1256          abort_message = 'Nudging: error -> no RH variable in file hur.nc'
     1257          CALL abort_gcm(modname, abort_message, 1)
     1258        ENDIF
     1259        WRITE(*, *) trim(modname) // ' ncidQ,varidQ', ncidQ, varidQ
     1260        IF (ncidpl==-99) ncidpl = ncidQ
     1261
     1262        status = nf90_inq_dimid(ncidQ, "LONV", dimid)
     1263        status = nf90_inquire_dimension(ncidQ, dimid, namedim, lendim)
     1264        IF (lendim /= iip1) THEN
     1265          abort_message = 'dimension LONV different from iip1 in hur.nc'
     1266          CALL abort_gcm(modname, abort_message, 1)
     1267        ENDIF
     1268
     1269        status = nf90_inq_dimid(ncidQ, "LATU", dimid)
     1270        status = nf90_inquire_dimension(ncidQ, dimid, namedim, lendim)
     1271        IF (lendim /= jjp1) THEN
     1272          abort_message = 'dimension LATU different from jjp1 in hur.nc'
     1273          CALL abort_gcm(modname, abort_message, 1)
     1274        ENDIF
     1275
     1276      endif
     1277
     1278      ! Pression de surface
     1279      IF ((guide_P).OR.(guide_modele)) THEN
     1280        rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
     1281        IF (rcode/=nf90_noerr) THEN
     1282          abort_message = 'Nudging: error -> no file ps.nc'
     1283          CALL abort_gcm(modname, abort_message, 1)
     1284        ENDIF
     1285        rcode = nf90_inq_varid(ncidps, 'SP', varidps)
     1286        IF (rcode/=nf90_noerr) THEN
     1287          abort_message = 'Nudging: error -> no SP variable in file ps.nc'
     1288          CALL abort_gcm(modname, abort_message, 1)
     1289        ENDIF
     1290        WRITE(*, *) trim(modname) // ' ncidps,varidps', ncidps, varidps
     1291      endif
     1292      ! Coordonnee verticale
     1293      IF (guide_plevs==0) THEN
     1294        rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
     1295        IF (rcode/=0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
     1296        WRITE(*, *) trim(modname) // ' ncidpl,varidpl', ncidpl, varidpl
     1297      endif
     1298      ! Coefs ap, bp pour calcul de la pression aux differents niveaux
     1299      IF (guide_plevs==1) THEN
     1300        status = nf90_get_var(ncidpl, varidap, apnc, [1], [nlevnc])
     1301        status = nf90_get_var(ncidpl, varidbp, bpnc, [1], [nlevnc])
     1302      ELSEIF (guide_plevs==0) THEN
     1303        status = nf90_get_var(ncidpl, varidpl, apnc, [1], [nlevnc])
     1304        !FC Pour les corrections la pression est deja en Pascals on commente la ligne ci-dessous
     1305        IF(convert_Pa) apnc = apnc * 100.! conversion en Pascals
     1306        bpnc(:) = 0.
     1307      endif
     1308      first = .FALSE.
     1309    endif ! (first)
     1310
     1311    ! -----------------------------------------------------------------
     1312    !   lecture des champs u, v, T, Q, ps
     1313    ! -----------------------------------------------------------------
     1314
     1315    !  dimensions pour les champs scalaires et le vent zonal
     1316    start(1) = 1
     1317    start(2) = 1
     1318    start(3) = 1
     1319    start(4) = timestep
     1320
     1321    count(1) = iip1
     1322    count(2) = jjp1
     1323    count(3) = nlevnc
     1324    count(4) = 1
     1325
     1326    ! Pression
     1327    IF (guide_plevs==2) THEN
     1328      status = nf90_get_var(ncidp, varidp, pnat2, start, count)
     1329      IF (invert_y) THEN
     1330        !           PRINT*,"Invertion impossible actuellement"
     1331        !           CALL abort_gcm(modname,abort_message,1)
     1332        CALL invert_lat(iip1, jjp1, nlevnc, pnat2)
     1333      ENDIF
     1334    endif
     1335
     1336    !  Vent zonal
     1337    IF (guide_u) THEN
     1338      status = nf90_get_var(ncidu, varidu, unat2, start, count)
     1339      IF (invert_y) THEN
     1340        CALL invert_lat(iip1, jjp1, nlevnc, unat2)
     1341      ENDIF
     1342    endif
     1343
     1344    !  Temperature
     1345    IF (guide_T) THEN
     1346      status = nf90_get_var(ncidt, varidt, tnat2, start, count)
     1347      IF (invert_y) THEN
     1348        CALL invert_lat(iip1, jjp1, nlevnc, tnat2)
     1349      ENDIF
     1350    endif
     1351
     1352    !  Humidite
     1353    IF (guide_Q) THEN
     1354      status = nf90_get_var(ncidQ, varidQ, qnat2, start, count)
     1355      IF (invert_y) THEN
     1356        CALL invert_lat(iip1, jjp1, nlevnc, qnat2)
     1357      ENDIF
     1358
     1359    endif
     1360
     1361    !  Vent meridien
     1362    IF (guide_v) THEN
     1363      count(2) = jjm
     1364      status = nf90_get_var(ncidv, varidv, vnat2, start, count)
     1365      IF (invert_y) THEN
     1366        CALL invert_lat(iip1, jjm, nlevnc, vnat2)
     1367      ENDIF
     1368    endif
     1369
     1370    !  Pression de surface
     1371    IF ((guide_P).OR.(guide_modele))  THEN
     1372      start(3) = timestep
     1373      start(4) = 0
     1374      count(2) = jjp1
     1375      count(3) = 1
     1376      count(4) = 0
     1377      status = nf90_get_var(ncidps, varidps, psnat2, start, count)
     1378      IF (invert_y) THEN
     1379        CALL invert_lat(iip1, jjp1, 1, psnat2)
     1380      ENDIF
     1381    endif
    13861382
    13871383  END SUBROUTINE guide_read
    13881384
    1389 !=======================================================================
     1385  !=======================================================================
    13901386  SUBROUTINE guide_read2D(timestep)
    13911387    IMPLICIT NONE
    … …  
    13941390    include "paramet.h"
    13951391
    1396     INTEGER, INTENT(IN)   :: timestep
    1397 
    1398     LOGICAL, SAVE         :: first=.TRUE.
    1399 ! Identification fichiers et variables NetCDF:
    1400     INTEGER, SAVE         :: ncidu,varidu,ncidv,varidv,ncidp,varidp
    1401     INTEGER, SAVE         :: ncidQ,varidQ,ncidt,varidt,ncidps,varidps
    1402     INTEGER               :: ncidpl,varidpl,varidap,varidbp
    1403 ! Variables auxiliaires NetCDF:
    1404     INTEGER, DIMENSION(4) :: start,count
    1405     INTEGER               :: status,rcode
    1406 ! Variables for 3D extension:
    1407     REAL, DIMENSION (jjp1,llm) :: zu
    1408     REAL, DIMENSION (jjm,llm) :: zv
    1409     INTEGER               :: i
    1410 
    1411     CHARACTER (len = 80)   :: abort_message
    1412     CHARACTER (len = 20)   :: modname = 'guide_read2D'
    1413 ! -----------------------------------------------------------------
    1414 ! Premier appel: initialisation de la lecture des fichiers
    1415 ! -----------------------------------------------------------------
     1392    INTEGER, INTENT(IN) :: timestep
     1393
     1394    LOGICAL, SAVE :: first = .TRUE.
     1395    ! Identification fichiers et variables NetCDF:
     1396    INTEGER, SAVE :: ncidu, varidu, ncidv, varidv, ncidp, varidp
     1397    INTEGER, SAVE :: ncidQ, varidQ, ncidt, varidt, ncidps, varidps
     1398    INTEGER :: ncidpl, varidpl, varidap, varidbp
     1399    ! Variables auxiliaires NetCDF:
     1400    INTEGER, DIMENSION(4) :: start, count
     1401    INTEGER :: status, rcode
     1402    ! Variables for 3D extension:
     1403    REAL, DIMENSION (jjp1, llm) :: zu
     1404    REAL, DIMENSION (jjm, llm) :: zv
     1405    INTEGER :: i
     1406
     1407    CHARACTER (len = 80) :: abort_message
     1408    CHARACTER (len = 20) :: modname = 'guide_read2D'
     1409    ! -----------------------------------------------------------------
     1410    ! Premier appel: initialisation de la lecture des fichiers
     1411    ! -----------------------------------------------------------------
    14161412    IF (first) THEN
    1417          ncidpl=-99
    1418          WRITE(*,*)trim(modname)//' : opening nudging files '
    1419 ! Ap et Bp si niveaux de pression hybrides
    1420          IF (guide_plevs==1) THEN
    1421            WRITE(*,*)trim(modname)//' Reading nudging on model levels'
    1422            rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
    1423            IF (rcode/=nf90_noerr) THEN
    1424              abort_message='Nudging: error -> no file apbp.nc'
    1425            CALL abort_gcm(modname,abort_message,1)
    1426            ENDIF
    1427            rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
    1428            IF (rcode/=nf90_noerr) THEN
    1429              abort_message='Nudging: error -> no AP variable in file apbp.nc'
    1430            CALL abort_gcm(modname,abort_message,1)
    1431            ENDIF
    1432            rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
    1433            IF (rcode/=nf90_noerr) THEN
    1434              abort_message='Nudging: error -> no BP variable in file apbp.nc'
    1435              CALL abort_gcm(modname,abort_message,1)
    1436            ENDIF
    1437            WRITE(*,*)trim(modname)//'ncidpl,varidap',ncidpl,varidap
    1438          endif
    1439 ! Pression
    1440          IF (guide_plevs==2) THEN
    1441            rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
    1442            IF (rcode/=nf90_noerr) THEN
    1443              abort_message='Nudging: error -> no file P.nc'
    1444              CALL abort_gcm(modname,abort_message,1)
    1445            ENDIF
    1446            rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
    1447            IF (rcode/=nf90_noerr) THEN
    1448              abort_message='Nudging: error -> no PRES variable in file P.nc'
    1449              CALL abort_gcm(modname,abort_message,1)
    1450            ENDIF
    1451            WRITE(*,*)trim(modname)//' ncidp,varidp',ncidp,varidp
    1452            IF (ncidpl==-99) ncidpl=ncidp
    1453          endif
    1454 ! Vent zonal
    1455          IF (guide_u) THEN
    1456            rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
    1457            IF (rcode/=nf90_noerr) THEN
    1458              abort_message='Nudging: error -> no file u.nc'
    1459              CALL abort_gcm(modname,abort_message,1)
    1460            ENDIF
    1461            rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
    1462            IF (rcode/=nf90_noerr) THEN
    1463              abort_message='Nudging: error -> no UWND variable in file u.nc'
    1464              CALL abort_gcm(modname,abort_message,1)
    1465            ENDIF
    1466            WRITE(*,*)trim(modname)//' ncidu,varidu',ncidu,varidu
    1467            IF (ncidpl==-99) ncidpl=ncidu
    1468          endif
    1469 ! Vent meridien
    1470          IF (guide_v) THEN
    1471            rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
    1472            IF (rcode/=nf90_noerr) THEN
    1473              abort_message='Nudging: error -> no file v.nc'
    1474              CALL abort_gcm(modname,abort_message,1)
    1475            ENDIF
    1476            rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
    1477            IF (rcode/=nf90_noerr) THEN
    1478              abort_message='Nudging: error -> no VWND variable in file v.nc'
    1479              CALL abort_gcm(modname,abort_message,1)
    1480            ENDIF
    1481            WRITE(*,*)trim(modname)//' ncidv,varidv',ncidv,varidv
    1482            IF (ncidpl==-99) ncidpl=ncidv
    1483          endif
    1484 ! Temperature
    1485          IF (guide_T) THEN
    1486            rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
    1487            IF (rcode/=nf90_noerr) THEN
    1488              abort_message='Nudging: error -> no file T.nc'
    1489              CALL abort_gcm(modname,abort_message,1)
    1490            ENDIF
    1491            rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
    1492            IF (rcode/=nf90_noerr) THEN
    1493              abort_message='Nudging: error -> no AIR variable in file T.nc'
    1494              CALL abort_gcm(modname,abort_message,1)
    1495            ENDIF
    1496            WRITE(*,*)trim(modname)//' ncidT,varidT',ncidt,varidt
    1497            IF (ncidpl==-99) ncidpl=ncidt
    1498          endif
    1499 ! Humidite
    1500          IF (guide_Q) THEN
    1501            rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
    1502            IF (rcode/=nf90_noerr) THEN
    1503              abort_message='Nudging: error -> no file hur.nc'
    1504              CALL abort_gcm(modname,abort_message,1)
    1505            ENDIF
    1506            rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
    1507            IF (rcode/=nf90_noerr) THEN
    1508              abort_message='Nudging: error -> no RH,variable in file hur.nc'
    1509              CALL abort_gcm(modname,abort_message,1)
    1510            ENDIF
    1511            WRITE(*,*)trim(modname)//' ncidQ,varidQ',ncidQ,varidQ
    1512            IF (ncidpl==-99) ncidpl=ncidQ
    1513          endif
    1514 ! Pression de surface
    1515          IF ((guide_P).OR.(guide_modele)) THEN
    1516            rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
    1517            IF (rcode/=nf90_noerr) THEN
    1518              abort_message='Nudging: error -> no file ps.nc'
    1519              CALL abort_gcm(modname,abort_message,1)
    1520            ENDIF
    1521            rcode = nf90_inq_varid(ncidps, 'SP', varidps)
    1522            IF (rcode/=nf90_noerr) THEN
    1523              abort_message='Nudging: error -> no SP variable in file ps.nc'
    1524              CALL abort_gcm(modname,abort_message,1)
    1525            ENDIF
    1526            WRITE(*,*)trim(modname)//' ncidps,varidps',ncidps,varidps
    1527          endif
    1528 ! Coordonnee verticale
    1529          IF (guide_plevs==0) THEN
    1530            rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
    1531            IF (rcode/=0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
    1532            WRITE(*,*)trim(modname)//' ncidpl,varidpl',ncidpl,varidpl
    1533          endif
    1534 ! Coefs ap, bp pour calcul de la pression aux differents niveaux
    1535          IF (guide_plevs==1) THEN
    1536              status=nf90_get_var(ncidpl,varidap,apnc,[1],[nlevnc])
    1537              status=nf90_get_var(ncidpl,varidbp,bpnc,[1],[nlevnc])
    1538          elseif (guide_plevs==0) THEN
    1539              status=nf90_get_var(ncidpl,varidpl,apnc,[1],[nlevnc])
    1540              apnc=apnc*100.! conversion en Pascals
    1541              bpnc(:)=0.
    1542          endif
    1543          first=.FALSE.
    1544      endif ! (first)
    1545 
    1546 ! -----------------------------------------------------------------
    1547 !   lecture des champs u, v, T, Q, ps
    1548 ! -----------------------------------------------------------------
    1549 
    1550 !  dimensions pour les champs scalaires et le vent zonal
    1551      start(1)=1
    1552      start(2)=1
    1553      start(3)=1
    1554      start(4)=timestep
    1555 
    1556      count(1)=1
    1557      count(2)=jjp1
    1558      count(3)=nlevnc
    1559      count(4)=1
    1560 
    1561 !  Pression
    1562      IF (guide_plevs==2) THEN
    1563          status=nf90_get_var(ncidp,varidp,zu,start,count)
    1564          DO i=1,iip1
    1565              pnat2(i,:,:)=zu(:,:)
    1566          ENDDO
    1567 
    1568          IF (invert_y) THEN
    1569 !           PRINT*,"Invertion impossible actuellement"
    1570 !           CALL abort_gcm(modname,abort_message,1)
    1571            CALL invert_lat(iip1,jjp1,nlevnc,pnat2)
    1572          ENDIF
    1573      endif
    1574 !  Vent zonal
    1575      IF (guide_u) THEN
    1576          status=nf90_get_var(ncidu,varidu,zu,start,count)
    1577          DO i=1,iip1
    1578              unat2(i,:,:)=zu(:,:)
    1579          ENDDO
    1580 
    1581          IF (invert_y) THEN
    1582            CALL invert_lat(iip1,jjp1,nlevnc,unat2)
    1583          ENDIF
    1584 
    1585      endif
    1586 
    1587 !  Temperature
    1588      IF (guide_T) THEN
    1589          status=nf90_get_var(ncidt,varidt,zu,start,count)
    1590          DO i=1,iip1
    1591              tnat2(i,:,:)=zu(:,:)
    1592          ENDDO
    1593 
    1594          IF (invert_y) THEN
    1595            CALL invert_lat(iip1,jjp1,nlevnc,tnat2)
    1596          ENDIF
    1597 
    1598      endif
    1599 
    1600 !  Humidite
    1601      IF (guide_Q) THEN
    1602          status=nf90_get_var(ncidQ,varidQ,zu,start,count)
    1603          DO i=1,iip1
    1604              qnat2(i,:,:)=zu(:,:)
    1605          ENDDO
    1606 
    1607          IF (invert_y) THEN
    1608            CALL invert_lat(iip1,jjp1,nlevnc,qnat2)
    1609          ENDIF
    1610 
    1611      endif
    1612 
    1613 !  Vent meridien
    1614      IF (guide_v) THEN
    1615          count(2)=jjm
    1616          status=nf90_get_var(ncidv,varidv,zv,start,count)
    1617          DO i=1,iip1
    1618              vnat2(i,:,:)=zv(:,:)
    1619          ENDDO
    1620 
    1621          IF (invert_y) THEN
    1622            CALL invert_lat(iip1,jjm,nlevnc,vnat2)
    1623          ENDIF
    1624 
    1625      endif
    1626 
    1627 !  Pression de surface
    1628      IF ((guide_P).OR.(guide_plevs==1))  THEN
    1629          start(3)=timestep
    1630          start(4)=0
    1631          count(2)=jjp1
    1632          count(3)=1
    1633          count(4)=0
    1634          status=nf90_get_var(ncidps,varidps,zu(:,1),start,count)
    1635          DO i=1,iip1
    1636              psnat2(i,:)=zu(:,1)
    1637          ENDDO
    1638 
    1639          IF (invert_y) THEN
    1640            CALL invert_lat(iip1,jjp1,1,psnat2)
    1641          ENDIF
    1642 
    1643      endif
     1413      ncidpl = -99
     1414      WRITE(*, *)trim(modname) // ' : opening nudging files '
     1415      ! Ap et Bp si niveaux de pression hybrides
     1416      IF (guide_plevs==1) THEN
     1417        WRITE(*, *)trim(modname) // ' Reading nudging on model levels'
     1418        rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
     1419        IF (rcode/=nf90_noerr) THEN
     1420          abort_message = 'Nudging: error -> no file apbp.nc'
     1421          CALL abort_gcm(modname, abort_message, 1)
     1422        ENDIF
     1423        rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
     1424        IF (rcode/=nf90_noerr) THEN
     1425          abort_message = 'Nudging: error -> no AP variable in file apbp.nc'
     1426          CALL abort_gcm(modname, abort_message, 1)
     1427        ENDIF
     1428        rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
     1429        IF (rcode/=nf90_noerr) THEN
     1430          abort_message = 'Nudging: error -> no BP variable in file apbp.nc'
     1431          CALL abort_gcm(modname, abort_message, 1)
     1432        ENDIF
     1433        WRITE(*, *)trim(modname) // 'ncidpl,varidap', ncidpl, varidap
     1434      endif
     1435      ! Pression
     1436      IF (guide_plevs==2) THEN
     1437        rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
     1438        IF (rcode/=nf90_noerr) THEN
     1439          abort_message = 'Nudging: error -> no file P.nc'
     1440          CALL abort_gcm(modname, abort_message, 1)
     1441        ENDIF
     1442        rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
     1443        IF (rcode/=nf90_noerr) THEN
     1444          abort_message = 'Nudging: error -> no PRES variable in file P.nc'
     1445          CALL abort_gcm(modname, abort_message, 1)
     1446        ENDIF
     1447        WRITE(*, *)trim(modname) // ' ncidp,varidp', ncidp, varidp
     1448        IF (ncidpl==-99) ncidpl = ncidp
     1449      endif
     1450      ! Vent zonal
     1451      IF (guide_u) THEN
     1452        rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
     1453        IF (rcode/=nf90_noerr) THEN
     1454          abort_message = 'Nudging: error -> no file u.nc'
     1455          CALL abort_gcm(modname, abort_message, 1)
     1456        ENDIF
     1457        rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
     1458        IF (rcode/=nf90_noerr) THEN
     1459          abort_message = 'Nudging: error -> no UWND variable in file u.nc'
     1460          CALL abort_gcm(modname, abort_message, 1)
     1461        ENDIF
     1462        WRITE(*, *)trim(modname) // ' ncidu,varidu', ncidu, varidu
     1463        IF (ncidpl==-99) ncidpl = ncidu
     1464      endif
     1465      ! Vent meridien
     1466      IF (guide_v) THEN
     1467        rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
     1468        IF (rcode/=nf90_noerr) THEN
     1469          abort_message = 'Nudging: error -> no file v.nc'
     1470          CALL abort_gcm(modname, abort_message, 1)
     1471        ENDIF
     1472        rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
     1473        IF (rcode/=nf90_noerr) THEN
     1474          abort_message = 'Nudging: error -> no VWND variable in file v.nc'
     1475          CALL abort_gcm(modname, abort_message, 1)
     1476        ENDIF
     1477        WRITE(*, *)trim(modname) // ' ncidv,varidv', ncidv, varidv
     1478        IF (ncidpl==-99) ncidpl = ncidv
     1479      endif
     1480      ! Temperature
     1481      IF (guide_T) THEN
     1482        rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
     1483        IF (rcode/=nf90_noerr) THEN
     1484          abort_message = 'Nudging: error -> no file T.nc'
     1485          CALL abort_gcm(modname, abort_message, 1)
     1486        ENDIF
     1487        rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
     1488        IF (rcode/=nf90_noerr) THEN
     1489          abort_message = 'Nudging: error -> no AIR variable in file T.nc'
     1490          CALL abort_gcm(modname, abort_message, 1)
     1491        ENDIF
     1492        WRITE(*, *)trim(modname) // ' ncidT,varidT', ncidt, varidt
     1493        IF (ncidpl==-99) ncidpl = ncidt
     1494      endif
     1495      ! Humidite
     1496      IF (guide_Q) THEN
     1497        rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
     1498        IF (rcode/=nf90_noerr) THEN
     1499          abort_message = 'Nudging: error -> no file hur.nc'
     1500          CALL abort_gcm(modname, abort_message, 1)
     1501        ENDIF
     1502        rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
     1503        IF (rcode/=nf90_noerr) THEN
     1504          abort_message = 'Nudging: error -> no RH,variable in file hur.nc'
     1505          CALL abort_gcm(modname, abort_message, 1)
     1506        ENDIF
     1507        WRITE(*, *)trim(modname) // ' ncidQ,varidQ', ncidQ, varidQ
     1508        IF (ncidpl==-99) ncidpl = ncidQ
     1509      endif
     1510      ! Pression de surface
     1511      IF ((guide_P).OR.(guide_modele)) THEN
     1512        rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
     1513        IF (rcode/=nf90_noerr) THEN
     1514          abort_message = 'Nudging: error -> no file ps.nc'
     1515          CALL abort_gcm(modname, abort_message, 1)
     1516        ENDIF
     1517        rcode = nf90_inq_varid(ncidps, 'SP', varidps)
     1518        IF (rcode/=nf90_noerr) THEN
     1519          abort_message = 'Nudging: error -> no SP variable in file ps.nc'
     1520          CALL abort_gcm(modname, abort_message, 1)
     1521        ENDIF
     1522        WRITE(*, *)trim(modname) // ' ncidps,varidps', ncidps, varidps
     1523      endif
     1524      ! Coordonnee verticale
     1525      IF (guide_plevs==0) THEN
     1526        rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
     1527        IF (rcode/=0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
     1528        WRITE(*, *)trim(modname) // ' ncidpl,varidpl', ncidpl, varidpl
     1529      endif
     1530      ! Coefs ap, bp pour calcul de la pression aux differents niveaux
     1531      IF (guide_plevs==1) THEN
     1532        status = nf90_get_var(ncidpl, varidap, apnc, [1], [nlevnc])
     1533        status = nf90_get_var(ncidpl, varidbp, bpnc, [1], [nlevnc])
     1534      elseif (guide_plevs==0) THEN
     1535        status = nf90_get_var(ncidpl, varidpl, apnc, [1], [nlevnc])
     1536        apnc = apnc * 100.! conversion en Pascals
     1537        bpnc(:) = 0.
     1538      endif
     1539      first = .FALSE.
     1540    endif ! (first)
     1541
     1542    ! -----------------------------------------------------------------
     1543    !   lecture des champs u, v, T, Q, ps
     1544    ! -----------------------------------------------------------------
     1545
     1546    !  dimensions pour les champs scalaires et le vent zonal
     1547    start(1) = 1
     1548    start(2) = 1
     1549    start(3) = 1
     1550    start(4) = timestep
     1551
     1552    count(1) = 1
     1553    count(2) = jjp1
     1554    count(3) = nlevnc
     1555    count(4) = 1
     1556
     1557    !  Pression
     1558    IF (guide_plevs==2) THEN
     1559      status = nf90_get_var(ncidp, varidp, zu, start, count)
     1560      DO i = 1, iip1
     1561        pnat2(i, :, :) = zu(:, :)
     1562      ENDDO
     1563
     1564      IF (invert_y) THEN
     1565        !           PRINT*,"Invertion impossible actuellement"
     1566        !           CALL abort_gcm(modname,abort_message,1)
     1567        CALL invert_lat(iip1, jjp1, nlevnc, pnat2)
     1568      ENDIF
     1569    endif
     1570    !  Vent zonal
     1571    IF (guide_u) THEN
     1572      status = nf90_get_var(ncidu, varidu, zu, start, count)
     1573      DO i = 1, iip1
     1574        unat2(i, :, :) = zu(:, :)
     1575      ENDDO
     1576
     1577      IF (invert_y) THEN
     1578        CALL invert_lat(iip1, jjp1, nlevnc, unat2)
     1579      ENDIF
     1580
     1581    endif
     1582
     1583    !  Temperature
     1584    IF (guide_T) THEN
     1585      status = nf90_get_var(ncidt, varidt, zu, start, count)
     1586      DO i = 1, iip1
     1587        tnat2(i, :, :) = zu(:, :)
     1588      ENDDO
     1589
     1590      IF (invert_y) THEN
     1591        CALL invert_lat(iip1, jjp1, nlevnc, tnat2)
     1592      ENDIF
     1593
     1594    endif
     1595
     1596    !  Humidite
     1597    IF (guide_Q) THEN
     1598      status = nf90_get_var(ncidQ, varidQ, zu, start, count)
     1599      DO i = 1, iip1
     1600        qnat2(i, :, :) = zu(:, :)
     1601      ENDDO
     1602
     1603      IF (invert_y) THEN
     1604        CALL invert_lat(iip1, jjp1, nlevnc, qnat2)
     1605      ENDIF
     1606
     1607    endif
     1608
     1609    !  Vent meridien
     1610    IF (guide_v) THEN
     1611      count(2) = jjm
     1612      status = nf90_get_var(ncidv, varidv, zv, start, count)
     1613      DO i = 1, iip1
     1614        vnat2(i, :, :) = zv(:, :)
     1615      ENDDO
     1616
     1617      IF (invert_y) THEN
     1618        CALL invert_lat(iip1, jjm, nlevnc, vnat2)
     1619      ENDIF
     1620
     1621    endif
     1622
     1623    !  Pression de surface
     1624    IF ((guide_P).OR.(guide_plevs==1))  THEN
     1625      start(3) = timestep
     1626      start(4) = 0
     1627      count(2) = jjp1
     1628      count(3) = 1
     1629      count(4) = 0
     1630      status = nf90_get_var(ncidps, varidps, zu(:, 1), start, count)
     1631      DO i = 1, iip1
     1632        psnat2(i, :) = zu(:, 1)
     1633      ENDDO
     1634
     1635      IF (invert_y) THEN
     1636        CALL invert_lat(iip1, jjp1, 1, psnat2)
     1637      ENDIF
     1638
     1639    endif
    16441640
    16451641  END SUBROUTINE guide_read2D
    1646  
    1647 !=======================================================================
    1648   SUBROUTINE guide_out(varname,hsize,vsize,field)
     1642
     1643  !=======================================================================
     1644  SUBROUTINE guide_out(varname, hsize, vsize, field)
    16491645
    16501646    USE comconst_mod, ONLY: pi
    … …  
    16571653    INCLUDE "paramet.h"
    16581654    INCLUDE "comgeom2.h"
    1659    
     1655
    16601656    ! Variables entree
    1661     CHARACTER*(*), INTENT(IN)                          :: varname
    1662     INTEGER,   INTENT (IN)                         :: hsize,vsize
    1663     REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field
     1657    CHARACTER*(*), INTENT(IN) :: varname
     1658    INTEGER, INTENT (IN) :: hsize, vsize
     1659    REAL, DIMENSION (iip1, hsize, vsize), INTENT(IN) :: field
    16641660
    16651661    ! Variables locales
    1666     INTEGER, SAVE :: timestep=0
     1662    INTEGER, SAVE :: timestep = 0
    16671663    ! Identites fichier netcdf
    1668     INTEGER       :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev
    1669     INTEGER       :: vid_lonu,vid_lonv,vid_latu,vid_latv,vid_cu,vid_cv,vid_lev
    1670     INTEGER       :: vid_au,vid_av, varid_alpha_t, varid_alpha_q
     1664    INTEGER :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev
     1665    INTEGER :: vid_lonu, vid_lonv, vid_latu, vid_latv, vid_cu, vid_cv, vid_lev
     1666    INTEGER :: vid_au, vid_av, varid_alpha_t, varid_alpha_q
    16711667    INTEGER, DIMENSION (3) :: dim3
    1672     INTEGER, DIMENSION (4) :: dim4,count,start
    1673     INTEGER                :: ierr, varid,l
    1674     REAL, DIMENSION (iip1,hsize,vsize) :: field2
    1675     CHARACTER(LEN=20),PARAMETER :: modname="guide_out"
    1676 
    1677     WRITE(*,*)trim(modname)//': output timestep',timestep,'var ',varname
     1668    INTEGER, DIMENSION (4) :: dim4, count, start
     1669    INTEGER :: ierr, varid, l
     1670    REAL, DIMENSION (iip1, hsize, vsize) :: field2
     1671    CHARACTER(LEN = 20), PARAMETER :: modname = "guide_out"
     1672
     1673    WRITE(*, *)trim(modname) // ': output timestep', timestep, 'var ', varname
    16781674    IF (timestep==0) THEN
    1679 ! ----------------------------------------------
    1680 ! initialisation fichier de sortie
    1681 ! ----------------------------------------------
    1682 ! Ouverture du fichier
    1683         ierr=nf90_create("guide_ins.nc",IOR(nf90_clobber,nf90_64bit_offset),nid)
    1684 ! Definition des dimensions
    1685         ierr=nf90_def_dim(nid,"LONU",iip1,id_lonu)
    1686         ierr=nf90_def_dim(nid,"LONV",iip1,id_lonv)
    1687         ierr=nf90_def_dim(nid,"LATU",jjp1,id_latu)
    1688         ierr=nf90_def_dim(nid,"LATV",jjm,id_latv)
    1689         ierr=nf90_def_dim(nid,"LEVEL",llm,id_lev)
    1690         ierr=nf90_def_dim(nid,"TIME",nf90_unlimited,id_tim)
    1691 
    1692 ! Creation des variables dimensions
    1693         ierr=nf90_def_var(nid,"LONU",nf90_float,id_lonu,vid_lonu)
    1694         ierr=nf90_def_var(nid,"LONV",nf90_float,id_lonv,vid_lonv)
    1695         ierr=nf90_def_var(nid,"LATU",nf90_float,id_latu,vid_latu)
    1696         ierr=nf90_def_var(nid,"LATV",nf90_float,id_latv,vid_latv)
    1697         ierr=nf90_def_var(nid,"LEVEL",nf90_float,id_lev,vid_lev)
    1698         ierr=nf90_def_var(nid,"cu",nf90_float,(/id_lonu,id_latu/),vid_cu)
    1699         ierr=nf90_def_var(nid,"cv",nf90_float,(/id_lonv,id_latv/),vid_cv)
    1700         ierr=nf90_def_var(nid,"au",nf90_float,(/id_lonu,id_latu/),vid_au)
    1701         ierr=nf90_def_var(nid,"av",nf90_float,(/id_lonv,id_latv/),vid_av)
    1702         CALL nf95_def_var(nid, "alpha_T", nf90_float, (/id_lonv, id_latu/), &
    1703              varid_alpha_t)
    1704         CALL nf95_def_var(nid, "alpha_q", nf90_float, (/id_lonv, id_latu/), &
    1705              varid_alpha_q)
    1706        
    1707         ierr=nf90_enddef(nid)
    1708 
    1709 ! Enregistrement des variables dimensions
    1710         ierr = nf90_put_var(nid,vid_lonu,rlonu*180./pi)
    1711         ierr = nf90_put_var(nid,vid_lonv,rlonv*180./pi)
    1712         ierr = nf90_put_var(nid,vid_latu,rlatu*180./pi)
    1713         ierr = nf90_put_var(nid,vid_latv,rlatv*180./pi)
    1714         ierr = nf90_put_var(nid,vid_lev,presnivs)
    1715         ierr = nf90_put_var(nid,vid_cu,cu)
    1716         ierr = nf90_put_var(nid,vid_cv,cv)
    1717         ierr = nf90_put_var(nid,vid_au,alpha_u)
    1718         ierr = nf90_put_var(nid,vid_av,alpha_v)
    1719         CALL nf95_put_var(nid, varid_alpha_t, alpha_t)
    1720         CALL nf95_put_var(nid, varid_alpha_q, alpha_q)
    1721 ! --------------------------------------------------------------------
    1722 ! Création des variables sauvegardées
    1723 ! --------------------------------------------------------------------
    1724         ierr = nf90_redef(nid)
    1725 ! Pressure (GCM)
    1726         dim4=(/id_lonv,id_latu,id_lev,id_tim/)
    1727         ierr = nf90_def_var(nid,"SP",nf90_float,dim4,varid)
    1728 ! Surface pressure (guidage)
    1729         IF (guide_P) THEN
    1730             dim3=(/id_lonv,id_latu,id_tim/)
    1731             ierr = nf90_def_var(nid,"ps",nf90_float,dim3,varid)
    1732         ENDIF
    1733 ! Zonal wind
    1734         IF (guide_u) THEN
    1735             dim4=(/id_lonu,id_latu,id_lev,id_tim/)
    1736             ierr = nf90_def_var(nid,"u",nf90_float,dim4,varid)
    1737             ierr = nf90_def_var(nid,"ua",nf90_float,dim4,varid)
    1738             ierr = nf90_def_var(nid,"ucov",nf90_float,dim4,varid)
    1739         ENDIF
    1740 ! Merid. wind
    1741         IF (guide_v) THEN
    1742             dim4=(/id_lonv,id_latv,id_lev,id_tim/)
    1743             ierr = nf90_def_var(nid,"v",nf90_float,dim4,varid)
    1744             ierr = nf90_def_var(nid,"va",nf90_float,dim4,varid)
    1745             ierr = nf90_def_var(nid,"vcov",nf90_float,dim4,varid)
    1746         ENDIF
    1747 ! Pot. Temperature
    1748         IF (guide_T) THEN
    1749             dim4=(/id_lonv,id_latu,id_lev,id_tim/)
    1750             ierr = nf90_def_var(nid,"teta",nf90_float,dim4,varid)
    1751         ENDIF
    1752 ! Specific Humidity
    1753         IF (guide_Q) THEN
    1754             dim4=(/id_lonv,id_latu,id_lev,id_tim/)
    1755             ierr = nf90_def_var(nid,"q",nf90_float,dim4,varid)
    1756         ENDIF
    1757        
    1758         ierr = nf90_enddef(nid)
    1759         ierr = nf90_close(nid)
     1675      ! ----------------------------------------------
     1676      ! initialisation fichier de sortie
     1677      ! ----------------------------------------------
     1678      ! Ouverture du fichier
     1679      ierr = nf90_create("guide_ins.nc", IOR(nf90_clobber, nf90_64bit_offset), nid)
     1680      ! Definition des dimensions
     1681      ierr = nf90_def_dim(nid, "LONU", iip1, id_lonu)
     1682      ierr = nf90_def_dim(nid, "LONV", iip1, id_lonv)
     1683      ierr = nf90_def_dim(nid, "LATU", jjp1, id_latu)
     1684      ierr = nf90_def_dim(nid, "LATV", jjm, id_latv)
     1685      ierr = nf90_def_dim(nid, "LEVEL", llm, id_lev)
     1686      ierr = nf90_def_dim(nid, "TIME", nf90_unlimited, id_tim)
     1687
     1688      ! Creation des variables dimensions
     1689      ierr = nf90_def_var(nid, "LONU", nf90_float, id_lonu, vid_lonu)
     1690      ierr = nf90_def_var(nid, "LONV", nf90_float, id_lonv, vid_lonv)
     1691      ierr = nf90_def_var(nid, "LATU", nf90_float, id_latu, vid_latu)
     1692      ierr = nf90_def_var(nid, "LATV", nf90_float, id_latv, vid_latv)
     1693      ierr = nf90_def_var(nid, "LEVEL", nf90_float, id_lev, vid_lev)
     1694      ierr = nf90_def_var(nid, "cu", nf90_float, (/id_lonu, id_latu/), vid_cu)
     1695      ierr = nf90_def_var(nid, "cv", nf90_float, (/id_lonv, id_latv/), vid_cv)
     1696      ierr = nf90_def_var(nid, "au", nf90_float, (/id_lonu, id_latu/), vid_au)
     1697      ierr = nf90_def_var(nid, "av", nf90_float, (/id_lonv, id_latv/), vid_av)
     1698      CALL nf95_def_var(nid, "alpha_T", nf90_float, (/id_lonv, id_latu/), &
     1699              varid_alpha_t)
     1700      CALL nf95_def_var(nid, "alpha_q", nf90_float, (/id_lonv, id_latu/), &
     1701              varid_alpha_q)
     1702
     1703      ierr = nf90_enddef(nid)
     1704
     1705      ! Enregistrement des variables dimensions
     1706      ierr = nf90_put_var(nid, vid_lonu, rlonu * 180. / pi)
     1707      ierr = nf90_put_var(nid, vid_lonv, rlonv * 180. / pi)
     1708      ierr = nf90_put_var(nid, vid_latu, rlatu * 180. / pi)
     1709      ierr = nf90_put_var(nid, vid_latv, rlatv * 180. / pi)
     1710      ierr = nf90_put_var(nid, vid_lev, presnivs)
     1711      ierr = nf90_put_var(nid, vid_cu, cu)
     1712      ierr = nf90_put_var(nid, vid_cv, cv)
     1713      ierr = nf90_put_var(nid, vid_au, alpha_u)
     1714      ierr = nf90_put_var(nid, vid_av, alpha_v)
     1715      CALL nf95_put_var(nid, varid_alpha_t, alpha_t)
     1716      CALL nf95_put_var(nid, varid_alpha_q, alpha_q)
     1717      ! --------------------------------------------------------------------
     1718      ! Création des variables sauvegardées
     1719      ! --------------------------------------------------------------------
     1720      ierr = nf90_redef(nid)
     1721      ! Pressure (GCM)
     1722      dim4 = (/id_lonv, id_latu, id_lev, id_tim/)
     1723      ierr = nf90_def_var(nid, "SP", nf90_float, dim4, varid)
     1724      ! Surface pressure (guidage)
     1725      IF (guide_P) THEN
     1726        dim3 = (/id_lonv, id_latu, id_tim/)
     1727        ierr = nf90_def_var(nid, "ps", nf90_float, dim3, varid)
     1728      ENDIF
     1729      ! Zonal wind
     1730      IF (guide_u) THEN
     1731        dim4 = (/id_lonu, id_latu, id_lev, id_tim/)
     1732        ierr = nf90_def_var(nid, "u", nf90_float, dim4, varid)
     1733        ierr = nf90_def_var(nid, "ua", nf90_float, dim4, varid)
     1734        ierr = nf90_def_var(nid, "ucov", nf90_float, dim4, varid)
     1735      ENDIF
     1736      ! Merid. wind
     1737      IF (guide_v) THEN
     1738        dim4 = (/id_lonv, id_latv, id_lev, id_tim/)
     1739        ierr = nf90_def_var(nid, "v", nf90_float, dim4, varid)
     1740        ierr = nf90_def_var(nid, "va", nf90_float, dim4, varid)
     1741        ierr = nf90_def_var(nid, "vcov", nf90_float, dim4, varid)
     1742      ENDIF
     1743      ! Pot. Temperature
     1744      IF (guide_T) THEN
     1745        dim4 = (/id_lonv, id_latu, id_lev, id_tim/)
     1746        ierr = nf90_def_var(nid, "teta", nf90_float, dim4, varid)
     1747      ENDIF
     1748      ! Specific Humidity
     1749      IF (guide_Q) THEN
     1750        dim4 = (/id_lonv, id_latu, id_lev, id_tim/)
     1751        ierr = nf90_def_var(nid, "q", nf90_float, dim4, varid)
     1752      ENDIF
     1753
     1754      ierr = nf90_enddef(nid)
     1755      ierr = nf90_close(nid)
    17601756    ENDIF ! timestep=0
    17611757
    1762 ! --------------------------------------------------------------------
    1763 ! Enregistrement du champ
    1764 ! --------------------------------------------------------------------
    1765     ierr=nf90_open("guide_ins.nc",nf90_write,nid)
    1766 
    1767     IF (varname=="SP") timestep=timestep+1
    1768 
    1769     ierr = nf90_inq_varid(nid,varname,varid)
     1758    ! --------------------------------------------------------------------
     1759    ! Enregistrement du champ
     1760    ! --------------------------------------------------------------------
     1761    ierr = nf90_open("guide_ins.nc", nf90_write, nid)
     1762
     1763    IF (varname=="SP") timestep = timestep + 1
     1764
     1765    ierr = nf90_inq_varid(nid, varname, varid)
    17701766    SELECT CASE (varname)
    1771     CASE ("SP","ps")
    1772         start=(/1,1,1,timestep/)
    1773         count=(/iip1,jjp1,llm,1/)
    1774     CASE ("v","va","vcov")
    1775         start=(/1,1,1,timestep/)
    1776         count=(/iip1,jjm,llm,1/)
     1767    CASE ("SP", "ps")
     1768      start = (/1, 1, 1, timestep/)
     1769      count = (/iip1, jjp1, llm, 1/)
     1770    CASE ("v", "va", "vcov")
     1771      start = (/1, 1, 1, timestep/)
     1772      count = (/iip1, jjm, llm, 1/)
    17771773    CASE DEFAULT
    1778         start=(/1,1,1,timestep/)
    1779         count=(/iip1,jjp1,llm,1/)
     1774      start = (/1, 1, 1, timestep/)
     1775      count = (/iip1, jjp1, llm, 1/)
    17801776    END SELECT
    17811777
    17821778    SELECT CASE (varname)
    1783     CASE("u","ua")
    1784         DO l=1,llm ; field2(:,2:jjm,l)=field(:,2:jjm,l)/cu(:,2:jjm) ; ENDDO
    1785         field2(:,1,:)=0. ; field2(:,jjp1,:)=0.
    1786     CASE("v","va")
    1787         DO l=1,llm ; field2(:,:,l)=field(:,:,l)/cv(:,:) ; ENDDO
     1779    CASE("u", "ua")
     1780      DO l = 1, llm ; field2(:, 2:jjm, l) = field(:, 2:jjm, l) / cu(:, 2:jjm) ;
     1781      ENDDO
     1782      field2(:, 1, :) = 0. ; field2(:, jjp1, :) = 0.
     1783    CASE("v", "va")
     1784      DO l = 1, llm ; field2(:, :, l) = field(:, :, l) / cv(:, :) ;
     1785      ENDDO
    17881786    CASE DEFAULT
    1789         field2=field
     1787      field2 = field
    17901788    END SELECT
    17911789
    1792     ierr = nf90_put_var(nid,varid,field2,start,count)
     1790    ierr = nf90_put_var(nid, varid, field2, start, count)
    17931791    ierr = nf90_close(nid)
    17941792
    17951793  END SUBROUTINE guide_out
    1796    
    1797  
    1798 !===========================================================================
    1799   SUBROUTINE correctbid(iim,nl,x)
    1800     INTEGER iim,nl
    1801     REAL x(iim+1,nl)
    1802     INTEGER i,l
     1794
     1795
     1796  !===========================================================================
     1797  SUBROUTINE correctbid(iim, nl, x)
     1798    INTEGER iim, nl
     1799    REAL x(iim + 1, nl)
     1800    INTEGER i, l
    18031801    REAL zz
    18041802
    1805     do l=1,nl
    1806         do i=2,iim-1
    1807             IF(abs(x(i,l))>1.e10) THEN
    1808                zz=0.5*(x(i-1,l)+x(i+1,l))
    1809               PRINT*,'correction ',i,l,x(i,l),zz
    1810                x(i,l)=zz
    1811             endif
    1812          enddo
    1813      enddo
     1803    do l = 1, nl
     1804      do i = 2, iim - 1
     1805        IF(abs(x(i, l))>1.e10) THEN
     1806          zz = 0.5 * (x(i - 1, l) + x(i + 1, l))
     1807          PRINT*, 'correction ', i, l, x(i, l), zz
     1808          x(i, l) = zz
     1809        endif
     1810      enddo
     1811    enddo
    18141812
    18151813  END SUBROUTINE  correctbid
    18161814
    1817 !===========================================================================
     1815  !===========================================================================
    18181816END MODULE guide_mod

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/iniacademic.F90

    r5117 r5118  
    1 
    21! $Id$
    32
    4 SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
     3SUBROUTINE iniacademic(vcov, ucov, teta, q, masse, ps, phis, time_0)
    54
    65  USE lmdz_filtreg, ONLY: inifilr
    7   USE infotrac,    ONLY: nqtot, niso, iqIsoPha, tracers, getKey, isoName
    8   USE control_mod, ONLY: day_step,planet_type
     6  USE infotrac, ONLY: nqtot, niso, iqIsoPha, tracers, getKey, isoName
     7  USE control_mod, ONLY: day_step, planet_type
    98  USE exner_hyb_m, ONLY: exner_hyb
    109  USE exner_milieu_m, ONLY: exner_milieu
    … …  
    1514  USE comvert_mod, ONLY: ap, bp, preff, pa, presnivs, pressure_exner
    1615  USE temps_mod, ONLY: annee_ref, day_ini, day_ref
    17   USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
     16  USE ener_mod, ONLY: etot0, ptot0, ztot0, stot0, ang0
    1817  USE lmdz_readTracFiles, ONLY: addPhase
    19   USE netcdf, ONLY: nf90_nowrite,nf90_open,nf90_noerr,nf90_inq_varid,nf90_close,nf90_get_var
     18  USE netcdf, ONLY: nf90_nowrite, nf90_open, nf90_noerr, nf90_inq_varid, nf90_close, nf90_get_var
    2019  USE lmdz_ran1, ONLY: ran1
     20  USE lmdz_iniprint, ONLY: lunout, prt_level
    2121
    2222  !   Author:    Frederic Hourdin      original: 15/01/93
    … …  
    3333  include "comgeom.h"
    3434  include "academic.h"
    35   include "iniprint.h"
    3635
    3736  !   Arguments:
    3837  !   ----------
    3938
    40   REAL,INTENT(OUT) :: time_0
     39  REAL, INTENT(OUT) :: time_0
    4140
    4241  !   fields
    43   REAL,INTENT(OUT) :: vcov(ip1jm,llm) ! meridional covariant wind
    44   REAL,INTENT(OUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind
    45   REAL,INTENT(OUT) :: teta(ip1jmp1,llm) ! potential temperature (K)
    46   REAL,INTENT(OUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers (.../kg_of_air)
    47   REAL,INTENT(OUT) :: ps(ip1jmp1) ! surface pressure (Pa)
    48   REAL,INTENT(OUT) :: masse(ip1jmp1,llm) ! air mass in grid cell (kg)
    49   REAL,INTENT(OUT) :: phis(ip1jmp1) ! surface geopotential
     42  REAL, INTENT(OUT) :: vcov(ip1jm, llm) ! meridional covariant wind
     43  REAL, INTENT(OUT) :: ucov(ip1jmp1, llm) ! zonal covariant wind
     44  REAL, INTENT(OUT) :: teta(ip1jmp1, llm) ! potential temperature (K)
     45  REAL, INTENT(OUT) :: q(ip1jmp1, llm, nqtot) ! advected tracers (.../kg_of_air)
     46  REAL, INTENT(OUT) :: ps(ip1jmp1) ! surface pressure (Pa)
     47  REAL, INTENT(OUT) :: masse(ip1jmp1, llm) ! air mass in grid cell (kg)
     48  REAL, INTENT(OUT) :: phis(ip1jmp1) ! surface geopotential
    5049
    5150  !   Local:
    5251  !   ------
    5352
    54   REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
     53  REAL p (ip1jmp1, llmp1)               ! pression aux interfac.des couches
    5554  REAL pks(ip1jmp1)                      ! exner au  sol
    56   REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
    57   REAL phi(ip1jmp1,llm)                  ! geopotentiel
    58   REAL ddsin,zsig,tetapv,w_pv  ! variables auxiliaires
     55  REAL pk(ip1jmp1, llm)                   ! exner au milieu des couches
     56  REAL phi(ip1jmp1, llm)                  ! geopotentiel
     57  REAL ddsin, zsig, tetapv, w_pv  ! variables auxiliaires
    5958  REAL tetastrat ! potential temperature in the stratosphere, in K
    60   REAL tetajl(jjp1,llm)
    61   INTEGER i,j,l,lsup,ij, iq, iName, iPhase, iqParent
    62 
    63   INTEGER :: nid_relief,varid,ierr
    64   REAL, DIMENSION(iip1,jjp1) :: relief
    65 
    66   REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
    67   REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
     59  REAL tetajl(jjp1, llm)
     60  INTEGER i, j, l, lsup, ij, iq, iName, iPhase, iqParent
     61
     62  INTEGER :: nid_relief, varid, ierr
     63  REAL, DIMENSION(iip1, jjp1) :: relief
     64
     65  REAL teta0, ttp, delt_y, delt_z, eps ! Constantes pour profil de T
     66  REAL k_f, k_c_a, k_c_s         ! Constantes de rappel
    6867  LOGICAL ok_geost             ! Initialisation vent geost. ou nul
    6968  LOGICAL ok_pv                ! Polar Vortex
    70   REAL phi_pv,dphi_pv,gam_pv,tetanoise   ! Constantes pour polar vortex
     69  REAL phi_pv, dphi_pv, gam_pv, tetanoise   ! Constantes pour polar vortex
    7170
    7271  REAL zz
    … …  
    7473
    7574  REAL zdtvr, tnat, alpha_ideal
    76   LOGICAL,PARAMETER :: tnat1=.TRUE.
    77  
    78   CHARACTER(LEN=*),parameter :: modname="iniacademic"
    79   CHARACTER(LEN=80) :: abort_message
     75  LOGICAL, PARAMETER :: tnat1 = .TRUE.
     76
     77  CHARACTER(LEN = *), parameter :: modname = "iniacademic"
     78  CHARACTER(LEN = 80) :: abort_message
    8079
    8180  ! Sanity check: verify that options selected by user are not incompatible
    8281  IF ((iflag_phys==1).AND. .NOT. read_start) THEN
    83     WRITE(lunout,*) trim(modname)," error: if read_start is set to ", &
    84     " false then iflag_phys should not be 1"
    85     WRITE(lunout,*) "You most likely want an aquaplanet initialisation", &
    86     " (iflag_phys >= 100)"
    87     CALL abort_gcm(modname,"incompatible iflag_phys==1 and read_start==.FALSE.",1)
     82    WRITE(lunout, *) trim(modname), " error: if read_start is set to ", &
     83            " false then iflag_phys should not be 1"
     84    WRITE(lunout, *) "You most likely want an aquaplanet initialisation", &
     85            " (iflag_phys >= 100)"
     86    CALL abort_gcm(modname, "incompatible iflag_phys==1 and read_start==.FALSE.", 1)
    8887  ENDIF
    89  
     88
    9089  !-----------------------------------------------------------------------
    9190  ! 1. Initializations for Earth-like case
    … …  
    9594  CALL conf_planete
    9695
    97   time_0=0.
    98   day_ref=1
    99   annee_ref=0
    100 
    101   im         = iim
    102   jm         = jjm
    103   day_ini    = 1
    104   dtvr    = daysec/REAL(day_step)
    105   zdtvr=dtvr
    106   etot0      = 0.
    107   ptot0      = 0.
    108   ztot0      = 0.
    109   stot0      = 0.
    110   ang0       = 0.
     96  time_0 = 0.
     97  day_ref = 1
     98  annee_ref = 0
     99
     100  im = iim
     101  jm = jjm
     102  day_ini = 1
     103  dtvr = daysec / REAL(day_step)
     104  zdtvr = dtvr
     105  etot0 = 0.
     106  ptot0 = 0.
     107  ztot0 = 0.
     108  stot0 = 0.
     109  ang0 = 0.
    111110
    112111  IF (llm == 1) THEN
    113      ! specific initializations for the shallow water case
    114      kappa=1
     112    ! specific initializations for the shallow water case
     113    kappa = 1
    115114  ENDIF
    116115
    … …  
    126125  IF (.NOT. read_start) THEN
    127126
    128      !------------------------------------------------------------------
    129      ! Lecture eventuelle d'un fichier de relief interpollee sur la grille
    130      ! du modele.
    131      ! On suppose que le fichier relief_in.nc est stoké sur une grille
    132      ! iim*jjp1
    133      ! Facile a créer à partir de la commande
    134      ! cdo remapcon,fichier_output_phys.nc Relief.nc relief_in.nc
    135      !------------------------------------------------------------------
    136 
    137      relief=0.
    138      ierr = nf90_open ('relief_in.nc', nf90_nowrite,nid_relief)
    139      IF (ierr==nf90_noerr) THEN
    140          ierr=nf90_inq_varid(nid_relief,'RELIEF',varid)
    141          IF (ierr==nf90_noerr) THEN
    142               ierr=nf90_get_var(nid_relief,varid,relief(1:iim,1:jjp1))
    143               relief(iip1,:)=relief(1,:)
    144          else
    145               CALL abort_gcm ('iniacademic','variable RELIEF pas la',1)
    146          endif
    147      endif
    148      ierr = nf90_close (nid_relief)
    149 
    150      !------------------------------------------------------------------
    151      ! Initialisation du geopotentiel au sol et de la pression
    152      !------------------------------------------------------------------
    153 
    154      PRINT*,'relief=',minval(relief),maxval(relief),'g=',g
    155      do j=1,jjp1
    156         do i=1,iip1
    157            phis((j-1)*iip1+i)=g*relief(i,j)
     127    !------------------------------------------------------------------
     128    ! Lecture eventuelle d'un fichier de relief interpollee sur la grille
     129    ! du modele.
     130    ! On suppose que le fichier relief_in.nc est stoké sur une grille
     131    ! iim*jjp1
     132    ! Facile a créer à partir de la commande
     133    ! cdo remapcon,fichier_output_phys.nc Relief.nc relief_in.nc
     134    !------------------------------------------------------------------
     135
     136    relief = 0.
     137    ierr = nf90_open ('relief_in.nc', nf90_nowrite, nid_relief)
     138    IF (ierr==nf90_noerr) THEN
     139      ierr = nf90_inq_varid(nid_relief, 'RELIEF', varid)
     140      IF (ierr==nf90_noerr) THEN
     141        ierr = nf90_get_var(nid_relief, varid, relief(1:iim, 1:jjp1))
     142        relief(iip1, :) = relief(1, :)
     143      else
     144        CALL abort_gcm ('iniacademic', 'variable RELIEF pas la', 1)
     145      endif
     146    endif
     147    ierr = nf90_close (nid_relief)
     148
     149    !------------------------------------------------------------------
     150    ! Initialisation du geopotentiel au sol et de la pression
     151    !------------------------------------------------------------------
     152
     153    PRINT*, 'relief=', minval(relief), maxval(relief), 'g=', g
     154    do j = 1, jjp1
     155      do i = 1, iip1
     156        phis((j - 1) * iip1 + i) = g * relief(i, j)
     157      enddo
     158    enddo
     159    PRINT*, 'phis=', minval(phis), maxval(phis), 'g=', g
     160
     161    ! ground geopotential
     162    !phis(:)=0.
     163    ps(:) = preff
     164    CALL pression (ip1jmp1, ap, bp, ps, p)
     165
     166    IF (pressure_exner) THEN
     167      CALL exner_hyb(ip1jmp1, ps, p, pks, pk)
     168    else
     169      CALL exner_milieu(ip1jmp1, ps, p, pks, pk)
     170    endif
     171    CALL massdair(p, masse)
     172  ENDIF
     173
     174  IF (llm == 1) THEN
     175    ! initialize fields for the shallow water case, if required
     176    IF (.NOT.read_start) THEN
     177      phis(:) = 0.
     178      q(:, :, :) = 0
     179      CALL sw_case_williamson91_6(vcov, ucov, teta, masse, ps)
     180    endif
     181  ENDIF
     182
     183  academic_case: if (iflag_phys >= 2) THEN
     184    ! initializations
     185
     186    ! 1. local parameters
     187    ! by convention, winter is in the southern hemisphere
     188    ! Geostrophic wind or no wind?
     189    ok_geost = .TRUE.
     190    CALL getin('ok_geost', ok_geost)
     191    ! Constants for Newtonian relaxation and friction
     192    k_f = 1.                !friction
     193    CALL getin('k_j', k_f)
     194    k_f = 1. / (daysec * k_f)
     195    k_c_s = 4.  !cooling surface
     196    CALL getin('k_c_s', k_c_s)
     197    k_c_s = 1. / (daysec * k_c_s)
     198    k_c_a = 40. !cooling free atm
     199    CALL getin('k_c_a', k_c_a)
     200    k_c_a = 1. / (daysec * k_c_a)
     201    ! Constants for Teta equilibrium profile
     202    teta0 = 315.     ! mean Teta (S.H. 315K)
     203    CALL getin('teta0', teta0)
     204    ttp = 200.       ! Tropopause temperature (S.H. 200K)
     205    CALL getin('ttp', ttp)
     206    eps = 0.         ! Deviation to N-S symmetry(~0-20K)
     207    CALL getin('eps', eps)
     208    delt_y = 60.     ! Merid Temp. Gradient (S.H. 60K)
     209    CALL getin('delt_y', delt_y)
     210    delt_z = 10.     ! Vertical Gradient (S.H. 10K)
     211    CALL getin('delt_z', delt_z)
     212    ! Polar vortex
     213    ok_pv = .FALSE.
     214    CALL getin('ok_pv', ok_pv)
     215    phi_pv = -50.            ! Latitude of edge of vortex
     216    CALL getin('phi_pv', phi_pv)
     217    phi_pv = phi_pv * pi / 180.
     218    dphi_pv = 5.             ! Width of the edge
     219    CALL getin('dphi_pv', dphi_pv)
     220    dphi_pv = dphi_pv * pi / 180.
     221    gam_pv = 4.              ! -dT/dz vortex (in K/km)
     222    CALL getin('gam_pv', gam_pv)
     223    tetanoise = 0.005
     224    CALL getin('tetanoise', tetanoise)
     225
     226    ! 2. Initialize fields towards which to relax
     227    ! Friction
     228    knewt_g = k_c_a
     229    DO l = 1, llm
     230      zsig = presnivs(l) / preff
     231      knewt_t(l) = (k_c_s - k_c_a) * MAX(0., (zsig - 0.7) / 0.3)
     232      kfrict(l) = k_f * MAX(0., (zsig - 0.7) / 0.3)
     233    ENDDO
     234    DO j = 1, jjp1
     235      clat4((j - 1) * iip1 + 1:j * iip1) = cos(rlatu(j))**4
     236    ENDDO
     237
     238    ! Potential temperature
     239    DO l = 1, llm
     240      zsig = presnivs(l) / preff
     241      tetastrat = ttp * zsig**(-kappa)
     242      tetapv = tetastrat
     243      IF ((ok_pv).AND.(zsig<0.1)) THEN
     244        tetapv = tetastrat * (zsig * 10.)**(kappa * cpp * gam_pv / 1000. / g)
     245      ENDIF
     246      DO j = 1, jjp1
     247        ! Troposphere
     248        ddsin = sin(rlatu(j))
     249        tetajl(j, l) = teta0 - delt_y * ddsin * ddsin + eps * ddsin &
     250                - delt_z * (1. - ddsin * ddsin) * log(zsig)
     251        IF (planet_type=="giant") THEN
     252          tetajl(j, l) = teta0 + (delt_y * &
     253                  ((sin(rlatu(j) * 3.14159 * eps + 0.0001))**2)   &
     254                  / ((rlatu(j) * 3.14159 * eps + 0.0001)**2))     &
     255                  - delt_z * log(zsig)
     256        endif
     257        ! Profil stratospherique isotherme (+vortex)
     258        w_pv = (1. - tanh((rlatu(j) - phi_pv) / dphi_pv)) / 2.
     259        tetastrat = tetastrat * (1. - w_pv) + tetapv * w_pv
     260        tetajl(j, l) = MAX(tetajl(j, l), tetastrat)
     261      ENDDO
     262    ENDDO
     263
     264    !          CALL writefield('theta_eq',tetajl)
     265
     266    do l = 1, llm
     267      do j = 1, jjp1
     268        do i = 1, iip1
     269          ij = (j - 1) * iip1 + i
     270          tetarappel(ij, l) = tetajl(j, l)
    158271        enddo
    159      enddo
    160      PRINT*,'phis=',minval(phis),maxval(phis),'g=',g
    161 
    162      ! ground geopotential
    163      !phis(:)=0.
    164      ps(:)=preff
    165      CALL pression ( ip1jmp1, ap, bp, ps, p       )
    166 
    167      IF (pressure_exner) THEN
    168        CALL exner_hyb( ip1jmp1, ps, p, pks, pk)
    169      else
    170        CALL exner_milieu(ip1jmp1,ps,p,pks,pk)
    171      endif
    172      CALL massdair(p,masse)
    173   ENDIF
    174 
    175   IF (llm == 1) THEN
    176      ! initialize fields for the shallow water case, if required
    177      IF (.NOT.read_start) THEN
    178         phis(:)=0.
    179         q(:,:,:)=0
    180         CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
    181      endif
    182   ENDIF
    183 
    184   academic_case: if (iflag_phys >= 2) THEN
    185      ! initializations
    186 
    187      ! 1. local parameters
    188      ! by convention, winter is in the southern hemisphere
    189      ! Geostrophic wind or no wind?
    190      ok_geost=.TRUE.
    191      CALL getin('ok_geost',ok_geost)
    192      ! Constants for Newtonian relaxation and friction
    193      k_f=1.                !friction
    194      CALL getin('k_j',k_f)
    195      k_f=1./(daysec*k_f)
    196      k_c_s=4.  !cooling surface
    197      CALL getin('k_c_s',k_c_s)
    198      k_c_s=1./(daysec*k_c_s)
    199      k_c_a=40. !cooling free atm
    200      CALL getin('k_c_a',k_c_a)
    201      k_c_a=1./(daysec*k_c_a)
    202      ! Constants for Teta equilibrium profile
    203      teta0=315.     ! mean Teta (S.H. 315K)
    204      CALL getin('teta0',teta0)
    205      ttp=200.       ! Tropopause temperature (S.H. 200K)
    206      CALL getin('ttp',ttp)
    207      eps=0.         ! Deviation to N-S symmetry(~0-20K)
    208      CALL getin('eps',eps)
    209      delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
    210      CALL getin('delt_y',delt_y)
    211      delt_z=10.     ! Vertical Gradient (S.H. 10K)
    212      CALL getin('delt_z',delt_z)
    213      ! Polar vortex
    214      ok_pv=.FALSE.
    215      CALL getin('ok_pv',ok_pv)
    216      phi_pv=-50.            ! Latitude of edge of vortex
    217      CALL getin('phi_pv',phi_pv)
    218      phi_pv=phi_pv*pi/180.
    219      dphi_pv=5.             ! Width of the edge
    220      CALL getin('dphi_pv',dphi_pv)
    221      dphi_pv=dphi_pv*pi/180.
    222      gam_pv=4.              ! -dT/dz vortex (in K/km)
    223      CALL getin('gam_pv',gam_pv)
    224      tetanoise=0.005
    225      CALL getin('tetanoise',tetanoise)
    226 
    227      ! 2. Initialize fields towards which to relax
    228      ! Friction
    229      knewt_g=k_c_a
    230      DO l=1,llm
    231         zsig=presnivs(l)/preff
    232         knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
    233         kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
    234      ENDDO
    235      DO j=1,jjp1
    236         clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
    237      ENDDO
    238 
    239      ! Potential temperature
    240      DO l=1,llm
    241         zsig=presnivs(l)/preff
    242         tetastrat=ttp*zsig**(-kappa)
    243         tetapv=tetastrat
    244         IF ((ok_pv).AND.(zsig<0.1)) THEN
    245            tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
    246         ENDIF
    247         DO j=1,jjp1
    248            ! Troposphere
    249            ddsin=sin(rlatu(j))
    250            tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin &
    251                 -delt_z*(1.-ddsin*ddsin)*log(zsig)
    252            IF (planet_type=="giant") THEN
    253              tetajl(j,l)=teta0+(delt_y*                   &
    254                 ((sin(rlatu(j)*3.14159*eps+0.0001))**2)   &
    255                 / ((rlatu(j)*3.14159*eps+0.0001)**2))     &
    256                 -delt_z*log(zsig)
    257            endif
    258            ! Profil stratospherique isotherme (+vortex)
    259            w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
    260            tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
    261            tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 
    262         ENDDO
    263      ENDDO
    264 
    265      !          CALL writefield('theta_eq',tetajl)
    266 
    267      do l=1,llm
    268         do j=1,jjp1
    269            do i=1,iip1
    270               ij=(j-1)*iip1+i
    271               tetarappel(ij,l)=tetajl(j,l)
    272            enddo
     272      enddo
     273    enddo
     274
     275    ! 3. Initialize fields (if necessary)
     276    IF (.NOT. read_start) THEN
     277      ! bulk initialization of temperature
     278      IF (iflag_phys>10000) THEN
     279        ! Particular case to impose a constant temperature T0=0.01*iflag_physx
     280        teta(:, :) = 0.01 * iflag_phys / (pk(:, :) / cpp)
     281      ELSE
     282        teta(:, :) = tetarappel(:, :)
     283      ENDIF
     284      ! geopotential
     285      CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
     286
     287      DO l = 1, llm
     288        PRINT*, 'presnivs,play,l', presnivs(l), (pk(1, l) / cpp)**(1. / kappa) * preff
     289        !pks(ij) = (cpp/preff) * ps(ij)
     290        !pk(ij,1) = .5*pks(ij)
     291        ! pk = cpp * (p/preff)^kappa
     292      ENDDO
     293
     294      ! winds
     295      IF (ok_geost) THEN
     296        CALL ugeostr(phi, ucov)
     297      else
     298        ucov(:, :) = 0.
     299      endif
     300      vcov(:, :) = 0.
     301
     302      ! bulk initialization of tracers
     303      IF (planet_type=="earth") THEN
     304        ! Earth: first two tracers will be water
     305        do iq = 1, nqtot
     306          q(:, :, iq) = 0.
     307          IF(tracers(iq)%name == addPhase('H2O', 'g')) q(:, :, iq) = 1.e-10
     308          IF(tracers(iq)%name == addPhase('H2O', 'l')) q(:, :, iq) = 1.e-15
     309
     310          ! CRisi: init des isotopes
     311          ! distill de Rayleigh très simplifiée
     312          iName = tracers(iq)%iso_iName
     313          IF (niso <= 0 .OR. iName <= 0) CYCLE
     314          iPhase = tracers(iq)%iso_iPhase
     315          iqParent = tracers(iq)%iqParent
     316          IF(tracers(iq)%iso_iZone == 0) THEN
     317            IF (tnat1) THEN
     318              tnat = 1.0
     319              alpha_ideal = 1.0
     320              WRITE(*, *) 'Attention dans iniacademic: alpha_ideal=1'
     321            else
     322              IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
     323                      CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
     324            endif
     325            q(:, :, iq) = q(:, :, iqParent) * tnat * (q(:, :, iqParent) / 30.e-3)**(alpha_ideal - 1.)
     326          ELSE !IF(tracers(iq)%iso_iZone == 0) THEN
     327            IF(tracers(iq)%iso_iZone == 1) THEN
     328              ! correction le 14 mai 2024 pour que tous les traceurs soient de la couleur 1.
     329              ! Sinon, on va avoir des porblèmes de conservation de masse de traceurs.
     330              q(:, :, iq) = q(:, :, iqIsoPha(iName, iPhase))
     331            else !IF(tracers(iq)%iso_iZone == 1) THEN
     332              q(:, :, iq) = 0.
     333            endif !IF(tracers(iq)%iso_iZone == 1) THEN
     334          END IF !IF(tracers(iq)%iso_iZone == 0) THEN
    273335        enddo
    274      enddo
    275 
    276      ! 3. Initialize fields (if necessary)
    277      IF (.NOT. read_start) THEN
    278         ! bulk initialization of temperature
    279         IF (iflag_phys>10000) THEN
    280         ! Particular case to impose a constant temperature T0=0.01*iflag_physx
    281            teta(:,:)= 0.01*iflag_phys/(pk(:,:)/cpp)
    282         ELSE
    283            teta(:,:)=tetarappel(:,:)
    284         ENDIF
    285         ! geopotential
    286         CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
    287 
    288         DO l=1,llm
    289           PRINT*,'presnivs,play,l',presnivs(l),(pk(1,l)/cpp)**(1./kappa)*preff
    290          !pks(ij) = (cpp/preff) * ps(ij)
    291          !pk(ij,1) = .5*pks(ij)
    292          ! pk = cpp * (p/preff)^kappa
    293         ENDDO
    294 
    295         ! winds
    296         IF (ok_geost) THEN
    297            CALL ugeostr(phi,ucov)
    298         else
    299            ucov(:,:)=0.
    300         endif
    301         vcov(:,:)=0.
    302 
    303         ! bulk initialization of tracers
    304         IF (planet_type=="earth") THEN
    305            ! Earth: first two tracers will be water
    306            do iq=1,nqtot
    307               q(:,:,iq)=0.
    308               IF(tracers(iq)%name == addPhase('H2O', 'g')) q(:,:,iq)=1.e-10
    309               IF(tracers(iq)%name == addPhase('H2O', 'l')) q(:,:,iq)=1.e-15
    310 
    311               ! CRisi: init des isotopes
    312               ! distill de Rayleigh très simplifiée
    313               iName    = tracers(iq)%iso_iName
    314               IF (niso <= 0 .OR. iName <= 0) CYCLE
    315               iPhase   = tracers(iq)%iso_iPhase
    316               iqParent = tracers(iq)%iqParent
    317               IF(tracers(iq)%iso_iZone == 0) THEN
    318                  IF (tnat1) THEN
    319                          tnat=1.0
    320                          alpha_ideal=1.0
    321                          WRITE(*,*) 'Attention dans iniacademic: alpha_ideal=1'
    322                  else
    323                     IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
    324                     CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
    325                  endif
    326                  q(:,:,iq) = q(:,:,iqParent)*tnat*(q(:,:,iqParent)/30.e-3)**(alpha_ideal-1.)
    327               ELSE !IF(tracers(iq)%iso_iZone == 0) THEN
    328                  IF(tracers(iq)%iso_iZone == 1) THEN
    329                     ! correction le 14 mai 2024 pour que tous les traceurs soient de la couleur 1.
    330                     ! Sinon, on va avoir des porblèmes de conservation de masse de traceurs.
    331                     q(:,:,iq) = q(:,:,iqIsoPha(iName,iPhase))
    332                  else !IF(tracers(iq)%iso_iZone == 1) THEN
    333                     q(:,:,iq) = 0.
    334                  endif !IF(tracers(iq)%iso_iZone == 1) THEN
    335               END IF !IF(tracers(iq)%iso_iZone == 0) THEN
    336            enddo
    337         else
    338            q(:,:,:)=0
    339         endif ! of if (planet_type=="earth")
    340 
    341         CALL check_isotopes_seq(q,1,ip1jmp1,'iniacademic_loc')
    342 
    343         ! add random perturbation to temperature
    344         idum  = -1
    345         zz = ran1(idum)
    346         idum  = 0
    347         do l=1,llm
    348            do ij=iip2,ip1jm
    349               teta(ij,l)=teta(ij,l)*(1.+tetanoise*ran1(idum))
    350            enddo
     336      else
     337        q(:, :, :) = 0
     338      endif ! of if (planet_type=="earth")
     339
     340      CALL check_isotopes_seq(q, 1, ip1jmp1, 'iniacademic_loc')
     341
     342      ! add random perturbation to temperature
     343      idum = -1
     344      zz = ran1(idum)
     345      idum = 0
     346      do l = 1, llm
     347        do ij = iip2, ip1jm
     348          teta(ij, l) = teta(ij, l) * (1. + tetanoise * ran1(idum))
    351349        enddo
    352 
    353         ! maintain periodicity in longitude
    354         do l=1,llm
    355            do ij=1,ip1jmp1,iip1
    356               teta(ij+iim,l)=teta(ij,l)
    357            enddo
     350      enddo
     351
     352      ! maintain periodicity in longitude
     353      do l = 1, llm
     354        do ij = 1, ip1jmp1, iip1
     355          teta(ij + iim, l) = teta(ij, l)
    358356        enddo
    359 
    360      ENDIF ! of IF (.NOT. read_start)
     357      enddo
     358
     359    ENDIF ! of IF (.NOT. read_start)
    361360  ENDIF academic_case
    362361

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/integrd.F90

    r5117 r5118  
    1111  USE comvert_mod, ONLY: ap, bp
    1212  USE temps_mod, ONLY: dt
     13  USE lmdz_iniprint, ONLY: lunout, prt_level
    1314
    1415  IMPLICIT NONE
    … …  
    3334  include "paramet.h"
    3435  include "comgeom.h"
    35   include "iniprint.h"
    3636
    3737  !   Arguments:
    … …  
    250250  END IF
    251251
    252 
    253252END SUBROUTINE integrd

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/leapfrog.F90

    r5117 r5118  
    2626  USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_PHYS
    2727  USE lmdz_description, ONLY: descript
     28  USE lmdz_iniprint, ONLY: lunout, prt_level
    2829
    2930  IMPLICIT NONE
    … …  
    6465  include "comdissnew.h"
    6566  include "comgeom.h"
    66   include "iniprint.h"
    6767  include "academic.h"
    6868
    … …  
    227227  !   -----
    228228
    229   jD_cur = jD_ref + day_ini - day_ref +                             &
    230   (itau+1)/day_step
    231   jH_cur = jH_ref + start_time +                                    &
    232           mod(itau+1, day_step)/float(day_step)
     229  jD_cur = jD_ref + day_ini - day_ref + &
     230          (itau + 1) / day_step
     231  jH_cur = jH_ref + start_time + &
     232          mod(itau + 1, day_step) / float(day_step)
    233233  jD_cur = jD_cur + int(jH_cur)
    234234  jH_cur = jH_cur - int(jH_cur)
    … …  
    237237
    238238  IF (ok_guide) THEN
    239     CALL guide_main(itau,ucov,vcov,teta,q,masse,ps)
     239    CALL guide_main(itau, ucov, vcov, teta, q, masse, ps)
    240240  ENDIF
    241241
    … …  
    346346      !maf stokage du flux de masse pour traceurs OFF-LINE
    347347
    348        CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis, &
    349              dtvr, itau)
    350 
     348      CALL fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, &
     349              dtvr, itau)
    351350
    352351    ENDIF ! of IF (offline)
    … …  
    405404    !       jH_cur = jH_ref +                                            &
    406405    ! &              (itau * dtvr / daysec - int(itau * dtvr / daysec))
    407     jD_cur = jD_ref + day_ini - day_ref +                        &
    408     (itau+1)/day_step
     406    jD_cur = jD_ref + day_ini - day_ref + &
     407            (itau + 1) / day_step
    409408
    410409    IF (planet_type =="generic") THEN
    … …  
    413412    ENDIF
    414413
    415     jH_cur = jH_ref + start_time +                               &
    416             mod(itau+1, day_step)/float(day_step)
     414    jH_cur = jH_ref + start_time + &
     415            mod(itau + 1, day_step) / float(day_step)
    417416    jD_cur = jD_cur + int(jH_cur)
    418417    jH_cur = jH_cur - int(jH_cur)
    … …  
    650649
    651650      IF (ok_dynzon) THEN
    652              CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav, &
    653                    ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
     651        CALL bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &
     652                ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
    654653      END IF
    655654      IF (ok_dyn_ave) THEN
    656              CALL writedynav(itau,vcov, &
    657                    ucov,teta,pk,phi,q,masse,ps,phis)
     655        CALL writedynav(itau, vcov, &
     656                ucov, teta, pk, phi, q, masse, ps, phis)
    658657      ENDIF
    659658
    … …  
    675674          vnat(:, l) = vcov(:, l) / cv(:)
    676675        enddo
    677           IF (ok_dyn_ins) THEN
    678             ! WRITE(lunout,*) "leapfrog: CALL writehist, itau=",itau
    679            CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
    680             ! CALL WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
    681             ! CALL WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
    682            ! CALL WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
    683            !  CALL WriteField('ps',reshape(ps,(/iip1,jmp1/)))
    684            !  CALL WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
    685           endif ! of if (ok_dyn_ins)
     676        IF (ok_dyn_ins) THEN
     677          ! WRITE(lunout,*) "leapfrog: CALL writehist, itau=",itau
     678          CALL writehist(itau, vcov, ucov, teta, phi, q, masse, ps, phis)
     679          ! CALL WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
     680          ! CALL WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
     681          ! CALL WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
     682          !  CALL WriteField('ps',reshape(ps,(/iip1,jmp1/)))
     683          !  CALL WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
     684        endif ! of if (ok_dyn_ins)
    686685        ! For some Grads outputs of fields
    687686        IF (output_grads_dyn) THEN
    … …  
    781780
    782781        IF (ok_dynzon) THEN
    783              CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav, &
    784                    ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
     782          CALL bilan_dyn(2, dtvr * iperiod, dtvr * day_step * periodav, &
     783                  ps, masse, pk, pbaru, pbarv, teta, phi, ucov, vcov, q)
    785784        ENDIF
    786785        IF (ok_dyn_ave) THEN
    787              CALL writedynav(itau,vcov, &
    788                    ucov,teta,pk,phi,q,masse,ps,phis)
     786          CALL writedynav(itau, vcov, &
     787                  ucov, teta, pk, phi, q, masse, ps, phis)
    789788        ENDIF
    790789
    … …  
    799798          vnat(:, l) = vcov(:, l) / cv(:)
    800799        enddo
    801           IF (ok_dyn_ins) THEN
    802              ! WRITE(lunout,*) "leapfrog: CALL writehist (b)",
    803   ! &                        itau,iecri
    804             CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
    805           endif ! of if (ok_dyn_ins)
     800        IF (ok_dyn_ins) THEN
     801          ! WRITE(lunout,*) "leapfrog: CALL writehist (b)",
     802          ! &                        itau,iecri
     803          CALL writehist(itau, vcov, ucov, teta, phi, q, masse, ps, phis)
     804        endif ! of if (ok_dyn_ins)
    806805        ! For some Grads outputs
    807806        IF (output_grads_dyn) THEN

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/sw_case_williamson91_6.F90

    r5116 r5118  
    2727  USE comconst_mod, ONLY: cpp, omeg, rad
    2828  USE comvert_mod, ONLY: ap, bp, preff
     29  USE lmdz_iniprint, ONLY: lunout, prt_level
    2930
    3031  IMPLICIT NONE
    … …  
    3637  include "paramet.h"
    3738  include "comgeom.h"
    38   include "iniprint.h"
    3939
    4040  !   Arguments:

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/top_bound.F90

    r5117 r5118  
    66          tau_top_bound
    77  USE comvert_mod, ONLY: presnivs, preff, scaleheight
     8  USE lmdz_iniprint, ONLY: lunout, prt_level
    89
    910  IMPLICIT NONE
    … …  
    5354
    5455  include "comdissipn.h"
    55   include "iniprint.h"
    5656
    5757  !   Arguments:

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/vlsplt.F90

    r5117 r5118  
    107107  USE infotrac, ONLY: nqtot,tracers, & ! CRisi
    108108        min_qParent,min_qMass,min_ratio ! MVals et CRisi
     109  USE lmdz_iniprint, ONLY: lunout, prt_level
    109110
    110111  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    121122  include "dimensions.h"
    122123  include "paramet.h"
    123   include "iniprint.h"
    124124  !
    125125  !
Note: See TracChangeset for help on using the changeset viewer.