Ignore:
Timestamp:
Sep 7, 2012, 2:49:58 PM (12 years ago)
Author:
emillour
Message:

Common dynamics: updates to keep up with LMDZ5 Earth (rev 1649)
See file "DOC/chantiers/commit_importants.log" for details.
EM

Location:
trunk/LMDZ.COMMON/libf/dyn3d
Files:
17 edited
1 moved

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.COMMON/libf/dyn3d/calfis.F

    r108 r776  
    170170      PARAMETER(ntetaSTD=3)
    171171      REAL rtetaSTD(ntetaSTD)
    172       DATA rtetaSTD/350., 380., 405./
     172      DATA rtetaSTD/350., 380., 405./ ! Earth-specific values, beware !!
    173173      REAL PVteta(ngridmx,ntetaSTD)
    174174c
     
    461461      if (planet_type=="earth") then
    462462#ifdef CPP_EARTH
     463! PVtheta calls tetalevel, which is in the (Earth) physics
    463464cIM calcul PV a teta=350, 380, 405K
    464465      CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
     
    482483! ne pose pas de probleme a priori.
    483484
    484 #ifdef CPP_PHYS
    485 
    486485!      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
    487486      zdt_split=dtphys/nsplit_phys
     
    490489      zdtfic(:,:)=0.
    491490      zdqfic(:,:,:)=0.
     491
     492#ifdef CPP_PHYS
    492493
    493494      do isplit=1,nsplit_phys
     
    563564         zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
    564565
    565       enddo
     566      enddo ! of do isplit=1,nsplit_phys
     567
     568#endif
     569! #endif of #ifdef CPP_PHYS
     570
    566571      zdufi(:,:)=zdufic(:,:)/nsplit_phys
    567572      zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
     
    569574      zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
    570575
    571 #endif
    572 ! #endif of #ifdef CPP_PHYS
    573576
    574577500   CONTINUE
  • trunk/LMDZ.COMMON/libf/dyn3d/ce0l.F90

    r492 r776  
    2828  IMPLICIT NONE
    2929#ifndef CPP_EARTH
    30   WRITE(*,*)'limit_netcdf: Earth-specific program, needs Earth physics'
     30#include "iniprint.h"
     31  WRITE(lunout,*)'limit_netcdf: Earth-specific program, needs Earth physics'
    3132#else
    3233!-------------------------------------------------------------------------------
  • trunk/LMDZ.COMMON/libf/dyn3d/comvert.h

    r127 r776  
    11!
    2 ! $Id: comvert.h 1520 2011-05-23 11:37:09Z emillour $
     2! $Id: comvert.h 1625 2012-05-09 13:14:48Z lguez $
    33!
    44!-----------------------------------------------------------------------
     
    99     &               aps(llm),bps(llm),scaleheight
    1010
    11       common/comverti/disvert_type
     11      common/comverti/disvert_type, pressure_exner
    1212
    1313      real ap     ! hybrid pressure contribution at interlayers
     
    3030                           !     using 'z2sig.def' (or 'esasig.def) file
    3131
     32      logical pressure_exner
     33!     compute pressure inside layers using Exner function, else use mean
     34!     of pressure values at interfaces
     35
    3236 !-----------------------------------------------------------------------
  • trunk/LMDZ.COMMON/libf/dyn3d/disvert.F90

    r128 r776  
    1 ! $Id: disvert.F90 1520 2011-05-23 11:37:09Z emillour $
     1! $Id: disvert.F90 1645 2012-07-30 16:01:50Z lguez $
    22
    33SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,scaleheight)
    44
    55  ! Auteur : P. Le Van
     6
     7  use new_unit_m, only: new_unit
     8  use ioipsl, only: getin
     9  use assert_m, only: assert
    610
    711  IMPLICIT NONE
     
    1822
    1923  real,intent(in) :: pa, preff
    20   real,intent(out) :: ap(llmp1), bp(llmp1)
     24  real,intent(out) :: ap(llmp1) ! in Pa
     25  real, intent(out):: bp(llmp1)
    2126  real,intent(out) :: dpres(llm), nivsigs(llm), nivsig(llmp1)
    2227  real,intent(out) :: presnivs(llm)
     
    2631  real zk, zkm1, dzk1, dzk2, k0, k1
    2732
    28   INTEGER l
     33  INTEGER l, unit
    2934  REAL dsigmin
    3035  REAL alpha, beta, deltaz
    31   INTEGER iostat
    3236  REAL x
    3337  character(len=*),parameter :: modname="disvert"
    3438
     39  character(len=6):: vert_sampling
     40  ! (allowed values are "param", "tropo", "strato" and "read")
     41
    3542  !-----------------------------------------------------------------------
     43
     44  print *, "Call sequence information: disvert"
    3645
    3746  ! default scaleheight is 8km for earth
    3847  scaleheight=8.
    3948
    40   OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
     49  vert_sampling = merge("strato", "tropo ", ok_strato) ! default value
     50  call getin('vert_sampling', vert_sampling)
     51  print *, 'vert_sampling = ' // vert_sampling
    4152
    42   IF (iostat == 0) THEN
    43      ! cas 1 on lit les options dans sigma.def:
     53  select case (vert_sampling)
     54  case ("param")
     55     ! On lit les options dans sigma.def:
     56     OPEN(99, file='sigma.def', status='old', form='formatted')
    4457     READ(99, *) scaleheight ! hauteur d'echelle 8.
    4558     READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
     
    6982     sig(llm+1)=0.
    7083
    71      DO l = 1, llm
    72         dsig(l) = sig(l)-sig(l+1)
    73      end DO
    74   ELSE
    75      if (ok_strato) then
    76         if (llm==39) then
    77            dsigmin=0.3
    78         else if (llm==50) then
    79            dsigmin=1.
    80         else
    81            write(lunout,*) trim(modname), &
    82            ' ATTENTION discretisation z a ajuster'
    83            dsigmin=1.
    84         endif
    85         write(lunout,*) trim(modname), &
    86         ' Discretisation verticale DSIGMIN=',dsigmin
    87      endif
     84     bp(: llm) = EXP(1. - 1. / sig(: llm)**2)
     85     bp(llmp1) = 0.
    8886
     87     ap = pa * (sig - bp)
     88  case("tropo")
    8989     DO l = 1, llm
    9090        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
    91 
    92         IF (ok_strato) THEN
    93            dsig(l) =(dsigmin + 7. * SIN(x)**2) &
    94                 *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
    95         ELSE
    96            dsig(l) = 1.0 + 7.0 * SIN(x)**2
    97         ENDIF
     91        dsig(l) = 1.0 + 7.0 * SIN(x)**2
    9892     ENDDO
    9993     dsig = dsig / sum(dsig)
     
    10296        sig(l) = sig(l+1) + dsig(l)
    10397     ENDDO
    104   ENDIF
     98
     99     bp(1)=1.
     100     bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
     101     bp(llmp1) = 0.
     102
     103     ap(1)=0.
     104     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
     105  case("strato")
     106     if (llm==39) then
     107        dsigmin=0.3
     108     else if (llm==50) then
     109        dsigmin=1.
     110     else
     111        write(lunout,*) trim(modname), ' ATTENTION discretisation z a ajuster'
     112        dsigmin=1.
     113     endif
     114     WRITE(LUNOUT,*) trim(modname), 'Discretisation verticale DSIGMIN=',dsigmin
     115
     116     DO l = 1, llm
     117        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
     118        dsig(l) =(dsigmin + 7. * SIN(x)**2) &
     119             *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
     120     ENDDO
     121     dsig = dsig / sum(dsig)
     122     sig(llm+1) = 0.
     123     DO l = llm, 1, -1
     124        sig(l) = sig(l+1) + dsig(l)
     125     ENDDO
     126
     127     bp(1)=1.
     128     bp(2: llm) = EXP(1. - 1. / sig(2: llm)**2)
     129     bp(llmp1) = 0.
     130
     131     ap(1)=0.
     132     ap(2: llm + 1) = pa * (sig(2: llm + 1) - bp(2: llm + 1))
     133  case("read")
     134     ! Read "ap" and "bp". First line is skipped (title line). "ap"
     135     ! should be in Pa. First couple of values should correspond to
     136     ! the surface, that is : "bp" should be in descending order.
     137     call new_unit(unit)
     138     open(unit, file="hybrid.txt", status="old", action="read", &
     139          position="rewind")
     140     read(unit, fmt=*) ! skip title line
     141     do l = 1, llm + 1
     142        read(unit, fmt=*) ap(l), bp(l)
     143     end do
     144     close(unit)
     145     call assert(ap(1) == 0., ap(llm + 1) == 0., bp(1) == 1., &
     146          bp(llm + 1) == 0., "disvert: bad ap or bp values")
     147  case default
     148     call abort_gcm("disvert", 'Wrong value for "vert_sampling"', 1)
     149  END select
    105150
    106151  DO l=1, llm
     
    111156     nivsig(l)= REAL(l)
    112157  ENDDO
    113 
    114   ! .... Calculs de ap(l) et de bp(l) ....
    115   ! ..... pa et preff sont lus sur les fichiers start par lectba .....
    116 
    117   bp(llmp1) = 0.
    118 
    119   DO l = 1, llm
    120      bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
    121      ap(l) = pa * ( sig(l) - bp(l) )
    122   ENDDO
    123 
    124   bp(1)=1.
    125   ap(1)=0.
    126 
    127   ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
    128158
    129159  write(lunout, *)  trim(modname),': BP '
  • trunk/LMDZ.COMMON/libf/dyn3d/dynetat0.F

    r492 r776  
    66
    77      USE infotrac
     8      use netcdf, only: nf90_get_var
    89      IMPLICIT NONE
    910
     
    2829#include "comconst.h"
    2930#include "comvert.h"
    30 #include "comgeom.h"
     31#include "comgeom2.h"
    3132#include "ener.h"
    3233#include "netcdf.inc"
     
    4041
    4142      CHARACTER*(*) fichnom
    42       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm)
    43       REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
    44       REAL ps(ip1jmp1),phis(ip1jmp1)
     43      REAL vcov(iip1, jjm,llm),ucov(iip1, jjp1,llm),teta(iip1, jjp1,llm)
     44      REAL q(iip1,jjp1,llm,nqtot),masse(iip1, jjp1,llm)
     45      REAL ps(iip1, jjp1),phis(iip1, jjp1)
    4546
    4647      REAL time
     
    7071         CALL abort
    7172      ENDIF
    72 #ifdef NC_DOUBLE
    73       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
    74 #else
    75       ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
    76 #endif
     73      ierr = nf90_get_var(nid, nvarid, tab_cntrl)
    7774      IF (ierr .NE. NF_NOERR) THEN
    7875         write(lunout,*)"dynetat0: Lecture echoue pour <controle>"
     
    142139         CALL abort
    143140      ENDIF
    144 #ifdef NC_DOUBLE
    145       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonu)
    146 #else
    147       ierr = NF_GET_VAR_REAL(nid, nvarid, rlonu)
    148 #endif
     141      ierr = nf90_get_var(nid, nvarid, rlonu)
    149142      IF (ierr .NE. NF_NOERR) THEN
    150143         write(lunout,*)"dynetat0: Lecture echouee pour <rlonu>"
     
    157150         CALL abort
    158151      ENDIF
    159 #ifdef NC_DOUBLE
    160       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatu)
    161 #else
    162       ierr = NF_GET_VAR_REAL(nid, nvarid, rlatu)
    163 #endif
     152      ierr = nf90_get_var(nid, nvarid, rlatu)
    164153      IF (ierr .NE. NF_NOERR) THEN
    165154         write(lunout,*)"dynetat0: Lecture echouee pour <rlatu>"
     
    172161         CALL abort
    173162      ENDIF
    174 #ifdef NC_DOUBLE
    175       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlonv)
    176 #else
    177       ierr = NF_GET_VAR_REAL(nid, nvarid, rlonv)
    178 #endif
     163      ierr = nf90_get_var(nid, nvarid, rlonv)
    179164      IF (ierr .NE. NF_NOERR) THEN
    180165         write(lunout,*)"dynetat0: Lecture echouee pour <rlonv>"
     
    187172         CALL abort
    188173      ENDIF
    189 #ifdef NC_DOUBLE
    190       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlatv)
    191 #else
    192       ierr = NF_GET_VAR_REAL(nid, nvarid, rlatv)
    193 #endif
     174      ierr = nf90_get_var(nid, nvarid, rlatv)
    194175      IF (ierr .NE. NF_NOERR) THEN
    195176         write(lunout,*)"dynetat0: Lecture echouee pour rlatv"
     
    202183         CALL abort
    203184      ENDIF
    204 #ifdef NC_DOUBLE
    205       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, cu)
    206 #else
    207       ierr = NF_GET_VAR_REAL(nid, nvarid, cu)
    208 #endif
     185      ierr = nf90_get_var(nid, nvarid, cu)
    209186      IF (ierr .NE. NF_NOERR) THEN
    210187         write(lunout,*)"dynetat0: Lecture echouee pour <cu>"
     
    217194         CALL abort
    218195      ENDIF
    219 #ifdef NC_DOUBLE
    220       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, cv)
    221 #else
    222       ierr = NF_GET_VAR_REAL(nid, nvarid, cv)
    223 #endif
     196      ierr = nf90_get_var(nid, nvarid, cv)
    224197      IF (ierr .NE. NF_NOERR) THEN
    225198         write(lunout,*)"dynetat0: Lecture echouee pour <cv>"
     
    232205         CALL abort
    233206      ENDIF
    234 #ifdef NC_DOUBLE
    235       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, aire)
    236 #else
    237       ierr = NF_GET_VAR_REAL(nid, nvarid, aire)
    238 #endif
     207      ierr = nf90_get_var(nid, nvarid, aire)
    239208      IF (ierr .NE. NF_NOERR) THEN
    240209         write(lunout,*)"dynetat0: Lecture echouee pour <aire>"
     
    247216         CALL abort
    248217      ENDIF
    249 #ifdef NC_DOUBLE
    250       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, phis)
    251 #else
    252       ierr = NF_GET_VAR_REAL(nid, nvarid, phis)
    253 #endif
     218      ierr = nf90_get_var(nid, nvarid, phis)
    254219      IF (ierr .NE. NF_NOERR) THEN
    255220         write(lunout,*)"dynetat0: Lecture echouee pour <phisinit>"
     
    262227         CALL abort
    263228      ENDIF
    264 #ifdef NC_DOUBLE
    265       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, time)
    266 #else
    267       ierr = NF_GET_VAR_REAL(nid, nvarid, time)
    268 #endif
     229      ierr = nf90_get_var(nid, nvarid, time)
    269230      IF (ierr .NE. NF_NOERR) THEN
    270231         write(lunout,*)"dynetat0: Lecture echouee <temps>"
     
    277238         CALL abort
    278239      ENDIF
    279 #ifdef NC_DOUBLE
    280       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ucov)
    281 #else
    282       ierr = NF_GET_VAR_REAL(nid, nvarid, ucov)
    283 #endif
     240      ierr = nf90_get_var(nid, nvarid, ucov)
    284241      IF (ierr .NE. NF_NOERR) THEN
    285242         write(lunout,*)"dynetat0: Lecture echouee pour <ucov>"
     
    292249         CALL abort
    293250      ENDIF
    294 #ifdef NC_DOUBLE
    295       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, vcov)
    296 #else
    297       ierr = NF_GET_VAR_REAL(nid, nvarid, vcov)
    298 #endif
     251      ierr = nf90_get_var(nid, nvarid, vcov)
    299252      IF (ierr .NE. NF_NOERR) THEN
    300253         write(lunout,*)"dynetat0: Lecture echouee pour <vcov>"
     
    307260         CALL abort
    308261      ENDIF
    309 #ifdef NC_DOUBLE
    310       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, teta)
    311 #else
    312       ierr = NF_GET_VAR_REAL(nid, nvarid, teta)
    313 #endif
     262      ierr = nf90_get_var(nid, nvarid, teta)
    314263      IF (ierr .NE. NF_NOERR) THEN
    315264         write(lunout,*)"dynetat0: Lecture echouee pour <teta>"
     
    325274     &                    "> est absent"
    326275           write(lunout,*)"          Il est donc initialise a zero"
    327            q(:,:,iq)=0.
     276           q(:,:,:,iq)=0.
    328277        ELSE
    329 #ifdef NC_DOUBLE
    330           ierr = NF_GET_VAR_DOUBLE(nid, nvarid, q(1,1,iq))
    331 #else
    332           ierr = NF_GET_VAR_REAL(nid, nvarid, q(1,1,iq))
    333 #endif
     278           ierr = NF90_GET_VAR(nid, nvarid, q(:,:,:,iq))
    334279          IF (ierr .NE. NF_NOERR) THEN
    335280            write(lunout,*)"dynetat0: Lecture echouee pour "//tname(iq)
     
    345290         CALL abort
    346291      ENDIF
    347 #ifdef NC_DOUBLE
    348       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, masse)
    349 #else
    350       ierr = NF_GET_VAR_REAL(nid, nvarid, masse)
    351 #endif
     292      ierr = nf90_get_var(nid, nvarid, masse)
    352293      IF (ierr .NE. NF_NOERR) THEN
    353294         write(lunout,*)"dynetat0: Lecture echouee pour <masse>"
     
    360301         CALL abort
    361302      ENDIF
    362 #ifdef NC_DOUBLE
    363       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ps)
    364 #else
    365       ierr = NF_GET_VAR_REAL(nid, nvarid, ps)
    366 #endif
     303      ierr = nf90_get_var(nid, nvarid, ps)
    367304      IF (ierr .NE. NF_NOERR) THEN
    368305         write(lunout,*)"dynetat0: Lecture echouee pour <ps>"
  • trunk/LMDZ.COMMON/libf/dyn3d/dynredem.F

    r492 r776  
    11!
    2 ! $Id: dynredem.F 1403 2010-07-01 09:02:53Z fairhead $
     2! $Id: dynredem.F 1635 2012-07-12 11:37:16Z lguez $
    33!
    44c
     
    88#endif
    99      USE infotrac
     10      use netcdf95, only: NF95_PUT_VAR
    1011 
    1112      IMPLICIT NONE
     
    1920#include "comconst.h"
    2021#include "comvert.h"
    21 #include "comgeom.h"
     22#include "comgeom2.h"
    2223#include "temps.h"
    2324#include "ener.h"
     
    3132c   ----------
    3233      INTEGER iday_end
    33       REAL phis(ip1jmp1)
     34      REAL phis(iip1, jjp1)
    3435      CHARACTER*(*) fichnom
    3536
     
    138139c
    139140      ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 27,
    140      .                       "Fichier demarrage dynamique")
     141     .                       "Fichier demmarage dynamique")
    141142c
    142143c Definir les dimensions du fichiers:
     
    166167     .                       "Parametres de controle")
    167168      ierr = NF_ENDDEF(nid)
    168 #ifdef NC_DOUBLE
    169       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
    170 #else
    171       ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
    172 #endif
     169      call NF95_PUT_VAR(nid,nvarid,tab_cntrl)
    173170c
    174171      ierr = NF_REDEF (nid)
     
    183180     .                       "Longitudes des points U")
    184181      ierr = NF_ENDDEF(nid)
    185 #ifdef NC_DOUBLE
    186       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonu)
    187 #else
    188       ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonu)
    189 #endif
     182      call NF95_PUT_VAR(nid,nvarid,rlonu)
    190183c
    191184      ierr = NF_REDEF (nid)
     
    200193     .                       "Latitudes des points U")
    201194      ierr = NF_ENDDEF(nid)
    202 #ifdef NC_DOUBLE
    203       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu)
    204 #else
    205       ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu)
    206 #endif
     195      call NF95_PUT_VAR (nid,nvarid,rlatu)
    207196c
    208197      ierr = NF_REDEF (nid)
     
    217206     .                       "Longitudes des points V")
    218207      ierr = NF_ENDDEF(nid)
    219 #ifdef NC_DOUBLE
    220       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv)
    221 #else
    222       ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv)
    223 #endif
     208      call NF95_PUT_VAR(nid,nvarid,rlonv)
    224209c
    225210      ierr = NF_REDEF (nid)
     
    234219     .                       "Latitudes des points V")
    235220      ierr = NF_ENDDEF(nid)
    236 #ifdef NC_DOUBLE
    237       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatv)
    238 #else
    239       ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatv)
    240 #endif
     221      call NF95_PUT_VAR(nid,nvarid,rlatv)
    241222c
    242223      ierr = NF_REDEF (nid)
     
    251232     .                       "Numero naturel des couches s")
    252233      ierr = NF_ENDDEF(nid)
    253 #ifdef NC_DOUBLE
    254       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,nivsigs)
    255 #else
    256       ierr = NF_PUT_VAR_REAL (nid,nvarid,nivsigs)
    257 #endif
     234      call NF95_PUT_VAR(nid,nvarid,nivsigs)
    258235c
    259236      ierr = NF_REDEF (nid)
     
    268245     .                       "Numero naturel des couches sigma")
    269246      ierr = NF_ENDDEF(nid)
    270 #ifdef NC_DOUBLE
    271       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,nivsig)
    272 #else
    273       ierr = NF_PUT_VAR_REAL (nid,nvarid,nivsig)
    274 #endif
     247      call NF95_PUT_VAR(nid,nvarid,nivsig)
    275248c
    276249      ierr = NF_REDEF (nid)
     
    285258     .                       "Coefficient A pour hybride")
    286259      ierr = NF_ENDDEF(nid)
    287 #ifdef NC_DOUBLE
    288       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,ap)
    289 #else
    290       ierr = NF_PUT_VAR_REAL (nid,nvarid,ap)
    291 #endif
     260      call NF95_PUT_VAR(nid,nvarid,ap)
    292261c
    293262      ierr = NF_REDEF (nid)
     
    302271     .                       "Coefficient B pour hybride")
    303272      ierr = NF_ENDDEF(nid)
    304 #ifdef NC_DOUBLE
    305       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bp)
    306 #else
    307       ierr = NF_PUT_VAR_REAL (nid,nvarid,bp)
    308 #endif
     273      call NF95_PUT_VAR(nid,nvarid,bp)
    309274c
    310275      ierr = NF_REDEF (nid)
     
    317282cIM 220306 END
    318283      ierr = NF_ENDDEF(nid)
    319 #ifdef NC_DOUBLE
    320       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,presnivs)
    321 #else
    322       ierr = NF_PUT_VAR_REAL (nid,nvarid,presnivs)
    323 #endif
     284      call NF95_PUT_VAR(nid,nvarid,presnivs)
    324285c
    325286c Coefficients de passage cov. <-> contra. <--> naturel
     
    338299     .                       "Coefficient de passage pour U")
    339300      ierr = NF_ENDDEF(nid)
    340 #ifdef NC_DOUBLE
    341       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cu)
    342 #else
    343       ierr = NF_PUT_VAR_REAL (nid,nvarid,cu)
    344 #endif
     301      call NF95_PUT_VAR(nid,nvarid,cu)
    345302c
    346303      ierr = NF_REDEF (nid)
     
    357314     .                       "Coefficient de passage pour V")
    358315      ierr = NF_ENDDEF(nid)
    359 #ifdef NC_DOUBLE
    360       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cv)
    361 #else
    362       ierr = NF_PUT_VAR_REAL (nid,nvarid,cv)
    363 #endif
     316      call NF95_PUT_VAR(nid,nvarid,cv)
    364317c
    365318c Aire de chaque maille:
     
    378331     .                       "Aires de chaque maille")
    379332      ierr = NF_ENDDEF(nid)
    380 #ifdef NC_DOUBLE
    381       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aire)
    382 #else
    383       ierr = NF_PUT_VAR_REAL (nid,nvarid,aire)
    384 #endif
     333      call NF95_PUT_VAR(nid,nvarid,aire)
    385334c
    386335c Geopentiel au sol:
     
    399348     .                       "Geopotentiel au sol")
    400349      ierr = NF_ENDDEF(nid)
    401 #ifdef NC_DOUBLE
    402       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,phis)
    403 #else
    404       ierr = NF_PUT_VAR_REAL (nid,nvarid,phis)
    405 #endif
     350      call NF95_PUT_VAR(nid,nvarid,phis)
    406351c
    407352c Definir les variables pour pouvoir les enregistrer plus tard:
     
    524469      USE infotrac
    525470      USE control_mod
     471      use netcdf, only: NF90_get_VAR
     472      use netcdf95, only: NF95_PUT_VAR
    526473 
    527474      IMPLICIT NONE
     
    538485#include "iniprint.h"
    539486
     487
    540488      INTEGER l
    541       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm)
    542       REAL teta(ip1jmp1,llm)                   
    543       REAL ps(ip1jmp1),masse(ip1jmp1,llm)                   
    544       REAL q(ip1jmp1,llm,nqtot)
     489      REAL vcov(iip1,jjm,llm),ucov(iip1, jjp1,llm)
     490      REAL teta(iip1, jjp1,llm)                   
     491      REAL ps(iip1, jjp1),masse(iip1, jjp1,llm)                   
     492      REAL q(iip1, jjp1, llm, nqtot)
    545493      CHARACTER*(*) fichnom
    546494     
     
    576524         CALL abort_gcm(modname,abort_message,ierr)
    577525      ENDIF
    578 #ifdef NC_DOUBLE
    579       ierr = NF_PUT_VAR1_DOUBLE (nid,nvarid,nb,time)
    580 #else
    581       ierr = NF_PUT_VAR1_REAL (nid,nvarid,nb,time)
    582 #endif
     526      call NF95_PUT_VAR(nid,nvarid,time,start=(/nb/))
    583527      write(lunout,*) "dynredem1: Enregistrement pour ", nb, time
    584528
     
    592536         CALL abort_gcm(modname,abort_message,ierr)
    593537      ENDIF
    594 #ifdef NC_DOUBLE
    595       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
    596 #else
    597       ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
    598 #endif
     538      ierr = NF90_GET_VAR(nid, nvarid, tab_cntrl)
    599539       tab_cntrl(31) = REAL(itau_dyn + itaufin)
    600 #ifdef NC_DOUBLE
    601       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
    602 #else
    603       ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
    604 #endif
     540      call NF95_PUT_VAR(nid,nvarid,tab_cntrl)
    605541
    606542c  Ecriture des champs
     
    612548         CALL abort_gcm(modname,abort_message,ierr)
    613549      ENDIF
    614 #ifdef NC_DOUBLE
    615       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,ucov)
    616 #else
    617       ierr = NF_PUT_VAR_REAL (nid,nvarid,ucov)
    618 #endif
     550      call NF95_PUT_VAR(nid,nvarid,ucov)
    619551
    620552      ierr = NF_INQ_VARID(nid, "vcov", nvarid)
     
    624556         CALL abort_gcm(modname,abort_message,ierr)
    625557      ENDIF
    626 #ifdef NC_DOUBLE
    627       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,vcov)
    628 #else
    629       ierr = NF_PUT_VAR_REAL (nid,nvarid,vcov)
    630 #endif
     558      call NF95_PUT_VAR(nid,nvarid,vcov)
    631559
    632560      ierr = NF_INQ_VARID(nid, "teta", nvarid)
     
    636564         CALL abort_gcm(modname,abort_message,ierr)
    637565      ENDIF
    638 #ifdef NC_DOUBLE
    639       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,teta)
    640 #else
    641       ierr = NF_PUT_VAR_REAL (nid,nvarid,teta)
    642 #endif
     566      call NF95_PUT_VAR(nid,nvarid,teta)
    643567
    644568      IF (type_trac == 'inca') THEN
     
    662586               CALL abort_gcm(modname,abort_message,ierr)
    663587            ENDIF
    664 #ifdef NC_DOUBLE
    665             ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,q(1,1,iq))
    666 #else
    667             ierr = NF_PUT_VAR_REAL (nid,nvarid,q(1,1,iq))
    668 #endif
     588            call NF95_PUT_VAR(nid,nvarid,q(:,:,:,iq))
    669589        ELSE ! type_trac = inca
    670590! lecture de la valeur du traceur dans start_trac.nc
     
    681601                   CALL abort_gcm(modname,abort_message,ierr)
    682602                ENDIF
    683 #ifdef NC_DOUBLE
    684                 ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,q(1,1,iq))
    685 #else
    686                 ierr = NF_PUT_VAR_REAL (nid,nvarid,q(1,1,iq))
    687 #endif
     603                call NF95_PUT_VAR(nid,nvarid,q(:,:,:,iq))
    688604               
    689605             ELSE
    690606                write(lunout,*) "dynredem1: ",trim(tname(iq)),
    691607     &              " est present dans start_trac.nc"
    692 #ifdef NC_DOUBLE
    693                ierr = NF_GET_VAR_DOUBLE(nid_trac, nvarid_trac, trac_tmp)
    694 #else
    695                ierr = NF_GET_VAR_REAL(nid_trac, nvarid_trac, trac_tmp)
    696 #endif
     608               ierr = NF90_GET_VAR(nid_trac, nvarid_trac, trac_tmp)
    697609                IF (ierr .NE. NF_NOERR) THEN
    698610                   abort_message="dynredem1: Lecture echouee pour"//
     
    708620                   CALL abort_gcm(modname,abort_message,ierr)
    709621                ENDIF
    710 #ifdef NC_DOUBLE
    711                 ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,trac_tmp)
    712 #else
    713                 ierr = NF_PUT_VAR_REAL (nid,nvarid,trac_tmp)
    714 #endif
     622                call NF95_PUT_VAR(nid, nvarid, trac_tmp)
    715623               
    716624             ENDIF ! IF (ierr .NE. NF_NOERR)
     
    725633                   CALL abort_gcm(modname,abort_message,ierr)
    726634             ENDIF
    727 #ifdef NC_DOUBLE
    728              ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,q(1,1,iq))
    729 #else
    730              ierr = NF_PUT_VAR_REAL (nid,nvarid,q(1,1,iq))
    731 #endif
     635             call NF95_PUT_VAR(nid,nvarid,q(:,:,:,iq))
    732636          ENDIF ! (ierr_file .ne. 2)
    733637       END IF   !type_trac
     
    742646         CALL abort_gcm(modname,abort_message,ierr)
    743647      ENDIF
    744 #ifdef NC_DOUBLE
    745       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,masse)
    746 #else
    747       ierr = NF_PUT_VAR_REAL (nid,nvarid,masse)
    748 #endif
     648      call NF95_PUT_VAR(nid,nvarid,masse)
    749649c
    750650      ierr = NF_INQ_VARID(nid, "ps", nvarid)
     
    754654         CALL abort_gcm(modname,abort_message,ierr)
    755655      ENDIF
    756 #ifdef NC_DOUBLE
    757       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,ps)
    758 #else
    759       ierr = NF_PUT_VAR_REAL (nid,nvarid,ps)
    760 #endif
     656      call NF95_PUT_VAR(nid,nvarid,ps)
    761657
    762658      ierr = NF_CLOSE(nid)
  • trunk/LMDZ.COMMON/libf/dyn3d/etat0_netcdf.F90

    r127 r776  
    11!
    2 ! $Id: etat0_netcdf.F90 1520 2011-05-23 11:37:09Z emillour $
     2! $Id: etat0_netcdf.F90 1625 2012-05-09 13:14:48Z lguez $
    33!
    44!-------------------------------------------------------------------------------
     
    251251!*******************************************************************************
    252252  CALL pression(ip1jmp1, ap, bp, psol, p3d)
    253   if (disvert_type.eq.1) then
     253  if (pressure_exner) then
    254254    CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y)
    255   else ! we assume that we are in the disvert_type==2 case
     255  else
    256256    CALL exner_milieu(ip1jmp1,psol,p3d,beta,pks,pk,y)
    257257  endif
  • trunk/LMDZ.COMMON/libf/dyn3d/exner_hyb.F

    r127 r776  
    5656      ! Sanity check
    5757      if (firstcall) then
    58         ! check that vertical discretization is compatible
    59         ! with this routine
    60         if (disvert_type.ne.1) then
    61           call abort_gcm(modname,
    62      &     "this routine should only be called if disvert_type==1",42)
    63         endif
    64        
    6558        ! sanity checks for Shallow Water case (1 vertical layer)
    6659        if (llm.eq.1) then
  • trunk/LMDZ.COMMON/libf/dyn3d/exner_milieu.F

    r127 r776  
    5353      ! Sanity check
    5454      if (firstcall) then
    55         ! check that vertical discretization is compatible
    56         ! with this routine
    57         if (disvert_type.ne.2) then
    58           call abort_gcm(modname,
    59      &     "this routine should only be called if disvert_type==2",42)
    60         endif
    61        
    6255        ! sanity checks for Shallow Water case (1 vertical layer)
    6356        if (llm.eq.1) then
  • trunk/LMDZ.COMMON/libf/dyn3d/gcm.F

    r492 r776  
    2121! A nettoyer. On ne veut qu'une ou deux routines d'interface
    2222! dynamique -> physique pour l'initialisation
    23 ! Ehouarn: for now these only apply to Earth:
     23! Ehouarn: the following are needed with (parallel) physics:
    2424#ifdef CPP_PHYS
    2525      USE dimphy
    2626      USE comgeomphy
    27 #endif
    28 #ifdef CPP_EARTH
    2927      USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
    3028#endif
     
    177175! A nettoyer. On ne veut qu'une ou deux routines d'interface
    178176! dynamique -> physique pour l'initialisation
    179 ! Ehouarn : temporarily (?) keep this only for Earth
    180 !      if (planet_type.eq."earth") then
    181 !#ifdef CPP_EARTH
    182177#ifdef CPP_PHYS
    183178      CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
    184179      call initcomgeomphy
    185180#endif
    186 !      endif
    187181!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    188182c
     
    465459         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys ,
    466460     ,                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp     )
    467 #endif ! CPP_PHYS
     461#endif
    468462         call_iniphys=.false.
    469463      ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
    470 !#endif
    471464
    472465c  numero de stockage pour les fichiers de redemarrage:
  • trunk/LMDZ.COMMON/libf/dyn3d/guide_mod.F90

    r127 r776  
    644644! -----------------------------------------------------------------
    645645    CALL pression( ip1jmp1, ap, bp, psi, p )
    646     if (disvert_type==1) then
     646    if (pressure_exner) then
    647647      CALL exner_hyb(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
    648     else ! we assume that we are in the disvert_type==2 case
     648    else
    649649      CALL exner_milieu(ip1jmp1,psi,p,beta,pks,pk,pkf)
    650650    endif
  • trunk/LMDZ.COMMON/libf/dyn3d/iniacademic.F90

    r492 r776  
    11!
    2 ! $Id: iniacademic.F90 1529 2011-05-26 15:17:33Z fairhead $
     2! $Id: iniacademic.F90 1625 2012-05-09 13:14:48Z lguez $
    33!
    44SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
     
    222222
    223223        CALL pression ( ip1jmp1, ap, bp, ps, p       )
    224         if (disvert_type.eq.1) then
     224        if (pressure_exner) then
    225225          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
    226         elseif (disvert_type.eq.2) then
     226        else
    227227          call exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf)
    228         else
    229           write(abort_message,*) "Wrong value for disvert_type: ", &
    230                               disvert_type
    231           call abort_gcm(modname,abort_message,0)
    232228        endif
    233229        CALL massdair(p,masse)
  • trunk/LMDZ.COMMON/libf/dyn3d/iniconst.F90

    r775 r776  
    11!
    2 ! $Id: iniconst.F 1520 2011-05-23 11:37:09Z emillour $
     2! $Id: iniconst.F90 1625 2012-05-09 13:14:48Z lguez $
    33!
    4       SUBROUTINE iniconst
     4SUBROUTINE iniconst
    55
    6       USE control_mod
     6  USE control_mod
    77#ifdef CPP_IOIPSL
    8       use IOIPSL
     8  use IOIPSL
    99#else
    10 ! if not using IOIPSL, we still need to use (a local version of) getin
    11       use ioipsl_getincom
     10  ! if not using IOIPSL, we still need to use (a local version of) getin
     11  use ioipsl_getincom
    1212#endif
    1313
    14       IMPLICIT NONE
    15 c
    16 c      P. Le Van
    17 c
    18 c-----------------------------------------------------------------------
    19 c   Declarations:
    20 c   -------------
    21 c
    22 #include "dimensions.h"
    23 #include "paramet.h"
    24 #include "comconst.h"
    25 #include "temps.h"
    26 #include "comvert.h"
    27 #include "iniprint.h"
     14  IMPLICIT NONE
     15  !
     16  !      P. Le Van
     17  !
     18  !   Declarations:
     19  !   -------------
     20  !
     21  include "dimensions.h"
     22  include "paramet.h"
     23  include "comconst.h"
     24  include "temps.h"
     25  include "comvert.h"
     26  include "iniprint.h"
    2827
     28  character(len=*),parameter :: modname="iniconst"
     29  character(len=80) :: abort_message
     30  !
     31  !
     32  !
     33  !-----------------------------------------------------------------------
     34  !   dimension des boucles:
     35  !   ----------------------
    2936
    30       character(len=*),parameter :: modname="iniconst"
    31       character(len=80) :: abort_message
    32 c
    33 c
    34 c
    35 c-----------------------------------------------------------------------
    36 c   dimension des boucles:
    37 c   ----------------------
     37  im      = iim
     38  jm      = jjm
     39  lllm    = llm
     40  imp1    = iim
     41  jmp1    = jjm + 1
     42  lllmm1  = llm - 1
     43  lllmp1  = llm + 1
    3844
    39       im      = iim
    40       jm      = jjm
    41       lllm    = llm
    42       imp1    = iim
    43       jmp1    = jjm + 1
    44       lllmm1  = llm - 1
    45       lllmp1  = llm + 1
     45  !-----------------------------------------------------------------------
    4646
    47 c-----------------------------------------------------------------------
     47  dtphys  = iphysiq * dtvr
     48  unsim   = 1./iim
     49  pi      = 2.*ASIN( 1. )
    4850
    49       dtphys  = iphysiq * dtvr
    50       unsim   = 1./iim
    51       pi      = 2.*ASIN( 1. )
     51  !-----------------------------------------------------------------------
     52  !
    5253
    53 c-----------------------------------------------------------------------
    54 c
     54  r       = cpp * kappa
    5555
    56       r       = cpp * kappa
     56  write(lunout,*) trim(modname),': R  CP  Kappa ',r,cpp,kappa
     57  !
     58  !-----------------------------------------------------------------------
    5759
    58       write(lunout,*) trim(modname),': R  CP  Kappa ',r,cpp,kappa
    59 c
    60 c-----------------------------------------------------------------------
     60  ! vertical discretization: default behavior depends on planet_type flag
     61  if (planet_type=="earth") then
     62     disvert_type=1
     63  else
     64     disvert_type=2
     65  endif
     66  ! but user can also specify using one or the other in run.def:
     67  call getin('disvert_type',disvert_type)
     68  write(lunout,*) trim(modname),': disvert_type=',disvert_type
    6169
    62 ! vertical discretization: default behavior depends on planet_type flag
    63       if (planet_type=="earth") then
    64         disvert_type=1
    65       else
    66         disvert_type=2
    67       endif
    68       ! but user can also specify using one or the other in run.def:
    69       call getin('disvert_type',disvert_type)
    70       write(lunout,*) trim(modname),': disvert_type=',disvert_type
    71      
    72       if (disvert_type==1) then
    73        ! standard case for Earth (automatic generation of levels)
    74        call disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,
    75      &              scaleheight)
    76       else if (disvert_type==2) then
    77         ! standard case for planets (levels generated using z2sig.def file)
    78         call disvert_noterre
    79       else
    80         write(abort_message,*) "Wrong value for disvert_type: ",
    81      &                        disvert_type
    82         call abort_gcm(modname,abort_message,0)
    83       endif
     70  pressure_exner = disvert_type == 1 ! default value
     71  call getin('pressure_exner', pressure_exner)
    8472
    85       END
     73  if (disvert_type==1) then
     74     ! standard case for Earth (automatic generation of levels)
     75     call disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig, scaleheight)
     76  else if (disvert_type==2) then
     77     ! standard case for planets (levels generated using z2sig.def file)
     78     call disvert_noterre
     79  else
     80     write(abort_message,*) "Wrong value for disvert_type: ", disvert_type
     81     call abort_gcm(modname,abort_message,0)
     82  endif
     83
     84END SUBROUTINE iniconst
  • trunk/LMDZ.COMMON/libf/dyn3d/inidissip.F90

    r270 r776  
    2828! Local variables:
    2929  REAL fact,zvert(llm),zz
    30   REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)
     30  REAL zh(ip1jmp1),zu(ip1jmp1), gx(ip1jmp1), divgra(ip1jmp1)
     31  real zv(ip1jm), gy(ip1jm), deltap(ip1jmp1,llm)
    3132  REAL ullm,vllm,umin,vmin,zhmin,zhmax
    32   REAL zllm,z1llm
     33  REAL zllm
    3334
    3435  INTEGER l,ij,idum,ii
     
    7879  DO l = 1,50
    7980     IF(lstardis) THEN
    80         CALL divgrad2(1,zh,deltap,niterh,zh)
     81        CALL divgrad2(1,zh,deltap,niterh,divgra)
    8182     ELSE
    82         CALL divgrad (1,zh,niterh,zh)
     83        CALL divgrad (1,zh,niterh,divgra)
    8384     ENDIF
    8485
    85      CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
    86 
    87      zllm  = ABS( zhmax )
    88      z1llm = 1./zllm
    89      DO ij = 1,ip1jmp1
    90         zh(ij) = zh(ij)* z1llm
    91      ENDDO
     86     zllm  = ABS(maxval(divgra))
     87     zh = divgra / zllm
    9288  ENDDO
    9389
     
    123119           !cccc             CALL covcont( 1,zu,zv,zu,zv )
    124120           IF(lstardis) THEN
    125               CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )
     121              CALL gradiv2( 1,zu,zv,nitergdiv,gx,gy )
    126122           ELSE
    127               CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )
     123              CALL gradiv ( 1,zu,zv,nitergdiv,gx,gy )
    128124           ENDIF
    129125        ELSE
    130126           IF(lstardis) THEN
    131               CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )
     127              CALL nxgraro2( 1,zu,zv,nitergrot,gx,gy )
    132128           ELSE
    133               CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )
     129              CALL nxgrarot( 1,zu,zv,nitergrot,gx,gy )
    134130           ENDIF
    135131        ENDIF
    136132
    137         CALL minmax(iip1*jjp1,zu,umin,ullm )
    138         CALL minmax(iip1*jjm, zv,vmin,vllm )
    139 
    140         ullm = ABS  ( ullm )
    141         vllm = ABS  ( vllm )
    142 
    143         zllm  = MAX( ullm,vllm )
    144         z1llm = 1./ zllm
    145         DO ij = 1, ip1jmp1
    146            zu(ij) = zu(ij)* z1llm
    147         ENDDO
    148         DO ij = 1, ip1jm
    149            zv(ij) = zv(ij)* z1llm
    150         ENDDO
     133        zllm = max(abs(maxval(gx)), abs(maxval(gy)))
     134        zu = gx / zllm
     135        zv = gy / zllm
    151136     end DO
    152137
  • trunk/LMDZ.COMMON/libf/dyn3d/inigrads.F

    r1 r776  
    99      implicit none
    1010
    11       integer if,im,jm,lm,i,j,l,lnblnk
     11      integer if,im,jm,lm,i,j,l
    1212      real x(im),y(jm),z(lm),fx,fy,fz,dt
    1313      real xmin,xmax,ymin,ymax
     
    4040      ivar(if)=0
    4141
    42       fichier(if)=file(1:lnblnk(file))
     42      fichier(if)=trim(file)
    4343
    4444      firsttime(if)=.true.
     
    7070
    7171      print*,4*(ifd(if)-iid(if))*(jfd(if)-jid(if))
    72       print*,file(1:lnblnk(file))//'.dat'
     72      print*,trim(file)//'.dat'
    7373
    74       OPEN (unit(if)+1,FILE=file(1:lnblnk(file))//'.dat'
     74      OPEN (unit(if)+1,FILE=trim(file)//'.dat'
    7575     s   ,FORM='unformatted',
    7676     s   ACCESS='direct'
  • trunk/LMDZ.COMMON/libf/dyn3d/integrd.F

    r270 r776  
    11!
    2 ! $Id: integrd.F 1550 2011-07-05 09:44:55Z lguez $
     2! $Id: integrd.F 1616 2012-02-17 11:59:00Z emillour $
    33!
    44      SUBROUTINE integrd
    55     $  (  nq,vcovm1,ucovm1,tetam1,psm1,massem1,
    6      $     dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis,finvmaold )
     6     $     dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis !,finvmaold
     7     &  )
    78
    89      use control_mod, only : planet_type
     
    3435#include "temps.h"
    3536#include "serre.h"
     37#include "iniprint.h"
    3638
    3739c   Arguments:
    3840c   ----------
    3941
    40       INTEGER nq
    41 
    42       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm)
    43       REAL q(ip1jmp1,llm,nq)
    44       REAL ps(ip1jmp1),masse(ip1jmp1,llm),phis(ip1jmp1)
    45 
    46       REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
    47       REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1),massem1(ip1jmp1,llm)
    48 
    49       REAL dv(ip1jm,llm),du(ip1jmp1,llm)
    50       REAL dteta(ip1jmp1,llm),dp(ip1jmp1)
    51       REAL dq(ip1jmp1,llm,nq), finvmaold(ip1jmp1,llm)
     42      integer,intent(in) :: nq ! number of tracers to handle in this routine
     43      real,intent(inout) :: vcov(ip1jm,llm) ! covariant meridional wind
     44      real,intent(inout) :: ucov(ip1jmp1,llm) ! covariant zonal wind
     45      real,intent(inout) :: teta(ip1jmp1,llm) ! potential temperature
     46      real,intent(inout) :: q(ip1jmp1,llm,nq) ! advected tracers
     47      real,intent(inout) :: ps(ip1jmp1) ! surface pressure
     48      real,intent(inout) :: masse(ip1jmp1,llm) ! atmospheric mass
     49      real,intent(in) :: phis(ip1jmp1) ! ground geopotential !!! unused
     50      ! values at previous time step
     51      real,intent(inout) :: vcovm1(ip1jm,llm)
     52      real,intent(inout) :: ucovm1(ip1jmp1,llm)
     53      real,intent(inout) :: tetam1(ip1jmp1,llm)
     54      real,intent(inout) :: psm1(ip1jmp1)
     55      real,intent(inout) :: massem1(ip1jmp1,llm)
     56      ! the tendencies to add
     57      real,intent(in) :: dv(ip1jm,llm)
     58      real,intent(in) :: du(ip1jmp1,llm)
     59      real,intent(in) :: dteta(ip1jmp1,llm)
     60      real,intent(in) :: dp(ip1jmp1)
     61      real,intent(in) :: dq(ip1jmp1,llm,nq) !!! unused
     62!      real,intent(out) :: finvmaold(ip1jmp1,llm) !!! unused
    5263
    5364c   Local:
     
    5566
    5667      REAL vscr( ip1jm ),uscr( ip1jmp1 ),hscr( ip1jmp1 ),pscr(ip1jmp1)
    57       REAL massescr( ip1jmp1,llm ), finvmasse(ip1jmp1,llm)
     68      REAL massescr( ip1jmp1,llm )
     69!      REAL finvmasse(ip1jmp1,llm)
    5870      REAL p(ip1jmp1,llmp1)
    5971      REAL tpn,tps,tppn(iim),tpps(iim)
     
    6173      REAL deltap( ip1jmp1,llm )
    6274
    63       INTEGER  l,ij,iq
     75      INTEGER  l,ij,iq,i,j
    6476
    6577      REAL SSUM
     
    88100      DO ij = 1,ip1jmp1
    89101        IF( ps(ij).LT.0. ) THEN
    90          PRINT*,' Au point ij = ',ij, ' , pression sol neg. ', ps(ij)
    91          print *, ' dans integrd'
    92          stop 1
     102         write(lunout,*) "integrd: negative surface pressure ",ps(ij)
     103         write(lunout,*) " at node ij =", ij
     104         ! since ij=j+(i-1)*jjp1 , we have
     105         j=modulo(ij,jjp1)
     106         i=1+(ij-j)/jjp1
     107         write(lunout,*) " lon = ",rlonv(i)*180./pi, " deg",
     108     &                   " lat = ",rlatu(j)*180./pi, " deg"
     109         stop
    93110        ENDIF
    94111      ENDDO
     
    110127      CALL massdair (     p  , masse         )
    111128
    112       CALL   SCOPY( ijp1llm  , masse, 1, finvmasse,  1      )
    113       CALL filtreg( finvmasse, jjp1, llm, -2, 2, .TRUE., 1  )
     129! Ehouarn : we don't use/need finvmaold and finvmasse,
     130!           so might as well not compute them
     131!      CALL   SCOPY( ijp1llm  , masse, 1, finvmasse,  1      )
     132!      CALL filtreg( finvmasse, jjp1, llm, -2, 2, .TRUE., 1  )
    114133c
    115134
     
    218237       ENDDO
    219238
    220 
    221       CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
     239! Ehouarn: forget about finvmaold
     240!      CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
    222241
    223242      endif ! of if (planet_type.eq."earth")
  • trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F

    r500 r776  
    124124
    125125      REAL  SSUM
    126       REAL time_0 , finvmaold(ip1jmp1,llm)
     126      REAL time_0
     127!     REAL finvmaold(ip1jmp1,llm)
    127128
    128129cym      LOGICAL  lafin
     
    243244      dq(:,:,:)=0.
    244245      CALL pression ( ip1jmp1, ap, bp, ps, p       )
    245       if (disvert_type==1) then
     246      if (pressure_exner) then
    246247        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
    247       else ! we assume that we are in the disvert_type==2 case
     248      else
    248249        CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
    249250      endif
     
    271272c   ----------------------------------
    272273
    273    1  CONTINUE
     274   1  CONTINUE ! Matsuno Forward step begins here
    274275
    275276      jD_cur = jD_ref + day_ini - day_ref +                             &
    276      &          int (itau * dtvr / daysec)
     277     &          itau/day_step
    277278      jH_cur = jH_ref + start_time +                                    &
    278      &          (itau * dtvr / daysec - int(itau * dtvr / daysec))
     279     &          mod(itau,day_step)/float(day_step)
    279280      jD_cur = jD_cur + int(jH_cur)
    280281      jH_cur = jH_cur - int(jH_cur)
     
    307308
    308309c   ...    P.Le Van .26/04/94  ....
    309 
    310       CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
    311       CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
    312 
    313    2  CONTINUE
     310! Ehouarn: finvmaold is actually not used
     311!      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
     312!      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
     313
     314   2  CONTINUE ! Matsuno backward or leapfrog step begins here
    314315
    315316c-----------------------------------------------------------------------
     
    357358      call tpot2t(ijp1llm,teta,temp,pk)
    358359      tsurpk = cpp*temp/pk
     360      ! compute geopotential phi()
    359361      CALL geopot  ( ip1jmp1, tsurpk  , pk , pks,  phis  , phi   )
    360362
     
    372374
    373375!      IF( forward. OR . leapf )  THEN
     376! Ehouarn: NB: at this point p with ps are not synchronized
     377!              (whereas mass and ps are...)
    374378      IF((.not.forward).OR. leapf )  THEN
    375379        ! Ehouarn: gather mass fluxes during backward Matsuno or LF step
     
    398402
    399403       CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
    400      $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
    401      $              finvmaold                                    )
     404     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
     405!     $              finvmaold                                    )
    402406
    403407       IF ((planet_type.eq."titan").and.(tidal)) then
     
    431435
    432436         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
    433          if (disvert_type==1) then
     437         if (pressure_exner) then
    434438           CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
    435          else ! we assume that we are in the disvert_type==2 case
     439         else
    436440           CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
    437441         endif
    438442
    439443           jD_cur = jD_ref + day_ini - day_ref +                        &
    440      &          int (itau * dtvr / daysec)
     444     &          itau/day_step
    441445           jH_cur = jH_ref + start_time +                               &
    442      &          (itau * dtvr / daysec - int(itau * dtvr / daysec))
     446     &          mod(itau,day_step)/float(day_step)
    443447           jD_cur = jD_cur + int(jH_cur)
    444448           jH_cur = jH_cur - int(jH_cur)
     
    545549
    546550        CALL pression ( ip1jmp1, ap, bp, ps, p                  )
    547         if (disvert_type==1) then
     551        if (pressure_exner) then
    548552          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
    549         else ! we assume that we are in the disvert_type==2 case
     553        else
    550554          CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
    551555        endif
     
    613617        ENDDO
    614618
    615         DO ij =  1,iim
    616           tppn(ij)  = aire(  ij    ) * ps (  ij    )
    617           tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
    618         ENDDO
    619           tpn  = SSUM(iim,tppn,1)/apoln
    620           tps  = SSUM(iim,tpps,1)/apols
    621 
    622         DO ij = 1, iip1
    623           ps(  ij    ) = tpn
    624           ps(ij+ip1jm) = tps
    625         ENDDO
    626 
     619        if (1 == 0) then
     620!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
     621!!!                     2) should probably not be here anyway
     622!!! but are kept for those who would want to revert to previous behaviour
     623           DO ij =  1,iim
     624             tppn(ij)  = aire(  ij    ) * ps (  ij    )
     625             tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
     626           ENDDO
     627             tpn  = SSUM(iim,tppn,1)/apoln
     628             tps  = SSUM(iim,tpps,1)/apols
     629
     630           DO ij = 1, iip1
     631             ps(  ij    ) = tpn
     632             ps(ij+ip1jm) = tps
     633           ENDDO
     634        endif ! of if (1 == 0)
    627635
    628636      END IF ! of IF(apdiss)
     
    749757
    750758              CLOSE(99)
     759              !!! Ehouarn: Why not stop here and now?
    751760            ENDIF ! of IF (itau.EQ.itaufin)
    752761
  • trunk/LMDZ.COMMON/libf/dyn3d/wrgrads.F

    r1 r776  
    2626c   local
    2727
    28       integer im,jm,lm,i,j,l,lnblnk,iv,iii,iji,iif,ijf
     28      integer im,jm,lm,i,j,l,iv,iii,iji,iif,ijf
    2929
    3030      logical writectl
     
    5959            nvar(if)=ivar(if)
    6060            var(ivar(if),if)=name
    61             tvar(ivar(if),if)=titlevar(1:lnblnk(titlevar))
     61            tvar(ivar(if),if)=trim(titlevar)
    6262            nld(ivar(if),if)=nl
    6363c           print*,'initialisation ecriture de ',var(ivar(if),if)
     
    101101      file=fichier(if)
    102102c   WARNING! on reecrase le fichier .ctl a chaque ecriture
    103       open(unit(if),file=file(1:lnblnk(file))//'.ctl'
     103      open(unit(if),file=trim(file)//'.ctl'
    104104     &         ,form='formatted',status='unknown')
    105105      write(unit(if),'(a5,1x,a40)')
    106      &       'DSET ','^'//file(1:lnblnk(file))//'.dat'
     106     &       'DSET ','^'//trim(file)//'.dat'
    107107
    108108      write(unit(if),'(a12)') 'UNDEF 1.0E30'
Note: See TracChangeset for help on using the changeset viewer.