Ignore:
Timestamp:
Oct 9, 2012, 3:35:26 PM (12 years ago)
Author:
Laurent Fairhead
Message:

Version testing basée sur la r1628

http://lmdz.lmd.jussieu.fr/utilisateurs/distribution-du-modele/versions-intermediaires

Testing release based on r1628

Location:
LMDZ5/branches/testing
Files:
6 deleted
81 edited
42 copied

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/testing

  • LMDZ5/branches/testing/libf/dyn3d/calfis.F

    r1407 r1665  
    165165      PARAMETER(ntetaSTD=3)
    166166      REAL rtetaSTD(ntetaSTD)
    167       DATA rtetaSTD/350., 380., 405./
     167      DATA rtetaSTD/350., 380., 405./ ! Earth-specific values, beware !!
    168168      REAL PVteta(ngridmx,ntetaSTD)
    169169c
    … …  
    434434c
    435435      if (planet_type=="earth") then
    436 #ifdef CPP_EARTH
     436#ifdef CPP_PHYS
     437! PVtheta calls tetalevel, which is in the physics
    437438cIM calcul PV a teta=350, 380, 405K
    438439      CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
    … …  
    450451
    451452
    452       if (planet_type=="earth") then
    453 #ifdef CPP_EARTH
    454453
    455454!      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
    … …  
    460459      zdqfic(:,:,:)=0.
    461460
    462       do isplit=1,nsplit_phys
     461      if (planet_type=="earth") then
     462#ifdef CPP_PHYS
     463
     464       do isplit=1,nsplit_phys
    463465
    464466         jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys)
    … …  
    503505         zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:)
    504506
    505       enddo
     507       enddo ! of do isplit=1,nsplit_phys
     508
     509#endif
     510! of #ifdef CPP_PHYS
     511      endif ! of if (planet_type=="earth")
     512
    506513      zdufi(:,:)=zdufic(:,:)/nsplit_phys
    507514      zdvfi(:,:)=zdvfic(:,:)/nsplit_phys
    … …  
    509516      zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys
    510517
    511 #endif
    512       endif !of if (planet_type=="earth")
    513518
    514519500   CONTINUE

  • LMDZ5/branches/testing/libf/dyn3d/ce0l.F90

    r1664 r1665  
    2828  IMPLICIT NONE
    2929#ifndef CPP_EARTH
     30#include "iniprint.h"
    3031  WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
    3132#else

  • LMDZ5/branches/testing/libf/dyn3d/comvert.h

    r1520 r1665  
    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 !-----------------------------------------------------------------------

  • LMDZ5/branches/testing/libf/dyn3d/conf_gcm.F

    r1664 r1665  
    155155      nday = 10
    156156      CALL getin('nday',nday)
     157
     158!Config  Key  = starttime
     159!Config  Desc = Heure de depart de la simulation
     160!Config  Def  = 0
     161!Config  Help = Heure de depart de la simulation
     162!Config         en jour
     163      starttime = 0
     164      CALL getin('starttime',starttime)
    157165
    158166!Config  Key  = day_step

  • LMDZ5/branches/testing/libf/dyn3d/control_mod.F90

    r1502 r1665  
    1010  IMPLICIT NONE
    1111
    12   REAL    :: periodav
     12  REAL    :: periodav, starttime
    1313  INTEGER :: nday,day_step,iperiod,iapp_tracvl,nsplit_phys
    1414  INTEGER :: iconser,iecri,dissip_period,iphysiq,iecrimoy

  • LMDZ5/branches/testing/libf/dyn3d/disvert.F90

    r1524 r1665  
    44
    55  ! Auteur : P. Le Van
     6
     7  use new_unit_m, only: new_unit
     8  use ioipsl, only: getin
    69
    710  IMPLICIT NONE
    … …  
    2629  real zk, zkm1, dzk1, dzk2, k0, k1
    2730
    28   INTEGER l
     31  INTEGER l, unit
    2932  REAL dsigmin
    3033  REAL alpha, beta, deltaz
    31   INTEGER iostat
    3234  REAL x
    3335  character(len=*),parameter :: modname="disvert"
    3436
     37  character(len=6):: vert_sampling
     38  ! (allowed values are "param", "tropo", "strato" and "read")
     39
    3540  !-----------------------------------------------------------------------
     41
     42  print *, "Call sequence information: disvert"
    3643
    3744  ! default scaleheight is 8km for earth
    3845  scaleheight=8.
    3946
    40   OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
     47  vert_sampling = merge("strato", "tropo ", ok_strato) ! default value
     48  call getin('vert_sampling', vert_sampling)
     49  print *, 'vert_sampling = ' // vert_sampling
    4150
    42   IF (iostat == 0) THEN
    43      ! cas 1 on lit les options dans sigma.def:
     51  select case (vert_sampling)
     52  case ("param")
     53     ! On lit les options dans sigma.def:
     54     OPEN(99, file='sigma.def', status='old', form='formatted')
    4455     READ(99, *) scaleheight ! hauteur d'echelle 8.
    4556     READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
    … …  
    6879
    6980     sig(llm+1)=0.
    70 
    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
    88 
     81  case("tropo")
    8982     DO l = 1, llm
    9083        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
     84        dsig(l) = 1.0 + 7.0 * SIN(x)**2
    9885     ENDDO
    9986     dsig = dsig / sum(dsig)
    … …  
    10289        sig(l) = sig(l+1) + dsig(l)
    10390     ENDDO
    104   ENDIF
     91  case("strato")
     92     if (llm==39) then
     93        dsigmin=0.3
     94     else if (llm==50) then
     95        dsigmin=1.
     96     else
     97        write(lunout,*) trim(modname), ' ATTENTION discretisation z a ajuster'
     98        dsigmin=1.
     99     endif
     100     WRITE(LUNOUT,*) trim(modname), 'Discretisation verticale DSIGMIN=',dsigmin
     101
     102     DO l = 1, llm
     103        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
     104        dsig(l) =(dsigmin + 7. * SIN(x)**2) &
     105             *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
     106     ENDDO
     107     dsig = dsig / sum(dsig)
     108     sig(llm+1) = 0.
     109     DO l = llm, 1, -1
     110        sig(l) = sig(l+1) + dsig(l)
     111     ENDDO
     112  case("read")
     113     call new_unit(unit)
     114     open(unit, file="hybrid.txt", status="old", action="read", &
     115          position="rewind")
     116     read(unit, fmt=*) ! skip title line
     117     do l = 1, llm + 1
     118        read(unit, fmt=*) sig(l)
     119     end do
     120     close(unit)
     121  case default
     122     call abort_gcm("disvert", 'Wrong value for "vert_sampling"', 1)
     123  END select
    105124
    106125  DO l=1, llm

  • LMDZ5/branches/testing/libf/dyn3d/dynetat0.F

    r1421 r1665  
    119119      day_ini = tab_cntrl(30)
    120120      itau_dyn = tab_cntrl(31)
     121      start_time = tab_cntrl(32)
    121122c   .................................................................
    122123c

  • LMDZ5/branches/testing/libf/dyn3d/dynredem.F

    r1664 r1665  
    120120       tab_cntrl(30) = REAL(iday_end)
    121121       tab_cntrl(31) = REAL(itau_dyn + itaufin)
     122c start_time: start_time of simulation (not necessarily 0.)
     123       tab_cntrl(32) = start_time
    122124c
    123125c    .........................................................

  • LMDZ5/branches/testing/libf/dyn3d/etat0_netcdf.F90

    r1520 r1665  
    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

  • LMDZ5/branches/testing/libf/dyn3d/exner_hyb.F

    r1520 r1665  
    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

  • LMDZ5/branches/testing/libf/dyn3d/exner_milieu.F

    r1520 r1665  
    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

  • LMDZ5/branches/testing/libf/dyn3d/gcm.F

    r1664 r1665  
    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:
    24 #ifdef CPP_EARTH
     23#ifdef CPP_PHYS
    2524      USE dimphy
    2625      USE comgeomphy
    … …  
    8887
    8988      REAL zdtvr
    90       INTEGER nbetatmoy, nbetatdem,nbetat
    9189
    9290c   variables dynamiques
    … …  
    181179! A nettoyer. On ne veut qu'une ou deux routines d'interface
    182180! dynamique -> physique pour l'initialisation
    183 ! Ehouarn : temporarily (?) keep this only for Earth
    184       if (planet_type.eq."earth") then
    185 #ifdef CPP_EARTH
     181#ifdef CPP_PHYS
    186182      CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
    187183      call InitComgeomphy
    188184#endif
    189       endif
    190185!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    191186c-----------------------------------------------------------------------
    … …  
    305300C on remet le calendrier à zero si demande
    306301c
     302      IF (start_time /= starttime) then
     303        WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le'
     304     &,' fichier restart ne correspond pas à celle lue dans le run.def'
     305        IF (raz_date == 1) then
     306          WRITE(lunout,*)'Je prends l''heure lue dans run.def'
     307          start_time = starttime
     308        ELSE
     309          WRITE(lunout,*)'Je m''arrete'
     310          CALL abort
     311        ENDIF
     312      ENDIF
    307313      IF (raz_date == 1) THEN
    308314        annee_ref = anneeref
    … …  
    373379#endif
    374380
    375 c  nombre d'etats dans les fichiers demarrage et histoire
    376       nbetatdem = nday / iecri
    377       nbetatmoy = nday / periodav + 1
    378381
    379382      if (iflag_phys.eq.1) then
    … …  
    428431         WRITE(lunout,*)
    429432     .       'GCM: WARNING!!! vitesse verticale nulle dans la physique'
    430 ! Earth:
    431          if (planet_type.eq."earth") then
    432 #ifdef CPP_EARTH
     433! Physics:
     434#ifdef CPP_PHYS
    433435         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys ,
    434436     ,                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp     )
    435437#endif
    436          endif ! of if (planet_type.eq."earth")
    437438         call_iniphys=.false.
    438439      ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
    439 !#endif
    440440
    441441c  numero de stockage pour les fichiers de redemarrage:
    … …  
    459459#endif
    460460
    461 #ifdef CPP_EARTH
     461#ifdef CPP_PHYS
    462462! Create start file (startphy.nc) and boundary conditions (limit.nc)
    463463! for the Earth verstion

  • LMDZ5/branches/testing/libf/dyn3d/guide_mod.F90

    r1520 r1665  
    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

  • LMDZ5/branches/testing/libf/dyn3d/iniacademic.F90

    r1664 r1665  
    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)

  • LMDZ5/branches/testing/libf/dyn3d/inidissip.F90

    r1502 r1665  
    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

  • LMDZ5/branches/testing/libf/dyn3d/inigrads.F

    r524 r1665  
    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'

  • LMDZ5/branches/testing/libf/dyn3d/integrd.F

    r1550 r1665  
    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")

  • LMDZ5/branches/testing/libf/dyn3d/leapfrog.F

    r1529 r1665  
    7272
    7373      real zqmin,zqmax
    74       INTEGER nbetatmoy, nbetatdem,nbetat
    7574
    7675c   variables dynamiques
    … …  
    118117
    119118      REAL  SSUM
    120       REAL time_0 , finvmaold(ip1jmp1,llm)
     119      REAL time_0
     120!     REAL finvmaold(ip1jmp1,llm)
    121121
    122122cym      LOGICAL  lafin
    … …  
    212212      dq(:,:,:)=0.
    213213      CALL pression ( ip1jmp1, ap, bp, ps, p       )
    214       if (disvert_type==1) then
     214      if (pressure_exner) then
    215215        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
    216       else ! we assume that we are in the disvert_type==2 case
     216      else
    217217        CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
    218218      endif
    … …  
    222222c   ----------------------------------
    223223
    224    1  CONTINUE
    225 
    226       jD_cur = jD_ref + day_ini - day_ref + int (itau * dtvr / daysec)
    227       jH_cur = jH_ref +                                                 &
    228      &          (itau * dtvr / daysec - int(itau * dtvr / daysec))
     224   1  CONTINUE ! Matsuno Forward step begins here
     225
     226      jD_cur = jD_ref + day_ini - day_ref +                             &
     227     &          itau/day_step
     228      jH_cur = jH_ref + start_time +                                    &
     229     &          mod(itau,day_step)/float(day_step)
     230      jD_cur = jD_cur + int(jH_cur)
     231      jH_cur = jH_cur - int(jH_cur)
    229232
    230233
    … …  
    255258
    256259c   ...    P.Le Van .26/04/94  ....
    257 
    258       CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
    259       CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
    260 
    261    2  CONTINUE
     260! Ehouarn: finvmaold is actually not used
     261!      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
     262!      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
     263
     264   2  CONTINUE ! Matsuno backward or leapfrog step begins here
    262265
    263266c-----------------------------------------------------------------------
    … …  
    302305c   --------------------------------
    303306
     307      ! compute geopotential phi()
    304308      CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
    305309
    … …  
    315319
    316320      IF( forward. OR . leapf )  THEN
    317 
     321! Ehouarn: NB: at this point p with ps are not synchronized
     322!              (whereas mass and ps are...)
    318323         CALL caladvtrac(q,pbaru,pbarv,
    319324     *        p, masse, dq,  teta,
    … …  
    340345
    341346       CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
    342      $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
    343      $              finvmaold                                    )
     347     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
     348!     $              finvmaold                                    )
    344349
    345350
    … …  
    364369
    365370         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
    366          if (disvert_type==1) then
     371         if (pressure_exner) then
    367372           CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
    368          else ! we assume that we are in the disvert_type==2 case
     373         else
    369374           CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
    370375         endif
    … …  
    372377!           rdaym_ini  = itau * dtvr / daysec
    373378!           rdayvrai   = rdaym_ini  + day_ini
    374            jD_cur = jD_ref + day_ini - day_ref
    375      $        + int (itau * dtvr / daysec)
    376            jH_cur = jH_ref +                                            &
    377      &              (itau * dtvr / daysec - int(itau * dtvr / daysec))
     379!           jD_cur = jD_ref + day_ini - day_ref
     380!     $        + int (itau * dtvr / daysec)
     381!           jH_cur = jH_ref +                                            &
     382!     &              (itau * dtvr / daysec - int(itau * dtvr / daysec))
     383           jD_cur = jD_ref + day_ini - day_ref +                        &
     384     &          itau/day_step
     385           jH_cur = jH_ref + start_time +                               &
     386     &              mod(itau,day_step)/float(day_step)
     387           jD_cur = jD_cur + int(jH_cur)
     388           jH_cur = jH_cur - int(jH_cur)
    378389!         write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
    379390!         call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
    … …  
    394405! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
    395406           IF (planet_type.eq."earth") THEN
     407#ifdef CPP_EARTH
    396408            CALL diagedyn(ztit,2,1,1,dtphys
    397409     &    , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
     410#endif
    398411           ENDIF
    399412         ENDIF ! of IF (ip_ebil_dyn.ge.1 )
    … …  
    472485
    473486        CALL pression ( ip1jmp1, ap, bp, ps, p                  )
    474         if (disvert_type==1) then
     487        if (pressure_exner) then
    475488          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
    476         else ! we assume that we are in the disvert_type==2 case
     489        else
    477490          CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
    478491        endif
    … …  
    529542        ENDDO
    530543
    531         DO ij =  1,iim
    532           tppn(ij)  = aire(  ij    ) * ps (  ij    )
    533           tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
    534         ENDDO
    535           tpn  = SSUM(iim,tppn,1)/apoln
    536           tps  = SSUM(iim,tpps,1)/apols
    537 
    538         DO ij = 1, iip1
    539           ps(  ij    ) = tpn
    540           ps(ij+ip1jm) = tps
    541         ENDDO
    542 
     544        if (1 == 0) then
     545!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
     546!!!                     2) should probably not be here anyway
     547!!! but are kept for those who would want to revert to previous behaviour
     548           DO ij =  1,iim
     549             tppn(ij)  = aire(  ij    ) * ps (  ij    )
     550             tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
     551           ENDDO
     552             tpn  = SSUM(iim,tppn,1)/apoln
     553             tps  = SSUM(iim,tpps,1)/apols
     554
     555           DO ij = 1, iip1
     556             ps(  ij    ) = tpn
     557             ps(ij+ip1jm) = tps
     558           ENDDO
     559        endif ! of if (1 == 0)
    543560
    544561      END IF ! of IF(apdiss)
    … …  
    622639             ! Ehouarn: output only during LF or Backward Matsuno
    623640             if (leapf.or.(.not.leapf.and.(.not.forward))) then
    624               nbetat = nbetatdem
    625641              CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
    626642              unat=0.
    … …  
    652668!              if (planet_type.eq."earth") then
    653669! Write an Earth-format restart file
    654                 CALL dynredem1("restart.nc",0.0,
     670                CALL dynredem1("restart.nc",start_time,
    655671     &                         vcov,ucov,teta,q,masse,ps)
    656672!              endif ! of if (planet_type.eq."earth")
    657673
    658674              CLOSE(99)
     675              !!! Ehouarn: Why not stop here and now?
    659676            ENDIF ! of IF (itau.EQ.itaufin)
    660677
    … …  
    740757              IF(MOD(itau,iecri         ).EQ.0) THEN
    741758c              IF(MOD(itau,iecri*day_step).EQ.0) THEN
    742                 nbetat = nbetatdem
    743759                CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
    744760                unat=0.
    … …  
    763779              IF(itau.EQ.itaufin) THEN
    764780!                if (planet_type.eq."earth") then
    765                   CALL dynredem1("restart.nc",0.0,
     781                  CALL dynredem1("restart.nc",start_time,
    766782     &                           vcov,ucov,teta,q,masse,ps)
    767783!                endif ! of if (planet_type.eq."earth")

  • LMDZ5/branches/testing/libf/dyn3d/temps.h

    r1279 r1665  
    1414
    1515      COMMON/temps/itaufin, dt, day_ini, day_end, annee_ref, day_ref,   &
    16      &             itau_dyn, itau_phy, jD_ref, jH_ref, calend
     16     &             itau_dyn, itau_phy, jD_ref, jH_ref, calend,          &
     17     &             start_time
     18
    1719
    1820      INTEGER   itaufin
    1921      INTEGER itau_dyn, itau_phy
    2022      INTEGER day_ini, day_end, annee_ref, day_ref
    21       REAL      dt, jD_ref, jH_ref
     23      REAL      dt, jD_ref, jH_ref, start_time
    2224      CHARACTER (len=10) :: calend
    2325

  • LMDZ5/branches/testing/libf/dyn3d/wrgrads.F

    r1025 r1665  
    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'

  • LMDZ5/branches/testing/libf/dyn3dpar/bands.F90

    r1279 r1665  
    9393   SUBROUTINE  Set_Bands
    9494     USE parallel
    95 #ifdef CPP_EARTH
    96 ! Ehouarn: what follows is only related to // physics; for now only for Earth
     95#ifdef CPP_PHYS
     96! Ehouarn: what follows is only related to // physics
    9797     USE mod_phys_lmdz_para, ONLY : jj_para_begin,jj_para_end
    9898#endif
    … …  
    106106      enddo
    107107         
    108 #ifdef CPP_EARTH
    109 ! Ehouarn: what follows is only related to // physics; for now only for Earth         
     108#ifdef CPP_PHYS
    110109      do i=0,MPI_Size-1
    111110        jj_Nb_physic(i)=jj_para_end(i)-jj_para_begin(i)+1
    … …  
    332331    subroutine AdjustBands_physic
    333332      use times
    334 #ifdef CPP_EARTH
    335 ! Ehouarn: what follows is only related to // physics; for now only for Earth
     333#ifdef CPP_PHYS
     334! Ehouarn: what follows is only related to // physics
    336335      USE mod_phys_lmdz_para, only : klon_mpi_para_nb
    337336#endif
    … …  
    359358      medium=medium/mpi_size     
    360359      NbTot=0
    361 #ifdef CPP_EARTH
    362 ! Ehouarn: what follows is only related to // physics; for now only for Earth
     360#ifdef CPP_PHYS
    363361      do i=0,mpi_size-1
    364362        Inc(i)=nint(klon_mpi_para_nb(i)*(medium-value(i))/value(i))

  • LMDZ5/branches/testing/libf/dyn3dpar/calfis_p.F

    r1407 r1665  
    2727     $                  pdqfi,
    2828     $                  pdpsfi)
    29 #ifdef CPP_EARTH
    30 ! Ehouarn: For now, calfis_p needs Earth physics
    31 c
    32 c    Auteur :  P. Le Van, F. Hourdin
    33 c   .........
     29#ifdef CPP_PHYS
     30! If using physics
    3431      USE dimphy
    3532      USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root
    … …  
    146143      REAL clesphy0( longcles )
    147144
    148 #ifdef CPP_EARTH
     145#ifdef CPP_PHYS
     146! Ehouarn: for now calfis_p needs some informations from physics to compile
    149147c    Local variables :
    150148c    -----------------
    … …  
    222220      PARAMETER(ntetaSTD=3)
    223221      REAL rtetaSTD(ntetaSTD)
    224       DATA rtetaSTD/350., 380., 405./
     222      DATA rtetaSTD/350., 380., 405./ ! Earth-specific values, beware !!
    225223      REAL PVteta(klon,ntetaSTD)
    226224     
    … …  
    489487
    490488
    491       IF (is_sequential) THEN
    492 c
     489      IF (is_sequential.and.(planet_type=="earth")) THEN
     490#ifdef CPP_PHYS
     491! PVtheta calls tetalevel, which is in the physics
    493492cIM calcul PV a teta=350, 380, 405K
    494493        CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
    … …  
    496495     $           ntetaSTD,rtetaSTD,PVteta)
    497496c
     497#endif
    498498      ENDIF
    499499
    … …  
    627627c$OMP BARRIER
    628628     
    629       if (planet_type=="earth") then
    630 #ifdef CPP_EARTH
    631 
    632629!$OMP MASTER
    633630!      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
    … …  
    639636      zdqfic_omp(:,:,:)=0.
    640637
     638      if (planet_type=="earth") then
     639#ifdef CPP_PHYS
    641640      do isplit=1,nsplit_phys
    642641
    … …  
    687686      enddo
    688687
     688#endif
     689! of #ifdef CPP_PHYS
     690      endif !of if (planet_type=="earth")
     691
    689692      zdufi_omp(:,:)=zdufic_omp(:,:)/nsplit_phys
    690693      zdvfi_omp(:,:)=zdvfic_omp(:,:)/nsplit_phys
    691694      zdtfi_omp(:,:)=zdtfic_omp(:,:)/nsplit_phys
    692695      zdqfi_omp(:,:,:)=zdqfic_omp(:,:,:)/nsplit_phys
    693 
    694 #endif
    695       endif !of if (planet_type=="earth")
    696696c$OMP BARRIER
    697697
    … …  
    11101110      firstcal = .FALSE.
    11111111
    1112 #else
    1113       write(lunout,*)
    1114      & "calfis_p: for now can only work with parallel physics"
    1115       stop
    1116 #endif
    1117 ! of #ifdef CPP_EARTH
     1112#else 
     1113      write(lunout,*) 
     1114     & "calfis_p: for now can only work with parallel physics" 
     1115      stop 
     1116#endif 
     1117! of #ifdef CPP_PHYS
    11181118      RETURN
    11191119      END

  • LMDZ5/branches/testing/libf/dyn3dpar/ce0l.F90

    r1664 r1665  
    1919  USE dimphy
    2020  USE comgeomphy
    21   USE mod_phys_lmdz_para
     21  USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root
    2222  USE mod_const_mpi
    2323  USE infotrac
    … …  
    3131  IMPLICIT NONE
    3232#ifndef CPP_EARTH
     33#include "iniprint.h"
    3334  WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
    3435#else
    … …  
    4243#include "temps.h"
    4344#include "logic.h"
     45#ifdef CPP_MPI
     46      include 'mpif.h'
     47#endif
     48
    4449  INTEGER, PARAMETER            :: longcles=20
     50  INTEGER                       :: ierr
    4551  REAL,    DIMENSION(longcles)  :: clesphy0
    4652  REAL,    DIMENSION(iip1,jjp1) :: masque
    … …  
    5056  CALL conf_gcm( 99, .TRUE. , clesphy0 )
    5157
     58#ifdef CPP_MPI
    5259  CALL init_mpi
     60#endif
    5361
    5462  CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
    … …  
    115123  END IF
    116124 
     125#ifdef CPP_MPI
    117126!$OMP MASTER
    118   CALL finalize_parallel
     127  CALL MPI_FINALIZE(ierr)
     128  IF (ierr /= 0) CALL abort_gcm('ce0l','Error in MPI_FINALIZE',1)
    119129!$OMP END MASTER
     130#endif
    120131
    121132#endif

  • LMDZ5/branches/testing/libf/dyn3dpar/comvert.h

    r1520 r1665  
    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 !-----------------------------------------------------------------------

  • LMDZ5/branches/testing/libf/dyn3dpar/conf_gcm.F

    r1664 r1665  
    167167      nday = 10
    168168      CALL getin('nday',nday)
     169
     170!Config  Key  = starttime
     171!Config  Desc = Heure de depart de la simulation
     172!Config  Def  = 0
     173!Config  Help = Heure de depart de la simulation
     174!Config         en jour
     175      starttime = 0
     176      CALL getin('starttime',starttime)
    169177
    170178!Config  Key  = day_step

  • LMDZ5/branches/testing/libf/dyn3dpar/control_mod.F90

    r1502 r1665  
    1010  IMPLICIT NONE
    1111
    12   REAL    :: periodav
     12  REAL    :: periodav, starttime
    1313  INTEGER :: nday,day_step,iperiod,iapp_tracvl,nsplit_phys
    1414  INTEGER :: iconser,iecri,dissip_period,iphysiq,iecrimoy

  • LMDZ5/branches/testing/libf/dyn3dpar/disvert.F90

    r1524 r1665  
    44
    55  ! Auteur : P. Le Van
     6
     7  use new_unit_m, only: new_unit
     8  use ioipsl, only: getin
    69
    710  IMPLICIT NONE
    … …  
    2629  real zk, zkm1, dzk1, dzk2, k0, k1
    2730
    28   INTEGER l
     31  INTEGER l, unit
    2932  REAL dsigmin
    3033  REAL alpha, beta, deltaz
    31   INTEGER iostat
    3234  REAL x
    3335  character(len=*),parameter :: modname="disvert"
    3436
     37  character(len=6):: vert_sampling
     38  ! (allowed values are "param", "tropo", "strato" and "read")
     39
    3540  !-----------------------------------------------------------------------
     41
     42  print *, "Call sequence information: disvert"
    3643
    3744  ! default scaleheight is 8km for earth
    3845  scaleheight=8.
    3946
    40   OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
     47  vert_sampling = merge("strato", "tropo ", ok_strato) ! default value
     48  call getin('vert_sampling', vert_sampling)
     49  print *, 'vert_sampling = ' // vert_sampling
    4150
    42   IF (iostat == 0) THEN
    43      ! cas 1 on lit les options dans sigma.def:
     51  select case (vert_sampling)
     52  case ("param")
     53     ! On lit les options dans sigma.def:
     54     OPEN(99, file='sigma.def', status='old', form='formatted')
    4455     READ(99, *) scaleheight ! hauteur d'echelle 8.
    4556     READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
    … …  
    6879
    6980     sig(llm+1)=0.
    70 
    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
    88 
     81  case("tropo")
    8982     DO l = 1, llm
    9083        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
     84        dsig(l) = 1.0 + 7.0 * SIN(x)**2
    9885     ENDDO
    9986     dsig = dsig / sum(dsig)
    … …  
    10289        sig(l) = sig(l+1) + dsig(l)
    10390     ENDDO
    104   ENDIF
     91  case("strato")
     92     if (llm==39) then
     93        dsigmin=0.3
     94     else if (llm==50) then
     95        dsigmin=1.
     96     else
     97        write(lunout,*) trim(modname), ' ATTENTION discretisation z a ajuster'
     98        dsigmin=1.
     99     endif
     100     WRITE(LUNOUT,*) trim(modname), 'Discretisation verticale DSIGMIN=',dsigmin
     101
     102     DO l = 1, llm
     103        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
     104        dsig(l) =(dsigmin + 7. * SIN(x)**2) &
     105             *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
     106     ENDDO
     107     dsig = dsig / sum(dsig)
     108     sig(llm+1) = 0.
     109     DO l = llm, 1, -1
     110        sig(l) = sig(l+1) + dsig(l)
     111     ENDDO
     112  case("read")
     113     call new_unit(unit)
     114     open(unit, file="hybrid.txt", status="old", action="read", &
     115          position="rewind")
     116     read(unit, fmt=*) ! skip title line
     117     do l = 1, llm + 1
     118        read(unit, fmt=*) sig(l)
     119     end do
     120     close(unit)
     121  case default
     122     call abort_gcm("disvert", 'Wrong value for "vert_sampling"', 1)
     123  END select
    105124
    106125  DO l=1, llm

  • LMDZ5/branches/testing/libf/dyn3dpar/dynetat0.F

    r1421 r1665  
    119119      day_ini = tab_cntrl(30)
    120120      itau_dyn = tab_cntrl(31)
     121      start_time = tab_cntrl(32)
    121122c   .................................................................
    122123c

  • LMDZ5/branches/testing/libf/dyn3dpar/dynredem.F

    r1664 r1665  
    120120       tab_cntrl(30) = REAL(iday_end)
    121121       tab_cntrl(31) = REAL(itau_dyn + itaufin)
     122c start_time: start_time of simulation (not necessarily 0.)
     123       tab_cntrl(32) = start_time
    122124c
    123125c    .........................................................

  • LMDZ5/branches/testing/libf/dyn3dpar/dynredem_p.F

    r1664 r1665  
    120120       tab_cntrl(30) =  REAL(iday_end)
    121121       tab_cntrl(31) =  REAL(itau_dyn + itaufin)
     122c start_time: start_time of simulation (not necessarily 0.)
     123       tab_cntrl(32) = start_time
    122124c
    123125c    .........................................................

  • LMDZ5/branches/testing/libf/dyn3dpar/etat0_netcdf.F90

    r1520 r1665  
    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

  • LMDZ5/branches/testing/libf/dyn3dpar/exner_hyb.F

    r1520 r1665  
    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

  • LMDZ5/branches/testing/libf/dyn3dpar/exner_hyb_p.F

    r1664 r1665  
    6060      ! Sanity check
    6161      if (firstcall) then
    62         ! check that vertical discretization is compatible
    63         ! with this routine
    64         if (disvert_type.ne.1) then
    65           call abort_gcm(modname,
    66      &     "this routine should only be called if disvert_type==1",42)
    67         endif
    68        
    6962        ! sanity checks for Shallow Water case (1 vertical layer)
    7063        if (llm.eq.1) then

  • LMDZ5/branches/testing/libf/dyn3dpar/exner_milieu.F

    r1520 r1665  
    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

  • LMDZ5/branches/testing/libf/dyn3dpar/exner_milieu_p.F

    r1664 r1665  
    5656      ! Sanity check
    5757      if (firstcall) then
    58         ! check that vertical discretization is compatible
    59         ! with this routine
    60         if (disvert_type.ne.2) then
    61           call abort_gcm(modname,
    62      &     "this routine should only be called if disvert_type==2",42)
    63         endif
    64        
    6558        ! sanity checks for Shallow Water case (1 vertical layer)
    6659        if (llm.eq.1) then

  • LMDZ5/branches/testing/libf/dyn3dpar/filtreg_p.F

    r1146 r1665  
    208208               IF( ifiltre.EQ.-2 )   THEN
    209209                  DO j = jdfil,jffil
     210#ifdef BLAS
    210211                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
    211212     &                    matrinvn(1,1,j), iim,
    212213     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
    213214     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
     215#else
     216                     champ_fft(:,j-jdfil+1,:)
     217     &                    =matmul(matrinvn(:,:,j),champ_loc(:iim,j,:))
     218#endif
    214219                  ENDDO
    215220                 
    216221               ELSE IF ( griscal )     THEN
    217222                  DO j = jdfil,jffil
     223#ifdef BLAS
    218224                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
    219225     &                    matriceun(1,1,j), iim,
    220226     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
    221227     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
     228#else
     229                     champ_fft(:,j-jdfil+1,:)
     230     &                    =matmul(matriceun(:,:,j),champ_loc(:iim,j,:))
     231#endif
    222232                  ENDDO
    223233                 
    224234               ELSE
    225235                  DO j = jdfil,jffil
     236#ifdef BLAS
    226237                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
    227238     &                    matricevn(1,1,j), iim,
    228239     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
    229240     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
     241#else
     242                     champ_fft(:,j-jdfil+1,:)
     243     &                    =matmul(matricevn(:,:,j),champ_loc(:iim,j,:))
     244#endif
    230245                  ENDDO
    231246                 
    … …  
    236251               IF( ifiltre.EQ.-2 )   THEN
    237252                  DO j = jdfil,jffil
     253#ifdef BLAS
    238254                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
    239255     &                    matrinvs(1,1,j-jfiltsu+1), iim,
    240256     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
    241257     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
     258#else
     259                     champ_fft(:,j-jdfil+1,:)
     260     &                    =matmul(matrinvs(:,:,j-jfiltsu+1),
     261     &                            champ_loc(:iim,j,:))
     262#endif
    242263                  ENDDO
    243264                 
    … …  
    245266                 
    246267                  DO j = jdfil,jffil
     268#ifdef BLAS
    247269                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
    248270     &                    matriceus(1,1,j-jfiltsu+1), iim,
    249271     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
    250272     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
     273#else
     274                     champ_fft(:,j-jdfil+1,:)
     275     &                    =matmul(matriceus(:,:,j-jfiltsu+1),
     276     &                            champ_loc(:iim,j,:))
     277#endif
    251278                  ENDDO
    252279                 
    … …  
    254281                 
    255282                  DO j = jdfil,jffil
     283#ifdef BLAS
    256284                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
    257285     &                    matricevs(1,1,j-jfiltsv+1), iim,
    258286     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
    259287     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
     288#else
     289                     champ_fft(:,j-jdfil+1,:)
     290     &                    =matmul(matricevs(:,:,j-jfiltsv+1),
     291     &                            champ_loc(:iim,j,:))
     292#endif
    260293                  ENDDO
    261294                 

  • LMDZ5/branches/testing/libf/dyn3dpar/gcm.F

    r1664 r1665  
    2020      USE control_mod
    2121
    22 ! Ehouarn: for now these only apply to Earth:
    23 #ifdef CPP_EARTH
     22#ifdef CPP_PHYS
    2423      USE mod_grid_phy_lmdz
    2524      USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
    … …  
    8786
    8887      REAL zdtvr
    89 c      INTEGER nbetatmoy, nbetatdem,nbetat
    90       INTEGER nbetatmoy, nbetatdem
    9188
    9289c   variables dynamiques
    … …  
    189186      call ini_getparam("out.def")
    190187      call Read_Distrib
    191 ! Ehouarn : temporarily (?) keep this only for Earth
    192       if (planet_type.eq."earth") then
    193 #ifdef CPP_EARTH
     188
     189#ifdef CPP_PHYS
    194190        CALL Init_Phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys)
    195191#endif
    196       endif ! of if (planet_type.eq."earth")
    197192      CALL set_bands
    198 #ifdef CPP_EARTH
    199 ! Ehouarn: For now only Earth physics is parallel
     193#ifdef CPP_PHYS
    200194      CALL Init_interface_dyn_phys
    201195#endif
    … …  
    209203c$OMP END PARALLEL
    210204
    211 ! Ehouarn : temporarily (?) keep this only for Earth
    212       if (planet_type.eq."earth") then
    213 #ifdef CPP_EARTH
     205#ifdef CPP_PHYS
    214206c$OMP PARALLEL
    215207      call InitComgeomphy
    216208c$OMP END PARALLEL
    217209#endif
    218       endif ! of if (planet_type.eq."earth")
    219210
    220211c-----------------------------------------------------------------------
    … …  
    323314C on remet le calendrier à zero si demande
    324315c
     316      IF (start_time /= starttime) then
     317        WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le'
     318     &,' fichier restart ne correspond pas à celle lue dans le run.def'
     319        IF (raz_date == 1) then
     320          WRITE(lunout,*)'Je prends l''heure lue dans run.def'
     321          start_time = starttime
     322        ELSE
     323          WRITE(lunout,*)'Je m''arrete'
     324          CALL abort
     325        ENDIF
     326      ENDIF
    325327      IF (raz_date == 1) THEN
    326328        annee_ref = anneeref
    … …  
    390392#endif
    391393
    392 c  nombre d'etats dans les fichiers demarrage et histoire
    393       nbetatdem = nday / iecri
    394       nbetatmoy = nday / periodav + 1
    395394
    396395      if (iflag_phys.eq.1) then
    … …  
    445444         WRITE(lunout,*)
    446445     .       'GCM: WARNING!!! vitesse verticale nulle dans la physique'
    447 ! Earth:
    448          if (planet_type.eq."earth") then
    449 #ifdef CPP_EARTH
     446! Physics:
     447#ifdef CPP_PHYS
    450448         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys ,
    451449     ,                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp     )
    452450#endif
    453          endif ! of if (planet_type.eq."earth")
    454451         call_iniphys=.false.
    455452      ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
    … …  
    484481 301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
    485482 302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
     483#endif
     484
     485#ifdef CPP_PHYS
     486! Create start file (startphy.nc) and boundary conditions (limit.nc)
     487! for the Earth verstion
     488       if (iflag_phys>=100) then
     489          call iniaqua(ngridmx,latfi,lonfi,iflag_phys)
     490       endif
    486491#endif
    487492

  • LMDZ5/branches/testing/libf/dyn3dpar/gr_dyn_fi_p.F

    r1279 r1665  
    33!
    44      SUBROUTINE gr_dyn_fi_p(nfield,im,jm,ngrid,pdyn,pfi)
    5 #ifdef CPP_EARTH
     5#ifdef CPP_PHYS
    66! Interface with parallel physics,
    7 ! for now this routine only works with Earth physics
    87      USE mod_interface_dyn_phys
    98      USE dimphy
    … …  
    4039      ENDDO
    4140c$OMP END DO NOWAIT
    42 #else
    43       write(lunout,*) "gr_fi_dyn_p : This routine should not be called",
    44      &   "without parallelized physics"
    45       stop
    4641#endif
    47 ! of #ifdef CPP_EARTH
     42! of #ifdef CPP_PHYS
    4843      RETURN
    4944      END

  • LMDZ5/branches/testing/libf/dyn3dpar/gr_fi_dyn_p.F

    r1279 r1665  
    33!
    44      SUBROUTINE gr_fi_dyn_p(nfield,ngrid,im,jm,pfi,pdyn)
    5 #ifdef CPP_EARTH
     5#ifdef CPP_PHYS
    66! Interface with parallel physics,
    7 ! for now this routine only works with Earth physics
    87      USE mod_interface_dyn_phys
    98      USE dimphy
    … …  
    5251      ENDDO
    5352c$OMP END DO NOWAIT
    54 #else
    55       write(lunout,*) "gr_fi_dyn_p : This routine should not be called",
    56      &   "without parallelized physics"
    57       stop
    5853#endif
    59 ! of #ifdef CPP_EARTH
     54! of #ifdef CPP_PHYS
    6055      RETURN
    6156      END

  • LMDZ5/branches/testing/libf/dyn3dpar/guide_p_mod.F90

    r1520 r1665  
    455455!       Calcul niveaux pression milieu de couches
    456456        CALL pression_p( ip1jmp1, ap, bp, ps, p )
    457         if (disvert_type==1) then
     457        if (pressure_exner) then
    458458          CALL exner_hyb_p(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
    459459        else
    … …  
    755755    ELSE
    756756        CALL pression_p( ip1jmp1, ap, bp, psi, p )
    757         if (disvert_type==1) then
     757        if (pressure_exner) then
    758758          CALL exner_hyb_p(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
    759         else ! we assume that we are in the disvert_type==2 case
     759        else
    760760          CALL exner_milieu_p(ip1jmp1,psi,p,beta,pks,pk,pkf)
    761761        endif

  • LMDZ5/branches/testing/libf/dyn3dpar/iniacademic.F90

    r1664 r1665  
    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)

  • LMDZ5/branches/testing/libf/dyn3dpar/inidissip.F90

    r1502 r1665  
    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

  • LMDZ5/branches/testing/libf/dyn3dpar/inigrads.F

    r774 r1665  
    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'

  • LMDZ5/branches/testing/libf/dyn3dpar/integrd_p.F

    r1550 r1665  
    44      SUBROUTINE integrd_p
    55     $  (  nq,vcovm1,ucovm1,tetam1,psm1,massem1,
    6      $     dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps0,masse,phis,finvmaold)
     6     $     dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps0,masse,phis) !,finvmaold)
    77      USE parallel
    88      USE control_mod, only : planet_type
    … …  
    3333#include "temps.h"
    3434#include "serre.h"
     35#include "iniprint.h"
    3536
    3637c   Arguments:
    3738c   ----------
    3839
    39       INTEGER nq
    40 
    41       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm)
    42       REAL q(ip1jmp1,llm,nq)
    43       REAL ps0(ip1jmp1),masse(ip1jmp1,llm),phis(ip1jmp1)
    44 
    45       REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
    46       REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1),massem1(ip1jmp1,llm)
    47 
    48       REAL dv(ip1jm,llm),du(ip1jmp1,llm)
    49       REAL dteta(ip1jmp1,llm),dp(ip1jmp1)
    50       REAL dq(ip1jmp1,llm,nq), finvmaold(ip1jmp1,llm)
     40      integer,intent(in) :: nq ! number of tracers to handle in this routine
     41      real,intent(inout) :: vcov(ip1jm,llm) ! covariant meridional wind
     42      real,intent(inout) :: ucov(ip1jmp1,llm) ! covariant zonal wind
     43      real,intent(inout) :: teta(ip1jmp1,llm) ! potential temperature
     44      real,intent(inout) :: q(ip1jmp1,llm,nq) ! advected tracers
     45      real,intent(inout) :: ps0(ip1jmp1) ! surface pressure
     46      real,intent(inout) :: masse(ip1jmp1,llm) ! atmospheric mass
     47      real,intent(in) :: phis(ip1jmp1) ! ground geopotential !!! unused
     48      ! values at previous time step
     49      real,intent(inout) :: vcovm1(ip1jm,llm)
     50      real,intent(inout) :: ucovm1(ip1jmp1,llm)
     51      real,intent(inout) :: tetam1(ip1jmp1,llm)
     52      real,intent(inout) :: psm1(ip1jmp1)
     53      real,intent(inout) :: massem1(ip1jmp1,llm)
     54      ! the tendencies to add
     55      real,intent(in) :: dv(ip1jm,llm)
     56      real,intent(in) :: du(ip1jmp1,llm)
     57      real,intent(in) :: dteta(ip1jmp1,llm)
     58      real,intent(in) :: dp(ip1jmp1)
     59      real,intent(in) :: dq(ip1jmp1,llm,nq) !!! unused
     60!      real,intent(out) :: finvmaold(ip1jmp1,llm) !!! unused
    5161
    5262c   Local:
    … …  
    5464
    5565      REAL vscr( ip1jm ),uscr( ip1jmp1 ),hscr( ip1jmp1 ),pscr(ip1jmp1)
    56       REAL massescr( ip1jmp1,llm ), finvmasse(ip1jmp1,llm)
     66      REAL massescr( ip1jmp1,llm )
     67!      REAL finvmasse(ip1jmp1,llm)
    5768      REAL,SAVE :: p(ip1jmp1,llmp1)
    5869      REAL tpn,tps,tppn(iim),tpps(iim)
    … …  
    6071      REAL,SAVE :: deltap( ip1jmp1,llm )
    6172
    62       INTEGER  l,ij,iq
     73      INTEGER  l,ij,iq,i,j
    6374
    6475      REAL SSUM
    … …  
    126137       
    127138        IF( .NOT. checksum ) THEN
    128          PRINT*,' Au point ij = ',stop_it, ' , pression sol neg. '
    129      &         , ps(stop_it)
    130          print *, ' dans integrd'
    131          stop 1
     139         write(lunout,*) "integrd: negative surface pressure ",
     140     &                                                ps(stop_it)
     141         write(lunout,*) " at node ij =", stop_it
     142         ! since ij=j+(i-1)*jjp1 , we have
     143         j=modulo(stop_it,jjp1)
     144         i=1+(stop_it-j)/jjp1
     145         write(lunout,*) " lon = ",rlonv(i)*180./pi, " deg",
     146     &                   " lat = ",rlatu(j)*180./pi, " deg"
    132147        ENDIF
    133148
    … …  
    167182      CALL massdair_p (     p  , masse         )
    168183
    169 c      CALL   SCOPY( ijp1llm  , masse, 1, finvmasse,  1      )
    170       ijb=ij_begin
    171       ije=ij_end
    172      
    173 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    174       DO  l = 1,llm
    175         finvmasse(ijb:ije,l)=masse(ijb:ije,l)
    176       ENDDO
    177 c$OMP END DO NOWAIT
    178 
    179       jjb=jj_begin
    180       jje=jj_end
    181       CALL filtreg_p( finvmasse,jjb,jje, jjp1, llm, -2, 2, .TRUE., 1  )
     184! Ehouarn : we don't use/need finvmaold and finvmasse,
     185!           so might as well not compute them
     186!c      CALL   SCOPY( ijp1llm  , masse, 1, finvmasse,  1      )
     187!      ijb=ij_begin
     188!      ije=ij_end
     189!     
     190!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     191!      DO  l = 1,llm
     192!        finvmasse(ijb:ije,l)=masse(ijb:ije,l)
     193!      ENDDO
     194!c$OMP END DO NOWAIT
     195!
     196!      jjb=jj_begin
     197!      jje=jj_end
     198!      CALL filtreg_p( finvmasse,jjb,jje, jjp1, llm, -2, 2, .TRUE., 1  )
    182199c
    183200
    … …  
    330347      ENDIF
    331348     
    332 c         CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
    333 
    334 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    335       DO l = 1, llm     
    336         finvmaold(ijb:ije,l)=finvmasse(ijb:ije,l)       
    337       ENDDO
    338 c$OMP END DO NOWAIT
     349! Ehouarn: forget about finvmaold
     350!c         CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
     351!
     352!c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     353!      DO l = 1, llm     
     354!        finvmaold(ijb:ije,l)=finvmasse(ijb:ije,l)       
     355!      ENDDO
     356!c$OMP END DO NOWAIT
    339357
    340358      endif ! of if (planet_type.eq."earth")

  • LMDZ5/branches/testing/libf/dyn3dpar/leapfrog_p.F

    r1664 r1665  
    7575
    7676      real zqmin,zqmax
    77       INTEGER nbetatmoy, nbetatdem,nbetat
    7877
    7978c   variables dynamiques
    … …  
    125124      REAL  SSUM
    126125      REAL time_0
    127       REAL,SAVE :: finvmaold(ip1jmp1,llm)
     126!      REAL,SAVE :: finvmaold(ip1jmp1,llm)
    128127
    129128cym      LOGICAL  lafin
    … …  
    234233      dq(:,:,:)=0.
    235234      CALL pression ( ip1jmp1, ap, bp, ps, p       )
    236       if (disvert_type==1) then
     235      if (pressure_exner) then
    237236        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
    238       else ! we assume that we are in the disvert_type==2 case
     237      else
    239238        CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
    240239      endif
    … …  
    245244c et du parallelisme !!
    246245
    247    1  CONTINUE
    248 
    249       jD_cur = jD_ref + day_ini - day_ref + int (itau * dtvr / daysec)
    250       jH_cur = jH_ref +                                                 &
    251      &          (itau * dtvr / daysec - int(itau * dtvr / daysec))
     246   1  CONTINUE ! Matsuno Forward step begins here
     247
     248      jD_cur = jD_ref + day_ini - day_ref +                             &
     249     &          itau/day_step
     250      jH_cur = jH_ref + start_time +                                    &
     251     &         mod(itau,day_step)/float(day_step)
     252      if (jH_cur > 1.0 ) then
     253        jD_cur = jD_cur +1.
     254        jH_cur = jH_cur -1.
     255      endif
    252256
    253257
    … …  
    280284         massem1= masse
    281285         psm1= ps
    282          
    283          finvmaold = masse
    284          CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
     286
     287! Ehouarn: finvmaold is actually not used       
     288!         finvmaold = masse
     289!         CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
    285290c$OMP END MASTER
    286291c$OMP BARRIER
    … …  
    300305           tetam1   (ijb:ije,l) = teta  (ijb:ije,l)
    301306           massem1  (ijb:ije,l) = masse (ijb:ije,l)
    302            finvmaold(ijb:ije,l)=masse(ijb:ije,l)
     307!           finvmaold(ijb:ije,l)=masse(ijb:ije,l)
    303308                 
    304309           if (pole_sud) ije=ij_end-iip1
    … …  
    309314c$OMP ENDDO 
    310315
    311 
    312           CALL filtreg_p ( finvmaold ,jj_begin,jj_end,jjp1,
    313      .                    llm, -2,2, .TRUE., 1 )
     316! Ehouarn: finvmaold not used
     317!          CALL filtreg_p ( finvmaold ,jj_begin,jj_end,jjp1,
     318!     .                    llm, -2,2, .TRUE., 1 )
    314319
    315320       endif ! of if (FirstCaldyn)
    … …  
    327332cym      call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
    328333
    329    2  CONTINUE
     334   2  CONTINUE ! Matsuno backward or leapfrog step begins here
    330335
    331336c$OMP MASTER
    … …  
    472477         call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
    473478     &                                jj_Nb_caldyn,0,0,TestRequest)
    474          call Register_SwapFieldHallo(finvmaold,finvmaold,ip1jmp1,llm,
    475      &                                jj_Nb_caldyn,0,0,TestRequest)
     479!         call Register_SwapFieldHallo(finvmaold,finvmaold,ip1jmp1,llm,
     480!     &                                jj_Nb_caldyn,0,0,TestRequest)
    476481 
    477482        do j=1,nqtot
    … …  
    555560      call start_timer(timer_caldyn)
    556561
     562      ! compute geopotential phi()
    557563      CALL geopot_p  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
    558564
    … …  
    629635
    630636       CALL integrd_p ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
    631      $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
    632      $              finvmaold                                    )
     637     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
     638!     $              finvmaold                                    )
    633639
    634640!       CALL FTRACE_REGION_END("integrd")
    … …  
    693699
    694700c$OMP BARRIER
    695          if (disvert_type==1) then
     701         if (pressure_exner) then
    696702           CALL exner_hyb_p(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
    697          else ! we assume that we are in the disvert_type==2 case
     703         else
    698704           CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
    699705         endif
    700706c$OMP BARRIER
    701707           jD_cur = jD_ref + day_ini - day_ref
    702      $        + int (itau * dtvr / daysec)
    703            jH_cur = jH_ref +                                            &
    704      &              (itau * dtvr / daysec - int(itau * dtvr / daysec))
     708     $        + itau/day_step
     709           jH_cur = jH_ref + start_time +                                &
     710     &              mod(itau,day_step)/float(day_step)
    705711!         call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
     712           if (jH_cur > 1.0 ) then
     713             jD_cur = jD_cur +1.
     714             jH_cur = jH_cur -1.
     715           endif
    706716
    707717c rajout debug
    … …  
    719729! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
    720730           IF (planet_type.eq."earth") THEN
     731#ifdef CPP_EARTH
    721732            CALL diagedyn(ztit,2,1,1,dtphys
    722733     &    , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
     734#endif
    723735           ENDIF
    724736      ENDIF
    … …  
    10361048        CALL pression_p ( ip1jmp1, ap, bp, ps, p                  )
    10371049c$OMP BARRIER
    1038         if (disvert_type==1) then
     1050        if (pressure_exner) then
    10391051          CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
    1040         else ! we assume that we are in the disvert_type==2 case
     1052        else
    10411053          CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
    10421054        endif
    … …  
    11851197c$OMP END DO NOWAIT
    11861198
     1199         if (1 == 0) then
     1200!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
     1201!!!                     2) should probably not be here anyway
     1202!!! but are kept for those who would want to revert to previous behaviour
    11871203c$OMP MASTER               
    11881204          DO ij =  1,iim
    … …  
    11951211          ENDDO
    11961212c$OMP END MASTER
    1197         endif
     1213         endif ! of if (1 == 0)
     1214        endif ! of of (pole_nord)
    11981215       
    11991216        if (pole_sud) then
    … …  
    12111228c$OMP END DO NOWAIT
    12121229
     1230         if (1 == 0) then
     1231!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
     1232!!!                     2) should probably not be here anyway
     1233!!! but are kept for those who would want to revert to previous behaviour
    12131234c$OMP MASTER               
    12141235          DO ij =  1,iim
    … …  
    12211242          ENDDO
    12221243c$OMP END MASTER
    1223         endif
     1244         endif ! of if (1 == 0)
     1245        endif ! of if (pole_sud)
    12241246
    12251247
    … …  
    14311453c$OMP BARRIER
    14321454c$OMP MASTER
    1433               nbetat = nbetatdem
    14341455              CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
    14351456       
    … …  
    16341655c$OMP BARRIER
    16351656c$OMP MASTER
    1636                 nbetat = nbetatdem
    16371657                CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
    16381658

  • LMDZ5/branches/testing/libf/dyn3dpar/mod_interface_dyn_phys.F90

    r1279 r1665  
    77 
    88 
    9 #ifdef CPP_EARTH
     9#ifdef CPP_PHYS
    1010! Interface with parallel physics,
    11 ! for now this routine only works with Earth physics
    1211CONTAINS
    1312 
    … …  
    5655  END SUBROUTINE Init_interface_dyn_phys
    5756#endif
    58 ! of #ifdef CPP_EARTH
     57! of #ifdef CPP_PHYS
    5958END MODULE mod_interface_dyn_phys

  • LMDZ5/branches/testing/libf/dyn3dpar/temps.h

    r1279 r1665  
    1414
    1515      COMMON/temps/itaufin, dt, day_ini, day_end, annee_ref, day_ref,   &
    16      &             itau_dyn, itau_phy, jD_ref, jH_ref, calend
     16     &             itau_dyn, itau_phy, jD_ref, jH_ref, calend,          &
     17     &             start_time
     18
    1719
    1820      INTEGER   itaufin
    1921      INTEGER itau_dyn, itau_phy
    2022      INTEGER day_ini, day_end, annee_ref, day_ref
    21       REAL      dt, jD_ref, jH_ref
     23      REAL      dt, jD_ref, jH_ref, start_time
    2224      CHARACTER (len=10) :: calend
    2325

  • LMDZ5/branches/testing/libf/dyn3dpar/wrgrads.F

    r774 r1665  
    2222c   local
    2323
    24       integer im,jm,lm,i,j,l,lnblnk,iv,iii,iji,iif,ijf
     24      integer im,jm,lm,i,j,l,iv,iii,iji,iif,ijf
    2525
    2626      logical writectl
    … …  
    5555            nvar(if)=ivar(if)
    5656            var(ivar(if),if)=name
    57             tvar(ivar(if),if)=titlevar(1:lnblnk(titlevar))
     57            tvar(ivar(if),if)=trim(titlevar)
    5858            nld(ivar(if),if)=nl
    5959            print*,'initialisation ecriture de ',var(ivar(if),if)
    … …  
    9696      file=fichier(if)
    9797c   WARNING! on reecrase le fichier .ctl a chaque ecriture
    98       open(unit(if),file=file(1:lnblnk(file))//'.ctl'
     98      open(unit(if),file=trim(file)//'.ctl'
    9999     &         ,form='formatted',status='unknown')
    100100      write(unit(if),'(a5,1x,a40)')
    101      &       'DSET ','^'//file(1:lnblnk(file))//'.dat'
     101     &       'DSET ','^'//trim(file)//'.dat'
    102102
    103103      write(unit(if),'(a12)') 'UNDEF 1.0E30'

  • LMDZ5/branches/testing/libf/filtrez/coefils.h

    r1146 r1665  
    11!
    2 ! $Header$
     2! $Id $
    33!
    44      COMMON/coefils/jfiltnu,jfiltsu,jfiltnv,jfiltsv,sddu(iim),sddv(iim)&
    … …  
    77     & ,coefilu2(iim,jjm),coefilv2(iim,jjm)
    88!c
    9       INTEGER jfiltnu,jfiltsu,jfiltnv,jfiltsv,modfrstu,modfrstv
     9      INTEGER jfiltnu ! index of the last lat line filtered in NH (U grid)
     10      INTEGER jfiltsu ! index of the first lat line filtered in SH (U grid)
     11      INTEGER jfiltnv ! index of the last lat line filtered in NH (V grid)
     12      INTEGER jfiltsv ! index of the first lat line filtered in SH (V grid)
     13      INTEGER modfrstu ! number of retained (ie: unfiltered) modes on U grid
     14      INTEGER modfrstv ! number of retained (ie: unfiltered) modes on V grid
    1015      REAL    sddu,sddv,unsddu,unsddv,coefilu,coefilv,eignfnu,eignfnv
    1116      REAL    coefilu2,coefilv2

  • LMDZ5/branches/testing/libf/filtrez/filtreg_mod.F90

    r1279 r1665  
     1!
     2! $Id $
     3!
    14MODULE filtreg_mod
    25
    … …  
    4245    INTEGER ixmineq
    4346#endif
    44     EXTERNAL  inifgn
    4547    !
    4648    ! ------------------------------------------------------------
    … …  
    7173    CALL inifgn(eignvl)
    7274    !
    73     PRINT *,' EIGNVL '
     75    PRINT *,'inifilr: EIGNVL '
    7476    PRINT 250,eignvl
    75 250 FORMAT( 1x,5e13.6)
     77250 FORMAT( 1x,5e14.6)
    7678    !
    7779    ! compute eigenvalues and eigenfunctions
    … …  
    113115#endif
    114116    !
     117    ! For a regular grid, we want the filter to start at latitudes
     118    ! corresponding to lengths dx of the same size as dy (in terms
     119    ! of angles: dx=2*dy) => at colat0=0.5 (i.e. colatitude=30 degrees
     120    !  <=> latitude=60 degrees).
     121    ! Same idea for the zoomed grid: start filtering polewards as soon
     122    ! as length dx becomes of the same size as dy
    115123    !
    116124    colat0  =  MIN( 0.5, dymin/dxmin )
    … …  
    158166    imx  = iim
    159167    !
    160     PRINT *,' TRUNCATION AT ',imx
    161     !
     168    PRINT *,'inifilr: TRUNCATION AT ',imx
     169    !
     170! Ehouarn: set up some defaults
     171    jfiltnu=2 ! avoid north pole
     172    jfiltsu=jjm ! avoid south pole (which is at jjm+1)
     173    jfiltnv=1 ! NB: no poles on the V grid
     174    jfiltsv=jjm
     175
    162176    DO j = 2, jjm/2+1
    163177       cof = COS( rlatu(j) )/ colat0
    164178       IF ( cof .LT. 1. ) THEN
    165           IF( rlamda(imx) * COS(rlatu(j) ).LT.1. ) jfiltnu= j
     179          IF( rlamda(imx) * COS(rlatu(j) ).LT.1. ) THEN
     180            jfiltnu= j
     181          ENDIF
    166182       ENDIF
    167183
    168184       cof = COS( rlatu(jjp1-j+1) )/ colat0
    169185       IF ( cof .LT. 1. ) THEN
    170           IF( rlamda(imx) * COS(rlatu(jjp1-j+1) ).LT.1. ) &
     186          IF( rlamda(imx) * COS(rlatu(jjp1-j+1) ).LT.1. ) THEN
    171187               jfiltsu= jjp1-j+1
     188          ENDIF
    172189       ENDIF
    173190    ENDDO
    … …  
    176193       cof = COS( rlatv(j) )/ colat0
    177194       IF ( cof .LT. 1. ) THEN
    178           IF( rlamda(imx) * COS(rlatv(j) ).LT.1. ) jfiltnv= j
     195          IF( rlamda(imx) * COS(rlatv(j) ).LT.1. ) THEN
     196            jfiltnv= j
     197          ENDIF
    179198       ENDIF
    180199
    181200       cof = COS( rlatv(jjm-j+1) )/ colat0
    182201       IF ( cof .LT. 1. ) THEN
    183           IF( rlamda(imx) * COS(rlatv(jjm-j+1) ).LT.1. ) &
     202          IF( rlamda(imx) * COS(rlatv(jjm-j+1) ).LT.1. ) THEN
    184203               jfiltsv= jjm-j+1
     204          ENDIF
    185205       ENDIF
    186206    ENDDO
    187207    !                                 
    188208
    189     IF ( jfiltnu.LE.0 ) jfiltnu=1
    190209    IF( jfiltnu.GT. jjm/2 +1 )  THEN
    191210       PRINT *,' jfiltnu en dehors des valeurs acceptables ' ,jfiltnu
    … …  
    193212    ENDIF
    194213
    195     IF( jfiltsu.LE.0) jfiltsu=1
    196214    IF( jfiltsu.GT.  jjm  +1 )  THEN
    197215       PRINT *,' jfiltsu en dehors des valeurs acceptables ' ,jfiltsu
    … …  
    199217    ENDIF
    200218
    201     IF( jfiltnv.LE.0) jfiltnv=1
    202219    IF( jfiltnv.GT. jjm/2    )  THEN
    203220       PRINT *,' jfiltnv en dehors des valeurs acceptables ' ,jfiltnv
    … …  
    205222    ENDIF
    206223
    207     IF( jfiltsv.LE.0) jfiltsv=1
    208224    IF( jfiltsv.GT.     jjm  )  THEN
    209225       PRINT *,' jfiltsv en dehors des valeurs acceptables ' ,jfiltsv
    … …  
    211227    ENDIF
    212228
    213     PRINT *,' jfiltnv jfiltsv jfiltnu jfiltsu ' , &
     229    PRINT *,'inifilr: jfiltnv jfiltsv jfiltnu jfiltsu ' , &
    214230         jfiltnv,jfiltsv,jfiltnu,jfiltsu
    215231
    216232    IF(first_call_inifilr) THEN
    217233       ALLOCATE(matriceun(iim,iim,jfiltnu))
    218        ALLOCATE(matriceus(iim,iim,jfiltsu))
     234       ALLOCATE(matriceus(iim,iim,jjm-jfiltsu+1))
    219235       ALLOCATE(matricevn(iim,iim,jfiltnv))
    220        ALLOCATE(matricevs(iim,iim,jfiltsv))
     236       ALLOCATE(matricevs(iim,iim,jjm-jfiltsv+1))
    221237       ALLOCATE( matrinvn(iim,iim,jfiltnu))
    222        ALLOCATE( matrinvs(iim,iim,jfiltsu))
     238       ALLOCATE( matrinvs(iim,iim,jjm-jfiltsu+1))
    223239       first_call_inifilr = .FALSE.
    224240    ENDIF
    … …  
    230246    !
    231247    DO j = 1,jjm
     248    !default initialization: all modes are retained (i.e. no filtering)
    232249       modfrstu( j ) = iim
    233250       modfrstv( j ) = iim
    … …  
    306323
    307324    IF(jfiltnv.GE.jjm/2 .OR. jfiltnu.GE.jjm/2)THEN
    308 
     325! Ehouarn: and what are these for??? Trying to handle a limit case
     326!          where filters extend to and meet at the equator?
    309327       IF(jfiltnv.EQ.jfiltsv)jfiltsv=1+jfiltnv
    310328       IF(jfiltnu.EQ.jfiltsu)jfiltsu=1+jfiltnu
    … …  
    334352             eignft(i,k) = eignfnv(k,i) * coff
    335353          ENDDO
    336        ENDDO
     354       ENDDO ! of DO i=1,iim
    337355#ifdef CRAY
    338356       CALL MXM( eignfnv,iim,eignft,iim,matriceun(1,1,j),iim )
    … …  
    350368             ENDDO
    351369          ENDDO
    352        ENDDO
    353 #endif
    354 #endif
    355 
    356     ENDDO
     370       ENDDO ! of DO k = 1, iim
     371#endif
     372#endif
     373
     374    ENDDO ! of DO j = 2, jfiltnu
    357375
    358376    DO j = jfiltsu, jjm
    … …  
    364382             eignft(i,k) = eignfnv(k,i) * coff
    365383          ENDDO
    366        ENDDO
     384       ENDDO ! of DO i=1,iim
    367385#ifdef CRAY
    368386       CALL MXM(eignfnv,iim,eignft,iim,matriceus(1,1,j-jfiltsu+1),iim)
    … …  
    381399             ENDDO
    382400          ENDDO
    383        ENDDO
    384 #endif
    385 #endif
    386 
    387     ENDDO
     401       ENDDO ! of DO k = 1, iim
     402#endif
     403#endif
     404
     405    ENDDO ! of DO j = jfiltsu, jjm
    388406
    389407    !   ...................................................................
    … …  
    421439#endif
    422440
    423     ENDDO
     441    ENDDO ! of DO j = 1, jfiltnv
    424442
    425443    DO j = jfiltsv, jjm
    … …  
    452470#endif
    453471
    454     ENDDO
     472    ENDDO ! of DO j = jfiltsv, jjm
    455473
    456474    !   ...................................................................
    … …  
    488506#endif
    489507
    490     ENDDO
     508    ENDDO ! of DO j = 2, jfiltnu
    491509
    492510    DO j = jfiltsu, jjm
    … …  
    518536#endif
    519537
    520     ENDDO
     538    ENDDO ! of DO j = jfiltsu, jjm
    521539
    522540    IF (use_filtre_fft) THEN

  • LMDZ5/branches/testing/libf/grid/dimension/makdim

    r1146 r1665  
    1 for i in $* ; do
    2    list=$list.$i
     1#!/bin/bash
     2#set -xv
     3
     4# sanity check: do we have the required argument ?
     5if (( $# < 1 )) || (( $# > 3 ))
     6then
     7 echo "Wrong number of parameters in $0 !!!"
     8 echo " Usage:"
     9 echo "  $0 [im] [jm] lm"
     10 echo " where im, jm and lm are the dimensions"
     11 exit
     12fi
     13
     14# build "fichnom", the relevant 'dimensions.im.jm.lm' file name
     15for i in $*
     16  do
     17  list=$list.$i
    318done
    419fichdim=dimensions${list}
    520
    6 if [ ! -f $fichdim ] ; then
    7 # si le fichier de dimensions n'existe pas, on le cree
     21if [ ! -f $fichdim ]
     22    then
     23#    echo "$fichdim does not exist"
    824
    9   if [ $# -ge 3 ] ; then
    10      im=$1
    11      jm=$2
    12      lm=$3
    13      n2=$1
    14      ndm=1
     25    # assign values of im, jm and lm
     26    if [ $# -ge 3 ]
     27        then
     28        im=$1
     29        jm=$2
     30        lm=$3
     31        ndm=1
     32    elif [ $# -ge 2 ]
     33        then
     34        im=1
     35        jm=$1
     36        lm=$2
     37        ndm=1
     38    elif [ $# -ge 1 ]
     39        then
     40        im=1
     41        jm=1
     42        lm=$1
     43        ndm=1
     44    fi
    1545
    16 # Le test suivant est commente car il est inutile avec le nouveau
    17 # filtre filtrez. Attention avec le "vieux" filtre (F. Forget,11/1994)
    18 #
    19 #     while [ "$n2" -gt 2 ]; do
    20 #       n2=`expr $n2 / 2`
    21 #       ndm=`expr $ndm + 1`
    22 #       echo $n2
    23 #    done
    24 #    if [ "$n2" != 2 ] ; then
    25 #       echo le nombre de longitude doit etre une puissance de 2
    26 #       exit
    27 #    fi
    28 
    29 
    30   else if [ $# -ge 2 ] ; then
    31       im=1
    32       jm=$1
    33       lm=$2
    34       ndm=1
    35   else if [ $# -ge 1 ] ; then
    36       im=1
    37       jm=1
    38       lm=$1
    39       ndm=1
    40   else
    41          echo il faut au moins une dimension
    42          exit
    43   fi
    44 fi
    45 fi
    46 
     46# since the file doesn't exist, we create it
    4747cat << EOF > $fichdim
    4848!-----------------------------------------------------------------------
    … …  
    6262fi
    6363
     64# remove 'old' dimensions.h file (if any) and replace it with new one
     65if [ -f ../dimensions.h ] ; then
    6466\rm ../dimensions.h
     67fi
    6568tar cf - $fichdim | ( cd .. ; tar xf - ; mv $fichdim dimensions.h )
     69# line above is a trick to preserve time of creation of dimensions.h files

  • LMDZ5/branches/testing/libf/phylmd/calwake.F

    r1403 r1665  
    6060      REAL wake_dtKE(klon,klev),wake_dqKE(klon,klev)
    6161      REAL wake_dtPBL(klon,klev),wake_dqPBL(klon,klev)
    62       REAL wake_omg(klon,klev+1),wake_dp_deltomg(klon,klev)
     62      REAL wake_omg(klon,klev),wake_dp_deltomg(klon,klev)
    6363      REAL wake_spread(klon,klev),wake_Cstar(klon)
    6464      REAL wake_ddeltat(klon,klev),wake_ddeltaq(klon,klev)
    … …  
    8484      REAL tu(klon,klev),qu(klon,klev)
    8585      REAL hw(klon),sigmaw(klon),wape(klon),fip(klon),gfl(klon)
    86       REAL omgbdth(klon,klev),dp_omgb(klon,klev)
     86      REAL omgbdth(klon,klev+1),dp_omgb(klon,klev)
    8787      REAL dtKE(klon,klev),dqKE(klon,klev)
    8888      REAL dtPBL(klon,klev),dqPBL(klon,klev)

  • LMDZ5/branches/testing/libf/phylmd/climb_hq_mod.F90

    r1084 r1665  
    252252       Bcoef(i) = -1. * RG / buf
    253253    END DO
     254    acoef(knon+1: klon) = 0.
     255    bcoef(knon+1: klon) = 0.
    254256
    255257  END SUBROUTINE calc_coef

  • LMDZ5/branches/testing/libf/phylmd/climb_wind_mod.F90

    r1067 r1665  
    209209       Bcoef(i) = -RG/buf
    210210    END DO
     211    acoef(knon+1: klon) = 0.
     212    bcoef(knon+1: klon) = 0.
    211213
    212214  END SUBROUTINE calc_coef

  • LMDZ5/branches/testing/libf/phylmd/coef_diff_turb_mod.F90

    r1067 r1665  
    389389!           calculer la fraction nuageuse (processus humide):
    390390!
    391           zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)
     391          if (zq /= 0.) then
     392             zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)
     393          else
     394             zfr = 0.
     395          end if
    392396          zfr = MAX(0.0,MIN(1.0,zfr))
    393397          IF (.NOT.richum) zfr = 0.0

  • LMDZ5/branches/testing/libf/phylmd/concvl.F

    r1664 r1665  
    248248         DO i = 1, klon
    249249          cbmf(i) = 0.
    250           plcl(i) = 0.
     250!          plcl(i) = 0.
    251251          sigd(i) = 0.
    252252         ENDDO
    … …  
    256256      plfc(:)  = 0.
    257257      wbeff(:) = 100.
     258      plcl(:) = 0.
    258259
    259260      DO k = 1, klev+1

  • LMDZ5/branches/testing/libf/phylmd/cpl_mod.F90

    r1454 r1665  
    345345       IF (is_sequential) THEN
    346346          ndexcs(:) = 0
    347           itau_w = itau_phy + itime
     347          itau_w = itau_phy + itime + start_time * day_step / iphysiq
    348348          DO i = 1, maxrecv
    349349            IF (inforecv(i)%action) THEN
    … …  
    12321232    IF (is_sequential) THEN
    12331233       ndexct(:) = 0
    1234        itau_w = itau_phy + itime
     1234       itau_w = itau_phy + itime + start_time * day_step / iphysiq
    12351235       CALL histwrite(nidct,'tauxe',itau_w,tmp_taux,iim*(jjm+1),ndexct)
    12361236       CALL histwrite(nidct,'tauyn',itau_w,tmp_tauy,iim*(jjm+1),ndexct)

  • LMDZ5/branches/testing/libf/phylmd/iostart.F90

    r1403 r1665  
    177177         ierr=NF90_GET_VAR(nid_start,varid,field_glo)
    178178         IF (ierr/=NF90_NOERR) THEN
     179           ! La variable exist dans le fichier mais la lecture a echouee.
    179180           PRINT*, 'phyetat0: Lecture echouee pour <'//field_name//'>'
    180            CALL abort
     181
     182           IF (field_name=='CLWCON' .OR. field_name=='RNEBCON' .OR. field_name=='RATQS') THEN
     183              ! Essaye de lire le variable sur surface uniqument, comme fait avant
     184              field_glo(:)=0.
     185              ierr=NF90_GET_VAR(nid_start,varid,field_glo(1:klon_glo))
     186              IF (ierr/=NF90_NOERR) THEN
     187                 PRINT*, 'phyetat0: Lecture echouee aussi en 2D pour <'//field_name//'>'
     188                 CALL abort
     189              ELSE
     190                 PRINT*, 'phyetat0: La variable <'//field_name//'> lu sur surface seulement'!, selon ancien format, le reste mis a zero'
     191              END IF
     192           ELSE
     193              CALL abort
     194           ENDIF
    181195         ENDIF
    182196

  • LMDZ5/branches/testing/libf/phylmd/limit_read_mod.F90

    r1001 r1665  
    150150   
    151151    INCLUDE "indicesol.h"
     152    INCLUDE "iniprint.h"
    152153
    153154! In- and ouput arguments
    … …  
    165166!$OMP THREADPRIVATE(lmt_pas)
    166167    LOGICAL, SAVE                             :: first_call=.TRUE.
    167 !$OMP THREADPRIVATE(first_call)   
     168!$OMP THREADPRIVATE(first_call) 
     169    INTEGER, SAVE                             :: jour_lu = -1
     170!$OMP THREADPRIVATE(jour_lu) 
    168171! Locals variables
    169172!****************************************************************************************
    … …  
    209212
    210213    is_modified = .FALSE.
    211     IF (MOD(itime-1, lmt_pas) == 0) THEN   ! time to read
     214    IF (MOD(itime-1, lmt_pas) == 0 .OR. jour_lu /= jour ) THEN   ! time to read
     215       jour_lu = jour
    212216       is_modified = .TRUE.
    213217!$OMP MASTER  ! Only master thread

  • LMDZ5/branches/testing/libf/phylmd/phyetat0.F

    r1458 r1665  
    749749      ENDIF
    750750
    751       u_ancien = 0.0   !AXC: We don't have u_ancien and v_ancien in the start
    752       v_ancien = 0.0   !AXC: files, therefore they have to be initialized.
    753 c
     751      CALL get_field("UANCIEN",u_ancien,found)
     752      IF (.NOT. found) THEN
     753         PRINT*, "phyetat0: Le champ <UANCIEN> est absent"
     754         PRINT*, "Depart legerement fausse. Mais je continue"
     755         ancien_ok = .FALSE.
     756      ENDIF
     757
     758      CALL get_field("VANCIEN",v_ancien,found)
     759      IF (.NOT. found) THEN
     760         PRINT*, "phyetat0: Le champ <VANCIEN> est absent"
     761         PRINT*, "Depart legerement fausse. Mais je continue"
     762         ancien_ok = .FALSE.
     763      ENDIF
    754764
    755765      clwcon=0.
    756       CALL get_field("CLWCON",clwcon(:,1),found)
     766      CALL get_field("CLWCON",clwcon,found)
    757767      IF (.NOT. found) THEN
    758768         PRINT*, "phyetat0: Le champ CLWCON est absent"
    … …  
    766776c
    767777      rnebcon = 0.
    768       CALL get_field("RNEBCON",rnebcon(:,1),found)
     778      CALL get_field("RNEBCON",rnebcon,found)
    769779      IF (.NOT. found) THEN
    770780         PRINT*, "phyetat0: Le champ RNEBCON est absent"
    … …  
    781791c
    782792      ratqs=0.
    783       CALL get_field("RATQS",ratqs(:,1),found)
     793      CALL get_field("RATQS",ratqs,found)
    784794      IF (.NOT. found) THEN
    785795         PRINT*, "phyetat0: Le champ <RATQS> est absent"

  • LMDZ5/branches/testing/libf/phylmd/phyredem.F

    r1458 r1665  
    267267     
    268268      CALL put_field("QANCIEN","QANCIEN",q_ancien)
    269      
     269
     270      CALL put_field("UANCIEN","",u_ancien)
     271
     272      CALL put_field("VANCIEN","",v_ancien)
     273
    270274      CALL put_field("RUGMER","Longueur de rugosite sur mer",
    271275     .               frugs(:,is_oce))
    272276     
    273       CALL put_field("CLWCON","Eau liquide convective",clwcon(:,1))
    274      
    275       CALL put_field("RNEBCON","Nebulosite convective",rnebcon(:,1))
    276      
    277       CALL put_field("RATQS", "Ratqs",ratqs(:,1))
     277      CALL put_field("CLWCON","Eau liquide convective",clwcon)
     278     
     279      CALL put_field("RNEBCON","Nebulosite convective",rnebcon)
     280     
     281      CALL put_field("RATQS", "Ratqs",ratqs)
    278282c
    279283c run_off_lic_0

  • LMDZ5/branches/testing/libf/phylmd/phys_output_mod.F90

    r1664 r1665  
    13181318          ! Couplage conv-CL
    13191319          IF (iflag_con.GE.3) THEN
    1320              IF (iflag_coupl>=1) THEN
    13211320                CALL histdef2d(iff,clef_stations(iff), &
    13221321                     o_ale_bl%flag,o_ale_bl%name, "ALE BL", "m2/s2")
    13231322                CALL histdef2d(iff,clef_stations(iff), &
    13241323                     o_alp_bl%flag,o_alp_bl%name, "ALP BL", "m2/s2")
    1325              ENDIF
    13261324          ENDIF !(iflag_con.GE.3)
    13271325
    … …  
    14911489                CALL histdef2d(iff,clef_stations(iff), &
    14921490                     o_alp_wk%flag,o_alp_wk%name, "ALP WK", "m2/s2")
    1493                 CALL histdef2d(iff,clef_stations(iff), &
    1494                      o_ale%flag,o_ale%name, "ALE", "m2/s2")
    1495                 CALL histdef2d(iff,clef_stations(iff), &
    1496                      o_alp%flag,o_alp%name, "ALP", "W/m2")
    1497                 CALL histdef2d(iff,clef_stations(iff),o_cin%flag,o_cin%name, "Convective INhibition", "m2/s2")
    1498                 CALL histdef2d(iff,clef_stations(iff),o_wape%flag,o_WAPE%name, "WAPE", "m2/s2")
    14991491                CALL histdef2d(iff,clef_stations(iff),o_wake_h%flag,o_wake_h%name, "wake_h", "-")
    15001492                CALL histdef2d(iff,clef_stations(iff),o_wake_s%flag,o_wake_s%name, "wake_s", "-")
    … …  
    15041496                CALL histdef3d(iff,clef_stations(iff),o_wake_deltaq%flag,o_wake_deltaq%name, "wake_deltaq", " ")
    15051497                CALL histdef3d(iff,clef_stations(iff),o_wake_omg%flag,o_wake_omg%name, "wake_omg", " ")
     1498                CALL histdef2d(iff,clef_stations(iff),o_wape%flag,o_WAPE%name, "WAPE", "m2/s2")
    15061499             ENDIF
     1500             CALL histdef2d(iff,clef_stations(iff), &
     1501                     o_ale%flag,o_ale%name, "ALE", "m2/s2")
     1502             CALL histdef2d(iff,clef_stations(iff), &
     1503                     o_alp%flag,o_alp%name, "ALP", "W/m2")
     1504             CALL histdef2d(iff,clef_stations(iff),o_cin%flag,o_cin%name, "Convective INhibition", "m2/s2")
    15071505             CALL histdef3d(iff,clef_stations(iff),o_Vprecip%flag,o_Vprecip%name, "precipitation vertical profile", "-")
    15081506             CALL histdef3d(iff,clef_stations(iff),o_ftd%flag,o_ftd%name, "tend temp due aux descentes precip", "-")
    … …  
    18591857    if ( type == 'day'.or.type == 'days'.or.type == 'jours'.or.type == 'jour' ) timestep = ttt * dayseconde
    18601858    if ( type == 'mounths'.or.type == 'mth'.or.type == 'mois' ) then
    1861        write(lunout,*)'annee_ref,day_ref mon_len',annee_ref,day_ref,ioget_mon_len(annee_ref,day_ref)
     1859       write(lunout,*)'annee_ref,day_ref mon_len',annee_ref,day_ref,mth_len
    18621860       timestep = ttt * dayseconde * mth_len
    18631861    endif

  • LMDZ5/branches/testing/libf/phylmd/phys_output_write.h

    r1664 r1665  
    1       itau_w = itau_phy + itap
     1      itau_w = itau_phy + itap + start_time * day_step / iphysiq
    22
    33      DO iff=1,nfiles
    … …  
    801801! Couplage convection-couche limite
    802802      IF (iflag_con.GE.3) THEN
    803       IF (iflag_coupl>=1) THEN
    804803       IF (o_ale_bl%flag(iff)<=lev_files(iff)) THEN
    805804       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
    … …  
    810809     $o_alp_bl%name,itau_w,alp_bl)
    811810       ENDIF
    812       ENDIF !iflag_coupl>=1
    813811      ENDIF !(iflag_con.GE.3)
    814812
    … …  
    825823       ENDIF
    826824
    827        IF (o_ale%flag(iff)<=lev_files(iff)) THEN
    828        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
    829      $o_ale%name,itau_w,ale)
    830        ENDIF
    831        IF (o_alp%flag(iff)<=lev_files(iff)) THEN
    832        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
    833      $o_alp%name,itau_w,alp)
    834        ENDIF
    835        IF (o_cin%flag(iff)<=lev_files(iff)) THEN
    836        CALL histwrite_phy(nid_files(iff),clef_stations(iff),
    837      $o_cin%name,itau_w,cin)
    838        ENDIF
    839825       IF (o_wape%flag(iff)<=lev_files(iff)) THEN
    840826       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
    … …  
    883869      ENDIF ! iflag_wake>=1
    884870
     871       IF (o_ale%flag(iff)<=lev_files(iff)) THEN
     872       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     873     $o_ale%name,itau_w,ale)
     874       ENDIF
     875       IF (o_alp%flag(iff)<=lev_files(iff)) THEN
     876       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     877     $o_alp%name,itau_w,alp)
     878       ENDIF
     879       IF (o_cin%flag(iff)<=lev_files(iff)) THEN
     880       CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     881     $o_cin%name,itau_w,cin)
     882       ENDIF
    885883        IF (o_Vprecip%flag(iff)<=lev_files(iff)) THEN
    886884       CALL histwrite_phy(nid_files(iff),clef_stations(iff),

  • LMDZ5/branches/testing/libf/phylmd/physiq.F

    r1664 r1665  
    614614      REAL dd_t(klon,klev),dd_q(klon,klev)
    615615
    616       real, save :: alp_bl_prescr=0.
    617       real, save :: ale_bl_prescr=0.
     616      real, save :: alp_bl_prescr=0.1
     617      real, save :: ale_bl_prescr=4.
    618618
    619619      real, save :: ale_max=1000.
    … …  
    791791cIM
    792792      EXTERNAL haut2bas  !variables de haut en bas
    793       INTEGER lnblnk1
    794       EXTERNAL lnblnk1   !enleve les blancs a la fin d'une variable de type
    795                          !caracter
    796793      EXTERNAL ini_undefSTD  !initialise a 0 une variable a 1 niveau de pression
    797794      EXTERNAL undefSTD      !somme les valeurs definies d'1 var a 1 niveau de pression
    … …  
    12201217      INTEGER :: nbtr_tmp ! Number of tracer inside concvl
    12211218      REAL, dimension(klon,klev) :: sh_in ! Specific humidity entering in phytrac
     1219      integer iostat
    12221220
    12231221cIM for NMC files
    … …  
    14921490cCR:04.12.07: initialisations poches froides
    14931491c Controle de ALE et ALP pour la fermeture convective (jyg)
    1494           if (iflag_wake>=1) then
    1495             CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
     1492         CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
    14961493     s                  ,alp_bl_prescr, ale_bl_prescr)
    14971494c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
    14981495c        print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
    1499           endif
    15001496
    15011497        do i = 1,klon
    … …  
    15081504      nCFMIP=npCFMIP
    15091505      OPEN(98,file='npCFMIP_param.data',status='old',
    1510      $          form='formatted',err=999)
     1506     $          form='formatted',iostat=iostat)
     1507            if (iostat == 0) then
    15111508      READ(98,*,end=998) nCFMIP
    15121509998   CONTINUE
    … …  
    15401537     $tabijGCM, lonGCM, latGCM, iGCM, jGCM)
    15411538c
    1542 999      CONTINUE
     1539            else
     1540               ALLOCATE(tabijGCM(0))
     1541               ALLOCATE(lonGCM(0), latGCM(0))
     1542               ALLOCATE(iGCM(0), jGCM(0))
     1543            end if
    15431544         ENDIF !debut
    15441545 
    … …  
    20412042c
    20422043
    2043       CALL pbl_surface(
    2044      e     dtime,     date0,     itap,    days_elapsed+1,
    2045      e     debut,     lafin,
    2046      e     rlon,      rlat,      rugoro,  rmu0,     
    2047      e     rain_fall, snow_fall, solsw,   sollw,   
    2048      e     t_seri,    q_seri,    u_seri,  v_seri,   
    2049      e     pplay,     paprs,     pctsrf,           
    2050      +     ftsol,     falb1,     falb2,   u10m,   v10m,
    2051      s     sollwdown, cdragh,    cdragm,  u1,    v1,
    2052      s     albsol1,   albsol2,   sens,    evap, 
    2053      s     zxtsol,    zxfluxlat, zt2m,    qsat2m,
    2054      s     d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf,
    2055      s     coefh,     coefm,     slab_wfbils,               
    2056      d     qsol,      zq2m,      s_pblh,  s_lcl,
    2057      d     s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
    2058      d     s_therm,   s_trmb1,   s_trmb2, s_trmb3,
    2059      d     zxrugs,    zu10m,     zv10m,   fder,
    2060      d     zxqsurf,   rh2m,      zxfluxu, zxfluxv,
    2061      d     frugs,     agesno,    fsollw,  fsolsw,
    2062      d     d_ts,      fevap,     fluxlat, t2m,
    2063      d     wfbils,    wfbilo,    fluxt,   fluxu,  fluxv,
    2064      -     dsens,     devap,     zxsnow,
    2065      -     zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke )
     2044      if (iflag_pbl/=0) then
     2045
     2046        CALL pbl_surface(
     2047     e       dtime,     date0,     itap,    days_elapsed+1,
     2048     e       debut,     lafin,
     2049     e       rlon,      rlat,      rugoro,  rmu0,     
     2050     e       rain_fall, snow_fall, solsw,   sollw,   
     2051     e       t_seri,    q_seri,    u_seri,  v_seri,   
     2052     e       pplay,     paprs,     pctsrf,           
     2053     +       ftsol,     falb1,     falb2,   u10m,   v10m,
     2054     s       sollwdown, cdragh,    cdragm,  u1,    v1,
     2055     s       albsol1,   albsol2,   sens,    evap, 
     2056     s       zxtsol,    zxfluxlat, zt2m,    qsat2m,
     2057     s       d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf,
     2058     s       coefh,     coefm,     slab_wfbils,               
     2059     d       qsol,      zq2m,      s_pblh,  s_lcl,
     2060     d       s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
     2061     d       s_therm,   s_trmb1,   s_trmb2, s_trmb3,
     2062     d       zxrugs,    zu10m,     zv10m,   fder,
     2063     d       zxqsurf,   rh2m,      zxfluxu, zxfluxv,
     2064     d       frugs,     agesno,    fsollw,  fsolsw,
     2065     d       d_ts,      fevap,     fluxlat, t2m,
     2066     d       wfbils,    wfbilo,    fluxt,   fluxu,  fluxv,
     2067     -       dsens,     devap,     zxsnow,
     2068     -       zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke )
    20662069
    20672070
    20682071!-----------------------------------------------------------------------------------------
    20692072! ajout des tendances de la diffusion turbulente
    2070       CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
     2073        CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
    20712074!-----------------------------------------------------------------------------------------
    20722075
    2073       if (mydebug) then
    2074         call writefield_phy('u_seri',u_seri,llm)
    2075         call writefield_phy('v_seri',v_seri,llm)
    2076         call writefield_phy('t_seri',t_seri,llm)
    2077         call writefield_phy('q_seri',q_seri,llm)
    2078       endif
    2079 
    2080 
    2081       IF (ip_ebil_phy.ge.2) THEN
    2082         ztit='after surface_main'
    2083         CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
     2076        if (mydebug) then
     2077          call writefield_phy('u_seri',u_seri,llm)
     2078          call writefield_phy('v_seri',v_seri,llm)
     2079          call writefield_phy('t_seri',t_seri,llm)
     2080          call writefield_phy('q_seri',q_seri,llm)
     2081        endif
     2082
     2083
     2084        IF (ip_ebil_phy.ge.2) THEN
     2085          ztit='after surface_main'
     2086          CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
    20842087     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
    20852088     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
    2086          call diagphy(airephy,ztit,ip_ebil_phy
     2089          call diagphy(airephy,ztit,ip_ebil_phy
    20872090     e      , zero_v, zero_v, zero_v, zero_v, sens
    20882091     e      , evap  , zero_v, zero_v, ztsol
    20892092     e      , d_h_vcol, d_qt, d_ec
    20902093     s      , fs_bound, fq_bound )
    2091       END IF
     2094        END IF
     2095
     2096      ENDIF
    20922097
    20932098c =================================================================== c
    … …  
    22392244cdans le thermique sinon
    22402245       if (iflag_coupl.eq.0) then
    2241           if (debut.and.prt_level.gt.9)
    2242      $                     WRITE(lunout,*)'ALE et ALP imposes'
    2243           do i = 1,klon
    2244 con ne couple que ale
    2245 c           ALE(i) = max(ale_wake(i),Ale_bl(i))
    2246             ALE(i) = max(ale_wake(i),ale_bl_prescr)
    2247 con ne couple que alp
    2248 c           ALP(i) = alp_wake(i) + Alp_bl(i)
    2249             ALP(i) = alp_wake(i) + alp_bl_prescr
    2250           enddo
     2246          if (debut.and.prt_level.gt.9)WRITE(lunout,*) 'ALE&ALP imposes'
     2247          Ale_bl(1:klon) = ale_bl_prescr
     2248          Alp_bl(1:klon) = alp_bl_prescr
    22512249       else
    22522250         IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
    2253 !         do i = 1,klon
    2254 !             ALE(i) = max(ale_wake(i),Ale_bl(i))
    2255 ! avant        ALP(i) = alp_wake(i) + Alp_bl(i)
    2256 !             ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
    2257 !         write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
    2258 !         write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
    2259 !         enddo
     2251       endif
    22602252
    22612253!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    22642256! w si <0
    22652257!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     2258
    22662259       do i = 1,klon
    22672260          ALE(i) = max(ale_wake(i),Ale_bl(i))
    … …  
    22762269          endif
    22772270       enddo
     2271
    22782272!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    22792273
    2280        endif
    22812274       do i=1,klon
    22822275          if (alp(i)>alp_max) then
    … …  
    26342627c  ==============
    26352628
    2636 ! Dans le cas où on active les thermiques, on fait partir l'ajustement
     2629! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
    26372630! a partir du sommet des thermiques.
    26382631! Dans le cas contraire, on demarre au niveau 1.
    … …  
    28212814! FH 22/09/2009
    28222815! La ligne ci-dessous faisait osciller le modele et donnait une solution
    2823 ! assymptotique bidon et dépendant fortement du pas de temps.
     2816! asymptotique bidon et d\'ependant fortement du pas de temps.
    28242817!        ratqs(:,:)=sqrt(ratqs(:,:)**2+ratqss(:,:)**2)
    28252818!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    35763569
    35773570       IF (ok_hines) then
    3578 
    35793571         CALL hines_gwd(klon,klev,dtime,paprs,pplay,
    35803572     i                  rlat,t_seri,u_seri,v_seri,
    … …  
    35843576c  ajout des tendances
    35853577        CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'hin')
    3586 
    35873578      ENDIF
    3588 c
    3589 
    3590 c
    3591 cIM cf. FLott BEG
     3579
    35923580C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
    35933581
    … …  
    36143602cIM calcul composantes axiales du moment angulaire et couple des montagnes
    36153603c
    3616       IF (is_sequential) THEN
     3604      IF (is_sequential .and. ok_orodr) THEN
    36173605     
    36183606        CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur,

  • LMDZ5/branches/testing/libf/phylmd/phytrac.F90

    r1664 r1665  
    214214     SELECT CASE(type_trac)
    215215     CASE('lmdz')
    216         CALL traclmdz_init(pctsrf, ftsol, tr_seri, t_seri, pplay, sh, pdtphys, aerosol, lessivage)
     216        CALL traclmdz_init(pctsrf, xlat, xlon, ftsol, tr_seri, t_seri, pplay, sh, pdtphys, aerosol, lessivage)
    217217     CASE('inca')
    218218        source(:,:)=0.

  • LMDZ5/branches/testing/libf/phylmd/readaerosol.F90

    r1492 r1665  
    247247 
    248248       WRITE(lunout,*) 'reading variable ',TRIM(varname),' in file ', TRIM(fname)
    249        CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, ncid) )
     249       CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, ncid), "pb open "//trim(varname) )
    250250
    251251! Test for equal longitudes and latitudes in file and model
    252252!****************************************************************************************
    253253       ! Read and test longitudes
    254        CALL check_err( nf90_inq_varid(ncid, 'lon', varid) )
    255        CALL check_err( nf90_get_var(ncid, varid, lon_src(:)) )
     254       CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" )
     255       CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" )
    256256       
    257257       IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN
    … …  
    264264
    265265       ! Read and test latitudes
    266        CALL check_err( nf90_inq_varid(ncid, 'lat', varid) )
    267        CALL check_err( nf90_get_var(ncid, varid, lat_src(:)) )
     266       CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" )
     267       CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" )
    268268
    269269       ! Invert source latitudes
    … …  
    311311! 2) Find vertical dimension klev_src
    312312!****************************************************************************************
    313        CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src) )
     313       CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src),"pb inq dim for PRESNIVS or lev" )
    314314       
    315315     ! Allocate variables depending on the number of vertical levels
    … …  
    330330!**************************************************************************************************
    331331       ierr = nf90_inq_dimid(ncid, 'TIME',dimid)
    332        CALL check_err( nf90_inquire_dimension(ncid, dimid, len = nbr_tsteps) )
     332       CALL check_err( nf90_inquire_dimension(ncid, dimid, len = nbr_tsteps),"pb inq dim TIME" )
    333333!       IF (nbr_tsteps /= 12 .AND. nbr_tsteps /= 14) THEN
    334334       IF (nbr_tsteps /= 12 ) THEN
    … …  
    339339!****************************************************************************************
    340340          ! Get variable id
    341           CALL check_err( nf90_inq_varid(ncid, TRIM(varname), varid) )
     341          CALL check_err( nf90_inq_varid(ncid, TRIM(varname), varid),"pb inq var "//TRIM(varname) )
    342342         
    343343          ! Get the variable
    344           CALL check_err( nf90_get_var(ncid, varid, varyear(:,:,:,:)) )
     344          CALL check_err( nf90_get_var(ncid, varid, varyear(:,:,:,:)),"pb get var "//TRIM(varname) )
    345345         
    346346! ++) Read surface pression, 12 month in one variable
    347347!****************************************************************************************
    348348          ! Get variable id
    349           CALL check_err( nf90_inq_varid(ncid, "ps", varid) )
     349          CALL check_err( nf90_inq_varid(ncid, "ps", varid),"pb inq var ps" )
    350350          ! Get the variable
    351           CALL check_err( nf90_get_var(ncid, varid, psurf_glo2D) )
     351          CALL check_err( nf90_get_var(ncid, varid, psurf_glo2D),"pb get var ps" )
    352352         
    353353! ++) Read mass load, 12 month in one variable
    354354!****************************************************************************************
    355355          ! Get variable id
    356           CALL check_err( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) )
     356          CALL check_err( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) ,"pb inq var load_"//TRIM(varname))
    357357          ! Get the variable
    358           CALL check_err( nf90_get_var(ncid, varid, load_glo2D) )
     358          CALL check_err( nf90_get_var(ncid, varid, load_glo2D),"pb get var load_"//TRIM(varname) )
    359359         
    360360! ++) Read ap
    361361!****************************************************************************************
    362362          ! Get variable id
    363           CALL check_err( nf90_inq_varid(ncid, "ap", varid) )
     363          CALL check_err( nf90_inq_varid(ncid, "ap", varid),"pb inq var ap" )
    364364          ! Get the variable
    365           CALL check_err( nf90_get_var(ncid, varid, pt_ap) )
     365          CALL check_err( nf90_get_var(ncid, varid, pt_ap),"pb get var ap" )
    366366
    367367! ++) Read b
    368368!****************************************************************************************
    369369          ! Get variable id
    370           CALL check_err( nf90_inq_varid(ncid, "b", varid) )
     370          CALL check_err( nf90_inq_varid(ncid, "b", varid),"pb inq var b" )
    371371          ! Get the variable
    372           CALL check_err( nf90_get_var(ncid, varid, pt_b) )
     372          CALL check_err( nf90_get_var(ncid, varid, pt_b),"pb get var b" )
    373373
    374374! ++) Read p0 : reference pressure
    375375!****************************************************************************************
    376376          ! Get variable id
    377           CALL check_err( nf90_inq_varid(ncid, "p0", varid) )
     377          CALL check_err( nf90_inq_varid(ncid, "p0", varid),"pb inq var p0" )
    378378          ! Get the variable
    379           CALL check_err( nf90_get_var(ncid, varid, p0) )
     379          CALL check_err( nf90_get_var(ncid, varid, p0),"pb get var p0" )
    380380         
    381381
    … …  
    412412             
    413413             ! Get variable id
    414              CALL check_err( nf90_inq_varid(ncid, TRIM(cvar), varid) )
     414             CALL check_err( nf90_inq_varid(ncid, TRIM(cvar), varid),"pb inq var "//TRIM(cvar) )
    415415             
    416416             ! Get the variable
    417              CALL check_err( nf90_get_var(ncid, varid, varmth) )
     417             CALL check_err( nf90_get_var(ncid, varid, varmth),"pb get var "//TRIM(cvar) )
    418418             
    419419             ! Store in variable for the whole year
    … …  
    432432! 4) Close file 
    433433!****************************************************************************************
    434        CALL check_err( nf90_close(ncid) )
     434       CALL check_err( nf90_close(ncid),"pb in close" )
    435435     
    436436
    … …  
    570570
    571571
    572   SUBROUTINE check_err(status)
     572  SUBROUTINE check_err(status,text)
    573573    USE netcdf
    574574    IMPLICIT NONE
    … …  
    576576    INCLUDE "iniprint.h"
    577577    INTEGER, INTENT (IN) :: status
     578    CHARACTER(len=*), INTENT (IN), OPTIONAL :: text
    578579
    579580    IF (status /= NF90_NOERR) THEN
    580        WRITE(lunout,*) 'Error in get_aero_fromfile ',status
     581       WRITE(lunout,*) 'Error in get_aero_fromfile, netcdf error code = ',status
     582       IF (PRESENT(text)) THEN
     583          WRITE(lunout,*) 'Error in get_aero_fromfile : ',text
     584       END IF
    581585       CALL abort_gcm('get_aero_fromfile',trim(nf90_strerror(status)),1)
    582586    END IF

  • LMDZ5/branches/testing/libf/phylmd/readaerosol_interp.F90

    r1492 r1665  
    175175     ! Reading values corresponding to the closest year taking into count the choice of aer_type.
    176176     ! For aer_type=scenario interpolation between 2 data sets is done in readaerosol.
    177      ! If aer_type=mix1 or mix2, the run type and file name depends on the aerosol.
     177     ! If aer_type=mix1, mix2 or mix3, the run type and file name depends on the aerosol.
    178178     IF (aer_type=='preind' .OR. aer_type=='actuel' .OR. aer_type=='annuel' .OR. aer_type=='scenario') THEN
    179179        ! Standard case
    … …  
    196196        ELSE
    197197           filename='aerosols'
     198           type='preind'
     199        END IF
     200     ELSE  IF (aer_type == 'mix3') THEN
     201        ! Special case using a mix of annual sulfate file and natrual aerosols
     202        IF (name_aero(id_aero) == 'SO4') THEN
     203           filename='aerosols'
     204           type='annuel'
     205        ELSE
     206           filename='aerosols'
    198207           type='preind'
    199208        END IF

  • LMDZ5/branches/testing/libf/phylmd/traclmdz_mod.F90

    r1454 r1665  
    8484
    8585
    86   SUBROUTINE traclmdz_init(pctsrf, ftsol, tr_seri, t_seri, pplay, sh, pdtphys, aerosol, lessivage)
     86  SUBROUTINE traclmdz_init(pctsrf, xlat, xlon, ftsol, tr_seri, t_seri, pplay, sh, pdtphys, aerosol, lessivage)
    8787    ! This subroutine allocates and initialize module variables and control variables.
    8888    ! Initialization of the tracers should be done here only for those not found in the restart file.
    … …  
    9999! Input variables
    100100    REAL,DIMENSION(klon,nbsrf),INTENT(IN)     :: pctsrf ! Pourcentage de sol f(nature du sol)
     101    REAL,DIMENSION(klon),INTENT(IN)           :: xlat   ! latitudes en degres pour chaque point
     102    REAL,DIMENSION(klon),INTENT(IN)           :: xlon   ! longitudes en degres pour chaque point
    101103    REAL,DIMENSION(klon,nbsrf),INTENT(IN)     :: ftsol  ! Temperature du sol (surf)(Kelvin)
    102104    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri! Concentration Traceur [U/KgA] 

  • LMDZ5/branches/testing/libf/phylmd/wrgradsfi.F

    r776 r1665  
    2323c   local
    2424
    25       integer lm,l,lnblnk
     25      integer lm,l
    2626
    2727

  • LMDZ5/branches/testing/libf/phylmd/write_bilKP_ave.h

    r776 r1665  
    99c Champs 2D:
    1010c
    11       itau_w = itau_phy + itap
     11      itau_w = itau_phy + itap + start_time * day_step / iphysiq
    1212c
    1313cym      CALL gr_fi_ecrit(klev, klon,iim,jjmp1, ue_lay,zx_tmp_3d)

  • LMDZ5/branches/testing/libf/phylmd/write_bilKP_ins.h

    r776 r1665  
    77      ndex3d = 0
    88c
    9       itau_w = itau_phy + itap
     9      itau_w = itau_phy + itap + start_time * day_step / iphysiq
    1010c
    1111c Champs 3D:

  • LMDZ5/branches/testing/libf/phylmd/write_histISCCP.h

    r1403 r1665  
    99       ndex3d = 0
    1010c
    11        itau_w = itau_phy + itap
     11       itau_w = itau_phy + itap + start_time * day_step / iphysiq
    1212c
    1313       IF(type_run.EQ."ENSP".OR.type_run.EQ."CLIM") THEN

  • LMDZ5/branches/testing/libf/phylmd/write_histREGDYN.h

    r776 r1665  
    88
    99      ndex3d = 0
    10       itau_w = itau_phy + itap
     10      itau_w = itau_phy + itap + start_time * day_step / iphysiq
    1111c
    1212       CALL histwrite(nid_regdyn,"hw1",itau_w,histoW(:,:,:,1),

  • LMDZ5/branches/testing/libf/phylmd/write_histdayNMC.h

    r1539 r1665  
    55c
    66       ndex3d = 0
    7        itau_w = itau_phy + itap
     7       itau_w = itau_phy + itap + start_time * day_step / iphysiq
    88ccc
    99c  Champs interpolles sur des niveaux de pression du NMC

  • LMDZ5/branches/testing/libf/phylmd/write_histday_seri.h

    r996 r1665  
    77c
    88      ndex2d = 0
    9       itau_w = itau_phy + itap
     9      itau_w = itau_phy + itap + start_time * day_step / iphysiq
    1010c
    1111c Champs 2D:

  • LMDZ5/branches/testing/libf/phylmd/write_histhf3d.h

    r776 r1665  
    77      ndex3d = 0
    88c
    9       itau_w = itau_phy + itap
     9      itau_w = itau_phy + itap + start_time * day_step / iphysiq
    1010c
    1111c Champs 3D:

  • LMDZ5/branches/testing/libf/phylmd/write_histhfNMC.h

    r1539 r1665  
    55c
    66       ndex3d = 0
    7        itau_w = itau_phy + itap
     7       itau_w = itau_phy + itap + start_time * day_step / iphysiq
    88ccc
    99c  Champs interpolles sur des niveaux de pression du NMC

  • LMDZ5/branches/testing/libf/phylmd/write_histmthNMC.h

    r1539 r1665  
    55c
    66       ndex3d = 0
    7        itau_w = itau_phy + itap
     7       itau_w = itau_phy + itap + start_time * day_step / iphysiq
    88ccc
    99c  Champs interpolles sur des niveaux de pression du NMC

  • LMDZ5/branches/testing/libf/phylmd/write_histrac.h

    r1664 r1665  
    55  IF (ecrit_tra > 0.) THEN
    66     
    7      itau_w = itau_phy + nstep
     7     itau_w = itau_phy + nstep + start_time * day_step / iphysiq
    88     
    99     CALL histwrite_phy(nid_tra,.FALSE.,"phis",itau_w,pphis)

  • LMDZ5/branches/testing/libf/phylmd/write_paramLMDZ_phy.h

    r1538 r1665  
    2727c
    2828      ndex2d = 0
    29       itau_w = itau_phy + itap
     29      itau_w = itau_phy + itap + start_time * day_step / iphysiq
    3030c
    3131c Variables globales
Note: See TracChangeset for help on using the changeset viewer.