Changeset 270 for trunk/LMDZ.COMMON/libf


Ignore:
Timestamp:
Aug 18, 2011, 12:09:27 PM (13 years ago)
Author:
emillour
Message:

Ehouarn: Mise a jour des dynamiques (seq et ) pour suivre
les changements et developpements de LMDZ5 terrestre
(mise a niveau avec LMDZ5 trunk, rev 1560). Ce qui ne devrais pas changer grand chose au fonctionnement de base du GCM).
Modification importante: correction du bug sur le cumul des flux de masse pour le transport des traceurs (mais dans les faits semble avoir peu d'impact).
(voir "commit_importants.log" pour les details des ajouts et modifications).

Location:
trunk/LMDZ.COMMON/libf
Files:
2 added
27 edited
4 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/dyn3d/advtrac.F90

    r269 r270  
    1 !
    2 ! $Id: advtrac.F 1446 2010-10-22 09:27:25Z emillour $
    3 !
    4 c
    5 c
    6       SUBROUTINE advtrac(pbaru,pbarv ,
    7      *                   p,  masse,q,iapptrac,teta,
    8      *                  flxw,
    9      *                  pk)
    10 c     Auteur :  F. Hourdin
    11 c
    12 c     Modif. P. Le Van     (20/12/97)
    13 c            F. Codron     (10/99)
    14 c            D. Le Croller (07/2001)
    15 c            M.A Filiberti (04/2002)
    16 c
    17       USE infotrac
    18       USE control_mod
    19  
    20 
    21       IMPLICIT NONE
    22 c
    23 #include "dimensions.h"
    24 #include "paramet.h"
    25 #include "comconst.h"
    26 #include "comvert.h"
    27 #include "comdissip.h"
    28 #include "comgeom2.h"
    29 #include "logic.h"
    30 #include "temps.h"
    31 #include "ener.h"
    32 #include "description.h"
    33 #include "iniprint.h"
    34 
    35 c-------------------------------------------------------------------
    36 c     Arguments
    37 c-------------------------------------------------------------------
    38 c     Ajout PPM
    39 c--------------------------------------------------------
    40       REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
    41 c--------------------------------------------------------
    42       INTEGER iapptrac
    43       REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
    44       REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
    45       REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
    46       REAL pk(ip1jmp1,llm)
    47       REAL flxw(ip1jmp1,llm)
    48 
    49 c-------------------------------------------------------------
    50 c     Variables locales
    51 c-------------------------------------------------------------
    52 
    53       REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
    54       REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
    55       REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
    56       REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
    57       INTEGER iadvtr
    58       INTEGER ij,l,iq,iiq
    59       REAL zdpmin, zdpmax
    60       EXTERNAL  minmax
    61       SAVE iadvtr, massem, pbaruc, pbarvc
    62       DATA iadvtr/0/
    63 c----------------------------------------------------------
    64 c     Rajouts pour PPM
    65 c----------------------------------------------------------
    66       INTEGER indice,n
    67       REAL dtbon ! Pas de temps adaptatif pour que CFL<1
    68       REAL CFLmaxz,aaa,bbb ! CFL maximum
    69       REAL psppm(iim,jjp1) ! pression  au sol
    70       REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
    71       REAL qppm(iim*jjp1,llm,nqtot)
    72       REAL fluxwppm(iim,jjp1,llm)
    73       REAL apppm(llmp1), bpppm(llmp1)
    74       LOGICAL dum,fill
    75       DATA fill/.true./
    76       DATA dum/.true./
    77 
    78       integer,save :: countcfl=0
    79       real cflx(ip1jmp1,llm)
    80       real cfly(ip1jm,llm)
    81       real cflz(ip1jmp1,llm)
    82       real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm)
    83 
    84       IF(iadvtr.EQ.0) THEN
    85          CALL initial0(ijp1llm,pbaruc)
    86          CALL initial0(ijmllm,pbarvc)
    87       ENDIF
    88 
    89 c   accumulation des flux de masse horizontaux
    90       DO l=1,llm
    91          DO ij = 1,ip1jmp1
    92             pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
    93          ENDDO
    94          DO ij = 1,ip1jm
    95             pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
    96          ENDDO
    97       ENDDO
    98 
    99 c   selection de la masse instantannee des mailles avant le transport.
    100       IF(iadvtr.EQ.0) THEN
    101 
    102          CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
    103 ccc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
    104 c
    105       ENDIF
    106 
    107       iadvtr   = iadvtr+1
    108       iapptrac = iadvtr
    109 
    110 
    111 c   Test pour savoir si on advecte a ce pas de temps
    112       IF ( iadvtr.EQ.iapp_tracvl ) THEN
    113 
    114 cc   ..  Modif P.Le Van  ( 20/12/97 )  ....
    115 cc
    116 
    117 c   traitement des flux de masse avant advection.
    118 c     1. calcul de w
    119 c     2. groupement des mailles pres du pole.
    120 
    121         CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
    122 
    123       ! ... Flux de masse diaganostiques traceurs
    124       flxw = wg / REAL(iapp_tracvl)
    125 
    126 c  test sur l'eventuelle creation de valeurs negatives de la masse
    127          DO l=1,llm-1
    128             DO ij = iip2+1,ip1jm
    129               zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l)
    130      s                  - pbarvg(ij-iip1,l) + pbarvg(ij,l)
    131      s                  +       wg(ij,l+1)  - wg(ij,l)
    132             ENDDO
    133             CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
    134             DO ij = iip2,ip1jm
    135                zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
    136             ENDDO
    137 
    138 
    139             CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
    140 
    141             IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
    142             PRINT*,'WARNING DP/P l=',l,'  MIN:',zdpmin,
    143      s        '   MAX:', zdpmax
    144             ENDIF
    145 
    146          ENDDO
    147 
    148 
    149 c-------------------------------------------------------------------
    150 ! Calcul des criteres CFL en X, Y et Z
    151 c-------------------------------------------------------------------
    152 
    153       if (countcfl == 0. ) then
    154           cflxmax(:)=0.
    155           cflymax(:)=0.
    156           cflzmax(:)=0.
    157       endif
    158 
    159       countcfl=countcfl+iapp_tracvl
    160       cflx(:,:)=0.
    161       cfly(:,:)=0.
    162       cflz(:,:)=0.
    163       do l=1,llm
    164          do ij=iip2,ip1jm-1
    165             if (pbarug(ij,l)>=0.) then
    166                 cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l)
    167             else
    168                 cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l)
    169             endif
    170          enddo
    171       enddo
    172       do l=1,llm
    173          do ij=iip2,ip1jm-1,iip1
    174             cflx(ij+iip1,l)=cflx(ij,l)
    175          enddo
    176       enddo
    177 
    178       do l=1,llm
    179          do ij=1,ip1jm
    180             if (pbarvg(ij,l)>=0.) then
    181                 cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l)
    182             else
    183                 cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l)
    184             endif
    185          enddo
    186       enddo
    187 
    188       do l=2,llm
    189          do ij=1,ip1jm
    190             if (wg(ij,l)>=0.) then
    191                 cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l)
    192             else
    193                 cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1)
    194             endif
    195          enddo
    196       enddo
    197 
    198       do l=1,llm
    199          cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l)))
    200          cflymax(l)=max(cflymax(l),maxval(cfly(:,l)))
    201          cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l)))
    202       enddo
     1! $Id: advtrac.F90 1549 2011-07-05 08:41:12Z lguez $
     2
     3SUBROUTINE advtrac(pbaru,pbarv , p,  masse,q,iapptrac,teta, flxw, pk)
     4  !     Auteur :  F. Hourdin
     5  !
     6  !     Modif. P. Le Van     (20/12/97)
     7  !            F. Codron     (10/99)
     8  !            D. Le Croller (07/2001)
     9  !            M.A Filiberti (04/2002)
     10  !
     11  USE infotrac
     12  USE control_mod
     13
     14
     15  IMPLICIT NONE
     16  !
     17  include "dimensions.h"
     18  include "paramet.h"
     19  include "comconst.h"
     20  include "comvert.h"
     21  include "comdissip.h"
     22  include "comgeom2.h"
     23  include "logic.h"
     24  include "temps.h"
     25  include "ener.h"
     26  include "description.h"
     27  include "iniprint.h"
     28
     29  !-------------------------------------------------------------------
     30  !     Arguments
     31  !-------------------------------------------------------------------
     32  !     Ajout PPM
     33  !--------------------------------------------------------
     34  REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
     35  !--------------------------------------------------------
     36  INTEGER iapptrac
     37  REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
     38  REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
     39  REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
     40  REAL pk(ip1jmp1,llm)
     41  REAL flxw(ip1jmp1,llm)
     42
     43  !-------------------------------------------------------------
     44  !     Variables locales
     45  !-------------------------------------------------------------
     46
     47  REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
     48  REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
     49  REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
     50  REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
     51  INTEGER iadvtr
     52  INTEGER ij,l,iq,iiq
     53  REAL zdpmin, zdpmax
     54  EXTERNAL  minmax
     55  SAVE iadvtr, massem, pbaruc, pbarvc
     56  DATA iadvtr/0/
     57  !----------------------------------------------------------
     58  !     Rajouts pour PPM
     59  !----------------------------------------------------------
     60  INTEGER indice,n
     61  REAL dtbon ! Pas de temps adaptatif pour que CFL<1
     62  REAL CFLmaxz,aaa,bbb ! CFL maximum
     63  REAL psppm(iim,jjp1) ! pression  au sol
     64  REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
     65  REAL qppm(iim*jjp1,llm,nqtot)
     66  REAL fluxwppm(iim,jjp1,llm)
     67  REAL apppm(llmp1), bpppm(llmp1)
     68  LOGICAL dum,fill
     69  DATA fill/.true./
     70  DATA dum/.true./
     71
     72  integer,save :: countcfl=0
     73  real cflx(ip1jmp1,llm)
     74  real cfly(ip1jm,llm)
     75  real cflz(ip1jmp1,llm)
     76  real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm)
     77
     78  IF(iadvtr.EQ.0) THEN
     79     CALL initial0(ijp1llm,pbaruc)
     80     CALL initial0(ijmllm,pbarvc)
     81  ENDIF
     82
     83  !   accumulation des flux de masse horizontaux
     84  DO l=1,llm
     85     DO ij = 1,ip1jmp1
     86        pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
     87     ENDDO
     88     DO ij = 1,ip1jm
     89        pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
     90     ENDDO
     91  ENDDO
     92
     93  !   selection de la masse instantannee des mailles avant le transport.
     94  IF(iadvtr.EQ.0) THEN
     95
     96     CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
     97     !cc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
     98     !
     99  ENDIF
     100
     101  iadvtr   = iadvtr+1
     102  iapptrac = iadvtr
     103
     104
     105  !   Test pour savoir si on advecte a ce pas de temps
     106  IF ( iadvtr.EQ.iapp_tracvl ) THEN
     107
     108     !c   ..  Modif P.Le Van  ( 20/12/97 )  ....
     109     !c
     110
     111     !   traitement des flux de masse avant advection.
     112     !     1. calcul de w
     113     !     2. groupement des mailles pres du pole.
     114
     115     CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
     116
     117     ! ... Flux de masse diaganostiques traceurs
     118     flxw = wg / REAL(iapp_tracvl)
     119
     120     !  test sur l'eventuelle creation de valeurs negatives de la masse
     121     DO l=1,llm-1
     122        DO ij = iip2+1,ip1jm
     123           zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l) &
     124                - pbarvg(ij-iip1,l) + pbarvg(ij,l) &
     125                +       wg(ij,l+1)  - wg(ij,l)
     126        ENDDO
     127        CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
     128        DO ij = iip2,ip1jm
     129           zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
     130        ENDDO
     131
     132
     133        CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
     134
     135        IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
     136           PRINT*,'WARNING DP/P l=',l,'  MIN:',zdpmin, &
     137                '   MAX:', zdpmax
     138        ENDIF
     139
     140     ENDDO
     141
     142
     143     !-------------------------------------------------------------------
     144     ! Calcul des criteres CFL en X, Y et Z
     145     !-------------------------------------------------------------------
     146
     147     if (countcfl == 0. ) then
     148        cflxmax(:)=0.
     149        cflymax(:)=0.
     150        cflzmax(:)=0.
     151     endif
     152
     153     countcfl=countcfl+iapp_tracvl
     154     cflx(:,:)=0.
     155     cfly(:,:)=0.
     156     cflz(:,:)=0.
     157     do l=1,llm
     158        do ij=iip2,ip1jm-1
     159           if (pbarug(ij,l)>=0.) then
     160              cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l)
     161           else
     162              cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l)
     163           endif
     164        enddo
     165     enddo
     166     do l=1,llm
     167        do ij=iip2,ip1jm-1,iip1
     168           cflx(ij+iip1,l)=cflx(ij,l)
     169        enddo
     170     enddo
     171
     172     do l=1,llm
     173        do ij=1,ip1jm
     174           if (pbarvg(ij,l)>=0.) then
     175              cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l)
     176           else
     177              cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l)
     178           endif
     179        enddo
     180     enddo
     181
     182     do l=2,llm
     183        do ij=1,ip1jm
     184           if (wg(ij,l)>=0.) then
     185              cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l)
     186           else
     187              cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1)
     188           endif
     189        enddo
     190     enddo
     191
     192     do l=1,llm
     193        cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l)))
     194        cflymax(l)=max(cflymax(l),maxval(cfly(:,l)))
     195        cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l)))
     196     enddo
    203197
    204198!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    205 ! Par defaut, on sort le diagnostic des CFL tous les jours.
    206 ! Si on veut le sortir a chaque pas d'advection en cas de plantage
    207 !     if (countcfl==iapp_tracvl) then
     199     ! Par defaut, on sort le diagnostic des CFL tous les jours.
     200     ! Si on veut le sortir a chaque pas d'advection en cas de plantage
     201     !     if (countcfl==iapp_tracvl) then
    208202!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    209       if (countcfl==day_step) then
    210          do l=1,llm
    211          write(lunout,*) 'L, CFLmax '
    212      s   ,l,maxval(cflx(:,l)),maxval(cfly(:,l)),maxval(cflz(:,l))
    213          enddo
    214          countcfl=0
    215       endif
    216    
    217 c-------------------------------------------------------------------
    218 c   Advection proprement dite (Modification Le Croller (07/2001)
    219 c-------------------------------------------------------------------
    220 
    221 c----------------------------------------------------
    222 c        Calcul des moyennes basées sur la masse
    223 c----------------------------------------------------
    224           call massbar(massem,massebx,masseby)         
    225 
    226 c-----------------------------------------------------------
    227 c     Appel des sous programmes d'advection
    228 c-----------------------------------------------------------
    229       do iq=1,nqtot
    230 c        call clock(t_initial)
     203     if (countcfl==day_step) then
     204        do l=1,llm
     205           write(lunout,*) 'L, CFLmax ' &
     206                ,l,maxval(cflx(:,l)),maxval(cfly(:,l)),maxval(cflz(:,l))
     207        enddo
     208        countcfl=0
     209     endif
     210
     211     !-------------------------------------------------------------------
     212     !   Advection proprement dite (Modification Le Croller (07/2001)
     213     !-------------------------------------------------------------------
     214
     215     !----------------------------------------------------
     216     !        Calcul des moyennes basées sur la masse
     217     !----------------------------------------------------
     218     call massbar(massem,massebx,masseby)         
     219
     220     !-----------------------------------------------------------
     221     !     Appel des sous programmes d'advection
     222     !-----------------------------------------------------------
     223     do iq=1,nqtot
     224        !        call clock(t_initial)
    231225        if(iadv(iq) == 0) cycle
    232 c   ----------------------------------------------------------------
    233 c   Schema de Van Leer I MUSCL
    234 c   ----------------------------------------------------------------
     226        !   ----------------------------------------------------------------
     227        !   Schema de Van Leer I MUSCL
     228        !   ----------------------------------------------------------------
    235229        if(iadv(iq).eq.10) THEN
    236             call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
    237 
    238 c   ----------------------------------------------------------------
    239 c   Schema "pseudo amont" + test sur humidite specifique
    240 C    pour la vapeur d'eau. F. Codron
    241 c   ----------------------------------------------------------------
     230           call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
     231
     232           !   ----------------------------------------------------------------
     233           !   Schema "pseudo amont" + test sur humidite specifique
     234           !    pour la vapeur d'eau. F. Codron
     235           !   ----------------------------------------------------------------
    242236        else if(iadv(iq).eq.14) then
    243 c
    244            CALL vlspltqs( q(1,1,1), 2., massem, wg ,
    245      *                 pbarug,pbarvg,dtvr,p,pk,teta )
    246 c   ----------------------------------------------------------------
    247 c   Schema de Frederic Hourdin
    248 c   ----------------------------------------------------------------
     237           !
     238           CALL vlspltqs( q(1,1,1), 2., massem, wg , &
     239                pbarug,pbarvg,dtvr,p,pk,teta )
     240           !   ----------------------------------------------------------------
     241           !   Schema de Frederic Hourdin
     242           !   ----------------------------------------------------------------
    249243        else if(iadv(iq).eq.12) then
    250 c            Pas de temps adaptatif
     244           !            Pas de temps adaptatif
    251245           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
    252246           if (n.GT.1) then
    253            write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
    254      s             dtvr,'n=',n
     247              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
     248                   dtvr,'n=',n
    255249           endif
    256250           do indice=1,n
    257             call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
     251              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
    258252           end do
    259253        else if(iadv(iq).eq.13) then
    260 c            Pas de temps adaptatif
     254           !            Pas de temps adaptatif
    261255           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
    262256           if (n.GT.1) then
    263            write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
    264      s             dtvr,'n=',n
    265            endif
    266           do indice=1,n
    267             call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
    268           end do
    269 c   ----------------------------------------------------------------
    270 c   Schema de pente SLOPES
    271 c   ----------------------------------------------------------------
     257              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
     258                   dtvr,'n=',n
     259           endif
     260           do indice=1,n
     261              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
     262           end do
     263           !   ----------------------------------------------------------------
     264           !   Schema de pente SLOPES
     265           !   ----------------------------------------------------------------
    272266        else if (iadv(iq).eq.20) then
    273             call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
    274 
    275 c   ----------------------------------------------------------------
    276 c   Schema de Prather
    277 c   ----------------------------------------------------------------
     267           call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
     268
     269           !   ----------------------------------------------------------------
     270           !   Schema de Prather
     271           !   ----------------------------------------------------------------
    278272        else if (iadv(iq).eq.30) then
    279 c            Pas de temps adaptatif
     273           !            Pas de temps adaptatif
    280274           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
    281275           if (n.GT.1) then
    282            write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
    283      s             dtvr,'n=',n
    284            endif
    285            call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
    286      s                     n,dtbon)
    287 
    288 c   ----------------------------------------------------------------
    289 c   Schemas PPM Lin et Rood
    290 c   ----------------------------------------------------------------
    291          else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
    292      s                     iadv(iq).LE.18)) then
    293 
    294 c        Test sur le flux horizontal
    295 c        Pas de temps adaptatif
    296          call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
    297          if (n.GT.1) then
    298          write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
    299      s             dtvr,'n=',n
    300          endif
    301 c        Test sur le flux vertical
    302          CFLmaxz=0.
    303          do l=2,llm
    304            do ij=iip2,ip1jm
    305             aaa=wg(ij,l)*dtvr/massem(ij,l)
    306             CFLmaxz=max(CFLmaxz,aaa)
    307             bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
    308             CFLmaxz=max(CFLmaxz,bbb)
     276              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
     277                   dtvr,'n=',n
     278           endif
     279           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
     280                n,dtbon)
     281
     282           !   ----------------------------------------------------------------
     283           !   Schemas PPM Lin et Rood
     284           !   ----------------------------------------------------------------
     285        else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
     286             iadv(iq).LE.18)) then
     287
     288           !        Test sur le flux horizontal
     289           !        Pas de temps adaptatif
     290           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
     291           if (n.GT.1) then
     292              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
     293                   dtvr,'n=',n
     294           endif
     295           !        Test sur le flux vertical
     296           CFLmaxz=0.
     297           do l=2,llm
     298              do ij=iip2,ip1jm
     299                 aaa=wg(ij,l)*dtvr/massem(ij,l)
     300                 CFLmaxz=max(CFLmaxz,aaa)
     301                 bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
     302                 CFLmaxz=max(CFLmaxz,bbb)
     303              enddo
    309304           enddo
    310          enddo
    311          if (CFLmaxz.GE.1) then
    312             write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
    313          endif
    314 
    315 c-----------------------------------------------------------
    316 c        Ss-prg interface LMDZ.4->PPM3d
    317 c-----------------------------------------------------------
    318 
    319           call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
    320      s                 apppm,bpppm,massebx,masseby,pbarug,pbarvg,
    321      s                 unatppm,vnatppm,psppm)
    322 
    323           do indice=1,n
    324 c---------------------------------------------------------------------
    325 c                         VL (version PPM) horiz. et PPM vert.
    326 c---------------------------------------------------------------------
    327                 if (iadv(iq).eq.11) then
    328 c                  Ss-prg PPM3d de Lin
    329                   call ppm3d(1,qppm(1,1,iq),
    330      s                       psppm,psppm,
    331      s                       unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
    332      s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
    333      s                       fill,dum,220.)
    334 
    335 c----------------------------------------------------------------------
    336 c                           Monotonic PPM
    337 c----------------------------------------------------------------------
    338                else if (iadv(iq).eq.16) then
    339 c                  Ss-prg PPM3d de Lin
    340                   call ppm3d(1,qppm(1,1,iq),
    341      s                       psppm,psppm,
    342      s                       unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
    343      s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
    344      s                       fill,dum,220.)
    345 c---------------------------------------------------------------------
    346 
    347 c---------------------------------------------------------------------
    348 c                           Semi Monotonic PPM
    349 c---------------------------------------------------------------------
    350                else if (iadv(iq).eq.17) then
    351 c                  Ss-prg PPM3d de Lin
    352                   call ppm3d(1,qppm(1,1,iq),
    353      s                       psppm,psppm,
    354      s                       unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
    355      s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
    356      s                       fill,dum,220.)
    357 c---------------------------------------------------------------------
    358 
    359 c---------------------------------------------------------------------
    360 c                         Positive Definite PPM
    361 c---------------------------------------------------------------------
    362                 else if (iadv(iq).eq.18) then
    363 c                  Ss-prg PPM3d de Lin
    364                   call ppm3d(1,qppm(1,1,iq),
    365      s                       psppm,psppm,
    366      s                       unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
    367      s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
    368      s                       fill,dum,220.)
    369 c---------------------------------------------------------------------
    370                 endif
    371             enddo
    372 c-----------------------------------------------------------------
    373 c               Ss-prg interface PPM3d-LMDZ.4
    374 c-----------------------------------------------------------------
    375                   call interpost(q(1,1,iq),qppm(1,1,iq))
    376             endif
    377 c----------------------------------------------------------------------
    378 
    379 c-----------------------------------------------------------------
    380 c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
    381 c et Nord j=1
    382 c-----------------------------------------------------------------
    383 
    384 c                  call traceurpole(q(1,1,iq),massem)
    385 
    386 c calcul du temps cpu pour un schema donne
    387 
    388 c                  call clock(t_final)
    389 cym                  tps_cpu=t_final-t_initial
    390 cym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
    391 
    392        end DO
    393 
    394 
    395 c------------------------------------------------------------------
    396 c   on reinitialise a zero les flux de masse cumules
    397 c---------------------------------------------------
    398           iadvtr=0
    399 
    400        ENDIF ! if iadvtr.EQ.iapp_tracvl
    401 
    402        RETURN
    403        END
    404 
     305           if (CFLmaxz.GE.1) then
     306              write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
     307           endif
     308
     309           !-----------------------------------------------------------
     310           !        Ss-prg interface LMDZ.4->PPM3d
     311           !-----------------------------------------------------------
     312
     313           call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
     314                apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
     315                unatppm,vnatppm,psppm)
     316
     317           do indice=1,n
     318              !----------------------------------------------------------------
     319              !                         VL (version PPM) horiz. et PPM vert.
     320              !----------------------------------------------------------------
     321              if (iadv(iq).eq.11) then
     322                 !                  Ss-prg PPM3d de Lin
     323                 call ppm3d(1,qppm(1,1,iq), &
     324                      psppm,psppm, &
     325                      unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
     326                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
     327                      fill,dum,220.)
     328
     329                 !-------------------------------------------------------------
     330                 !                           Monotonic PPM
     331                 !-------------------------------------------------------------
     332              else if (iadv(iq).eq.16) then
     333                 !                  Ss-prg PPM3d de Lin
     334                 call ppm3d(1,qppm(1,1,iq), &
     335                      psppm,psppm, &
     336                      unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
     337                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
     338                      fill,dum,220.)
     339                 !-------------------------------------------------------------
     340
     341                 !-------------------------------------------------------------
     342                 !                           Semi Monotonic PPM
     343                 !-------------------------------------------------------------
     344              else if (iadv(iq).eq.17) then
     345                 !                  Ss-prg PPM3d de Lin
     346                 call ppm3d(1,qppm(1,1,iq), &
     347                      psppm,psppm, &
     348                      unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
     349                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
     350                      fill,dum,220.)
     351                 !-------------------------------------------------------------
     352
     353                 !-------------------------------------------------------------
     354                 !                         Positive Definite PPM
     355                 !-------------------------------------------------------------
     356              else if (iadv(iq).eq.18) then
     357                 !                  Ss-prg PPM3d de Lin
     358                 call ppm3d(1,qppm(1,1,iq), &
     359                      psppm,psppm, &
     360                      unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
     361                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
     362                      fill,dum,220.)
     363                 !-------------------------------------------------------------
     364              endif
     365           enddo
     366           !-----------------------------------------------------------------
     367           !               Ss-prg interface PPM3d-LMDZ.4
     368           !-----------------------------------------------------------------
     369           call interpost(q(1,1,iq),qppm(1,1,iq))
     370        endif
     371        !----------------------------------------------------------------------
     372
     373        !-----------------------------------------------------------------
     374        ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
     375        ! et Nord j=1
     376        !-----------------------------------------------------------------
     377
     378        !                  call traceurpole(q(1,1,iq),massem)
     379
     380        ! calcul du temps cpu pour un schema donne
     381
     382        !                  call clock(t_final)
     383        !ym                  tps_cpu=t_final-t_initial
     384        !ym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
     385
     386     end DO
     387
     388
     389     !------------------------------------------------------------------
     390     !   on reinitialise a zero les flux de masse cumules
     391     !---------------------------------------------------
     392     iadvtr=0
     393
     394  ENDIF ! if iadvtr.EQ.iapp_tracvl
     395
     396END SUBROUTINE advtrac

  • trunk/LMDZ.COMMON/libf/dyn3d/ce0l.F90

    r101 r270  
    11!
    2 ! $Id: ce0l.F90 1425 2010-09-02 13:45:23Z lguez $
     2! $Id: ce0l.F90 1511 2011-04-28 15:21:47Z jghattas $
    33!
    44!-------------------------------------------------------------------------------
    … …  
    2727#endif
    2828  IMPLICIT NONE
    29 #include "iniprint.h"
    3029#ifndef CPP_EARTH
    3130  WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
    … …  
    3736#include "paramet.h"
    3837#include "indicesol.h"
     38#include "iniprint.h"
    3939#include "temps.h"
    4040#include "logic.h"
     41  INTEGER, PARAMETER            :: longcles=20
     42  REAL,    DIMENSION(longcles)  :: clesphy0
    4143  REAL,    DIMENSION(iip1,jjp1) :: masque
    4244  CHARACTER(LEN=15)             :: calnd
     45  REAL,    DIMENSION(iip1,jjp1) :: phis ! geopotentiel au sol
    4346!-------------------------------------------------------------------------------
    44   CALL conf_gcm( 99, .TRUE. )
     47  CALL conf_gcm( 99, .TRUE. , clesphy0 )
    4548
    4649  CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
    … …  
    7881  WRITE(lunout,'(//)')
    7982  WRITE(lunout,*) ' interbar = ',interbar
    80   CALL etat0_netcdf(interbar,masque,ok_etat0)
     83  CALL etat0_netcdf(interbar,masque,phis,ok_etat0)
    8184
    8285  IF(ok_limit) THEN
    … …  
    8992  END IF
    9093
     94  IF (grilles_gcm_netcdf) THEN
     95     WRITE(lunout,'(//)')
     96     WRITE(lunout,*) '  ***************************  '
     97     WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
     98     WRITE(lunout,*) '  ***************************  '
     99     WRITE(lunout,'(//)')
     100     CALL grilles_gcm_netcdf_sub(masque,phis)
     101  END IF
    91102#endif
    92103! of #ifndef CPP_EARTH #else

  • trunk/LMDZ.COMMON/libf/dyn3d/comdissipn.h

    r1 r270  
    22! $Header$
    33!
    4 c-----------------------------------------------------------------------
    5 c INCLUDE comdissipn.h
     4!  Attention : ce fichier include est compatible format fixe/format libre
     5!                 veillez à n'utiliser que des ! pour les commentaires
     6!                 et à bien positionner les & des lignes de continuation
     7!                 (les placer en colonne 6 et en colonne 73)
     8!-----------------------------------------------------------------------
     9! INCLUDE comdissipn.h
    610
    711      REAL  tetaudiv, tetaurot, tetah, cdivu, crot, cdivh
    8 c
    9       COMMON/comdissipn/ tetaudiv(llm),tetaurot(llm),tetah(llm)   ,
    10      1                        cdivu,      crot,         cdivh
     12!
     13      COMMON/comdissipn/ tetaudiv(llm),tetaurot(llm),tetah(llm)   ,     &
     14     &                        cdivu,      crot,         cdivh
    1115
    12 c
    13 c    Les parametres de ce common proviennent des calculs effectues dans
    14 c             Inidissip  .
    15 c
    16 c-----------------------------------------------------------------------
     16!
     17!    Les parametres de ce common proviennent des calculs effectues dans
     18!             Inidissip  .
     19!
     20!-----------------------------------------------------------------------

  • trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F

    r119 r270  
    230230       CALL getin('output_grads_dyn',output_grads_dyn)
    231231
    232 !Config  Key  = idissip
     232!Config  Key  = dissip_period
    233233!Config  Desc = periode de la dissipation
    234 !Config  Def  = 10
     234!Config  Def  = 0
    235235!Config  Help = periode de la dissipation
    236 !Config         (en pas) ... a completer !
    237        idissip = 10
    238        CALL getin('idissip',idissip)
     236!Config  dissip_period=0 => la valeur sera calcule dans inidissip       
     237!Config  dissip_period>0 => on prend cette valeur
     238       dissip_period = 0
     239       CALL getin('dissip_period',dissip_period)
    239240
    240241ccc  ....   P. Le Van , modif le 29/04/97 .pour la dissipation  ...
    … …  
    652653      write(lunout,*)' periodav = ', periodav
    653654      write(lunout,*)' output_grads_dyn = ', output_grads_dyn
    654       write(lunout,*)' idissip = ', idissip
     655      write(lunout,*)' dissip_period = ', dissip_period
    655656      write(lunout,*)' lstardis = ', lstardis
    656657      write(lunout,*)' nitergdiv = ', nitergdiv
    … …  
    876877      CALL getin('ok_etat0',ok_etat0)
    877878
     879!Config  Key  = grilles_gcm_netcdf
     880!Config  Desc = creation de fichier grilles_gcm.nc dans create_etat0_limit
     881!Config  Def  = n
     882      grilles_gcm_netcdf = .FALSE.
     883      CALL getin('grilles_gcm_netcdf',grilles_gcm_netcdf)
     884
    878885      write(lunout,*)' #########################################'
    879886      write(lunout,*)' Configuration des parametres de cel0'
    … …  
    890897      write(lunout,*)' periodav = ', periodav
    891898      write(lunout,*)' output_grads_dyn = ', output_grads_dyn
    892       write(lunout,*)' idissip = ', idissip
     899      write(lunout,*)' dissip_period = ', dissip_period
    893900      write(lunout,*)' lstardis = ', lstardis
    894901      write(lunout,*)' nitergdiv = ', nitergdiv
    … …  
    922929      write(lunout,*)' ok_limit = ', ok_limit
    923930      write(lunout,*)' ok_etat0 = ', ok_etat0
     931      write(lunout,*)' grilles_gcm_netcdf = ', grilles_gcm_netcdf
    924932c
    925933      RETURN

  • trunk/LMDZ.COMMON/libf/dyn3d/control_mod.F90

    r97 r270  
    1212  REAL    :: periodav
    1313  INTEGER :: nday,day_step,iperiod,iapp_tracvl,nsplit_phys
    14   INTEGER :: iconser,iecri,idissip,iphysiq,iecrimoy
     14  INTEGER :: iconser,iecri,dissip_period,iphysiq,iecrimoy
    1515  INTEGER :: dayref,anneeref, raz_date, ip_ebil_dyn
    1616  LOGICAL :: offline

  • trunk/LMDZ.COMMON/libf/dyn3d/defrun.F

    r6 r270  
    139139
    140140      READ (tapedef,9001) ch1,ch4
    141       READ (tapedef,*)    idissip
    142       WRITE(tapeout,9001) ch1,'idissip'
    143       WRITE(tapeout,*)    idissip
     141      READ (tapedef,*)    dissip_period
     142      WRITE(tapeout,9001) ch1,'dissip_period'
     143      WRITE(tapeout,*)    dissip_period
    144144
    145145ccc  ....   P. Le Van , modif le 29/04/97 .pour la dissipation  ...

  • trunk/LMDZ.COMMON/libf/dyn3d/gcm.F

    r130 r270  
    427427c   -------------------------------
    428428
    429       IF (call_iniphys.and.(iflag_phys.eq.1)) THEN
     429      IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN
    430430         latfi(1)=rlatu(1)
    431431         lonfi(1)=0.
    … …  
    486486#endif
    487487
     488#ifdef CPP_EARTH
     489! Create start file (startphy.nc) and boundary conditions (limit.nc)
     490! for the Earth verstion
     491       if (iflag_phys>=100) then
     492          call iniaqua(ngridmx,latfi,lonfi,iflag_phys)
     493       endif
     494#endif
     495
    488496      if (planet_type.eq."mars") then
    489497! POUR MARS, METTRE UNE FONCTION A PART, genre dynredem0_mars

  • trunk/LMDZ.COMMON/libf/dyn3d/ini_paramLMDZ_dyn.h

    r1 r270  
    7878     .                "once", dt_cum,dt_cum)
    7979c
    80          CALL histdef(nid_ctesGCM, "idissip",
     80         CALL histdef(nid_ctesGCM, "dissip_period",
    8181     .  "periode de la dissipation (en pas) ... a completer",
    8282     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32,

  • trunk/LMDZ.COMMON/libf/dyn3d/iniacademic.F90

    r127 r270  
    11!
    2 ! $Id: iniacademic.F90 1520 2011-05-23 11:37:09Z emillour $
     2! $Id: iniacademic.F90 1529 2011-05-26 15:17:33Z fairhead $
    33!
    44SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
    … …  
    115115  endif
    116116
    117   academic_case: if (iflag_phys == 2) then
     117  academic_case: if (iflag_phys >= 2) then
    118118     ! initializations
    119119
    … …  
    208208     IF (.NOT. read_start) THEN
    209209        ! surface pressure
    210         ps(:)=preff
     210        if (iflag_phys>2) then
     211           ps(:)=108080.  ! Earth aqua/terra planets
     212        else
     213           ps(:)=preff
     214        endif
    211215        ! ground geopotential
    212216        phis(:)=0.

  • trunk/LMDZ.COMMON/libf/dyn3d/iniconst.F

    r127 r270  
    4747c-----------------------------------------------------------------------
    4848
    49       dtdiss  = idissip * dtvr
    5049      dtphys  = iphysiq * dtvr
    5150      unsim   = 1./iim

  • trunk/LMDZ.COMMON/libf/dyn3d/inidissip.F90

    r269 r270  
    11!
    2 ! $Id: inidissip.F 1403 2010-07-01 09:02:53Z fairhead $
     2! $Id: inidissip.F90 1502 2011-03-21 16:07:54Z jghattas $
    33!
    4       SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh  ,
    5      *                       tetagdiv,tetagrot,tetatemp             )
    6 c=======================================================================
    7 c   initialisation de la dissipation horizontale
    8 c=======================================================================
    9 c-----------------------------------------------------------------------
    10 c   declarations:
    11 c   -------------
    12 
    13       USE control_mod
    14 
    15       IMPLICIT NONE
    16 #include "dimensions.h"
    17 #include "paramet.h"
    18 #include "comdissipn.h"
    19 #include "comconst.h"
    20 #include "comvert.h"
    21 #include "logic.h"
    22 
    23       LOGICAL lstardis
    24       INTEGER nitergdiv,nitergrot,niterh
    25       REAL    tetagdiv,tetagrot,tetatemp
    26       REAL fact,zvert(llm),zz
    27       REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)
    28       REAL ullm,vllm,umin,vmin,zhmin,zhmax
    29       REAL zllm,z1llm
    30 
    31       INTEGER l,ij,idum,ii
    32       REAL tetamin
    33       REAL Pup
    34 
    35       REAL ran1
    36 
    37 
    38 c-----------------------------------------------------------------------
    39 c
    40 c   calcul des valeurs propres des operateurs par methode iterrative:
    41 c   -----------------------------------------------------------------
    42 
    43       crot     = -1.
    44       cdivu    = -1.
    45       cdivh    = -1.
    46 
    47 c   calcul de la valeur propre de divgrad:
    48 c   --------------------------------------
    49       idum = 0
    50       DO l = 1, llm
    51        DO ij = 1, ip1jmp1
     4SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh  , &
     5     tetagdiv,tetagrot,tetatemp             )
     6  !=======================================================================
     7  !   initialisation de la dissipation horizontale
     8  !=======================================================================
     9  !-----------------------------------------------------------------------
     10  !   declarations:
     11  !   -------------
     12
     13  USE control_mod, only : dissip_period,iperiod
     14
     15  IMPLICIT NONE
     16  include "dimensions.h"
     17  include "paramet.h"
     18  include "comdissipn.h"
     19  include "comconst.h"
     20  include "comvert.h"
     21  include "logic.h"
     22  include "iniprint.h"
     23
     24  LOGICAL,INTENT(in) :: lstardis
     25  INTEGER,INTENT(in) :: nitergdiv,nitergrot,niterh
     26  REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp
     27
     28! Local variables:
     29  REAL fact,zvert(llm),zz
     30  REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)
     31  REAL ullm,vllm,umin,vmin,zhmin,zhmax
     32  REAL zllm,z1llm
     33
     34  INTEGER l,ij,idum,ii
     35  REAL tetamin
     36  REAL Pup
     37  character (len=80) :: abort_message
     38
     39  REAL ran1
     40
     41
     42  !-----------------------------------------------------------------------
     43  !
     44  !   calcul des valeurs propres des operateurs par methode iterrative:
     45  !   -----------------------------------------------------------------
     46
     47  crot     = -1.
     48  cdivu    = -1.
     49  cdivh    = -1.
     50
     51  !   calcul de la valeur propre de divgrad:
     52  !   --------------------------------------
     53  idum = 0
     54  DO l = 1, llm
     55     DO ij = 1, ip1jmp1
    5256        deltap(ij,l) = 1.
    53        ENDDO
    54       ENDDO
    55 
    56       idum  = -1
    57       zh(1) = RAN1(idum)-.5
    58       idum  = 0
    59       DO ij = 2, ip1jmp1
    60         zh(ij) = RAN1(idum) -.5
    61       ENDDO
    62 
    63       CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)
    64 
    65       CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
    66 
    67       IF ( zhmin .GE. zhmax  )     THEN
    68          PRINT*,'  Inidissip  zh min max  ',zhmin,zhmax
    69          STOP'probleme generateur alleatoire dans inidissip'
    70       ENDIF
    71 
    72       zllm = ABS( zhmax )
    73       DO l = 1,50
    74          IF(lstardis) THEN
    75             CALL divgrad2(1,zh,deltap,niterh,zh)
    76          ELSE
    77             CALL divgrad (1,zh,niterh,zh)
    78          ENDIF
    79 
    80         CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
    81 
    82          zllm  = ABS( zhmax )
    83          z1llm = 1./zllm
    84          DO ij = 1,ip1jmp1
    85             zh(ij) = zh(ij)* z1llm
    86          ENDDO
    87       ENDDO
    88 
    89       IF(lstardis) THEN
    90          cdivh = 1./ zllm
    91       ELSE
    92          cdivh = zllm ** ( -1./niterh )
    93       ENDIF
    94 
    95 c   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)
    96 c   -----------------------------------------------------------------
    97       print*,'calcul des valeurs propres'
    98 
    99       DO  20  ii = 1, 2
    100 c
    101          DO ij = 1, ip1jmp1
    102            zu(ij)  = RAN1(idum) -.5
    103          ENDDO
    104          CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)
    105          DO ij = 1, ip1jm
    106             zv(ij) = RAN1(idum) -.5
    107          ENDDO
    108          CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)
    109 
    110          CALL minmax(iip1*jjp1,zu,umin,ullm )
    111          CALL minmax(iip1*jjm, zv,vmin,vllm )
    112 
    113          ullm = ABS ( ullm )
    114          vllm = ABS ( vllm )
    115 
    116          DO  5  l = 1, 50
    117             IF(ii.EQ.1) THEN
    118 ccccc             CALL covcont( 1,zu,zv,zu,zv )
    119                IF(lstardis) THEN
    120                   CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )
    121                ELSE
    122                   CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )
    123                ENDIF
    124             ELSE
    125                IF(lstardis) THEN
    126                   CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )
    127                ELSE
    128                   CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )
    129                ENDIF
    130             ENDIF
    131 
    132             CALL minmax(iip1*jjp1,zu,umin,ullm )
    133             CALL minmax(iip1*jjm, zv,vmin,vllm )
    134 
    135             ullm = ABS  ( ullm )
    136             vllm = ABS  ( vllm )
    137 
    138             zllm  = MAX( ullm,vllm )
    139             z1llm = 1./ zllm
    140             DO ij = 1, ip1jmp1
    141               zu(ij) = zu(ij)* z1llm
    142             ENDDO
    143             DO ij = 1, ip1jm
    144                zv(ij) = zv(ij)* z1llm
    145             ENDDO
    146  5       CONTINUE
    147 
    148          IF ( ii.EQ.1 ) THEN
    149             IF(lstardis) THEN
    150                cdivu  = 1./zllm
    151             ELSE
    152                cdivu  = zllm **( -1./nitergdiv )
    153             ENDIF
    154          ELSE
    155             IF(lstardis) THEN
    156                crot   = 1./ zllm
    157             ELSE
    158                crot   = zllm **( -1./nitergrot )
    159             ENDIF
    160          ENDIF
    161 
    162  20   CONTINUE
    163 
    164 c   petit test pour les operateurs non star:
    165 c   ----------------------------------------
    166 
    167 c     IF(.NOT.lstardis) THEN
    168          fact    = rad*24./REAL(jjm)
    169          fact    = fact*fact
    170          PRINT*,'coef u ', fact/cdivu, 1./cdivu
    171          PRINT*,'coef r ', fact/crot , 1./crot
    172          PRINT*,'coef h ', fact/cdivh, 1./cdivh
    173 c     ENDIF
    174 
    175 c-----------------------------------------------------------------------
    176 c   variation verticale du coefficient de dissipation:
    177 c   --------------------------------------------------
    178 
    179 c First step: going from 1 to dissip_fac_mid (in gcm.def)
    180 c============
    181       DO l=1,llm
    182          zz      = 1. - preff/presnivs(l)
    183          zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
    184       ENDDO
    185 
    186          write(*,*) 'Dissipation : '
    187          write(*,*) 'Multiplication de la dissipation en altitude :'
    188          write(*,*) '  dissip_fac_mid =', dissip_fac_mid
    189 
    190 c Second step if ok_strato:  from dissip_fac_mid to dissip_fac_up (in gcm.def)
    191 c==========================
    192 c Utilisation de la fonction tangente hyperbolique pour augmenter
    193 c arbitrairement la dissipation et donc la stabilite du modele a
    194 c partir d'une certaine altitude.
    195 
    196 c   Le facteur multiplicatif de basse atmosphere etant deja pris
    197 c   en compte, il faut diviser le facteur multiplicatif de haute
    198 c   atmosphere par celui-ci.
    199 
    200       if (ok_strato) then
    201 
    202        Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
    203        do l=1,llm
    204          zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1)
    205      &   *(1-(0.5*(1+tanh(-6./dissip_deltaz*
    206      &    (-dissip_hdelta*log(presnivs(l)/Pup))  ))))  ))
    207        enddo
    208 
    209          write(*,*) '  dissip_fac_up =', dissip_fac_up
    210          write(*,*) 'Transition mid /up:  Pupstart,delta =',
    211      &             dissip_pupstart,'Pa', dissip_deltaz , '(km)'
    212 
    213       endif
    214 
    215 
    216       PRINT*,'Constantes de temps de la diffusion horizontale'
    217 
    218       tetamin =  1.e+6
    219 
    220       DO l=1,llm
    221         tetaudiv(l)   = zvert(l)/tetagdiv
    222         tetaurot(l)   = zvert(l)/tetagrot
    223         tetah(l)      = zvert(l)/tetatemp
    224 
    225         IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l)
    226         IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l)
    227         IF( tetamin.GT. (1./   tetah(l)) ) tetamin = 1./    tetah(l)
    228       ENDDO
    229 
    230       PRINT *,' INIDI tetamin dtvr ',tetamin,dtvr,iperiod
    231       idissip = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod
    232       PRINT *,' INIDI tetamin idissip ',tetamin,idissip
    233       idissip = MAX(iperiod,idissip)
    234       dtdiss  = idissip * dtvr
    235       PRINT *,' INIDI idissip dtdiss ',idissip,dtdiss
    236 
    237       DO l = 1,llm
    238          PRINT*,zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l),
    239      *                   dtdiss*tetah(l)
    240       ENDDO
    241 
    242 c
    243       RETURN
    244       END
     57     ENDDO
     58  ENDDO
     59
     60  idum  = -1
     61  zh(1) = RAN1(idum)-.5
     62  idum  = 0
     63  DO ij = 2, ip1jmp1
     64     zh(ij) = RAN1(idum) -.5
     65  ENDDO
     66
     67  CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)
     68
     69  CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
     70
     71  IF ( zhmin .GE. zhmax  )     THEN
     72     write(lunout,*)'  Inidissip  zh min max  ',zhmin,zhmax
     73     abort_message='probleme generateur alleatoire dans inidissip'
     74     call abort_gcm('inidissip',abort_message,1)
     75  ENDIF
     76
     77  zllm = ABS( zhmax )
     78  DO l = 1,50
     79     IF(lstardis) THEN
     80        CALL divgrad2(1,zh,deltap,niterh,zh)
     81     ELSE
     82        CALL divgrad (1,zh,niterh,zh)
     83     ENDIF
     84
     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
     92  ENDDO
     93
     94  IF(lstardis) THEN
     95     cdivh = 1./ zllm
     96  ELSE
     97     cdivh = zllm ** ( -1./niterh )
     98  ENDIF
     99
     100  !   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)
     101  !   -----------------------------------------------------------------
     102  write(lunout,*)'inidissip: calcul des valeurs propres'
     103
     104  DO    ii = 1, 2
     105     !
     106     DO ij = 1, ip1jmp1
     107        zu(ij)  = RAN1(idum) -.5
     108     ENDDO
     109     CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)
     110     DO ij = 1, ip1jm
     111        zv(ij) = RAN1(idum) -.5
     112     ENDDO
     113     CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)
     114
     115     CALL minmax(iip1*jjp1,zu,umin,ullm )
     116     CALL minmax(iip1*jjm, zv,vmin,vllm )
     117
     118     ullm = ABS ( ullm )
     119     vllm = ABS ( vllm )
     120
     121     DO    l = 1, 50
     122        IF(ii.EQ.1) THEN
     123           !cccc             CALL covcont( 1,zu,zv,zu,zv )
     124           IF(lstardis) THEN
     125              CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )
     126           ELSE
     127              CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )
     128           ENDIF
     129        ELSE
     130           IF(lstardis) THEN
     131              CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )
     132           ELSE
     133              CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )
     134           ENDIF
     135        ENDIF
     136
     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
     151     end DO
     152
     153     IF ( ii.EQ.1 ) THEN
     154        IF(lstardis) THEN
     155           cdivu  = 1./zllm
     156        ELSE
     157           cdivu  = zllm **( -1./nitergdiv )
     158        ENDIF
     159     ELSE
     160        IF(lstardis) THEN
     161           crot   = 1./ zllm
     162        ELSE
     163           crot   = zllm **( -1./nitergrot )
     164        ENDIF
     165     ENDIF
     166
     167  end DO
     168
     169  !   petit test pour les operateurs non star:
     170  !   ----------------------------------------
     171
     172  !     IF(.NOT.lstardis) THEN
     173  fact    = rad*24./REAL(jjm)
     174  fact    = fact*fact
     175  write(lunout,*)'inidissip: coef u ', fact/cdivu, 1./cdivu
     176  write(lunout,*)'inidissip: coef r ', fact/crot , 1./crot
     177  write(lunout,*)'inidissip: coef h ', fact/cdivh, 1./cdivh
     178  !     ENDIF
     179
     180  !-----------------------------------------------------------------------
     181  !   variation verticale du coefficient de dissipation:
     182  !   --------------------------------------------------
     183
     184! First step: going from 1 to dissip_fac_mid (in gcm.def)
     185!============
     186  DO l=1,llm
     187     zz      = 1. - preff/presnivs(l)
     188     zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
     189  ENDDO
     190
     191  write(lunout,*) 'Dissipation : '
     192  write(lunout,*) 'Multiplication de la dissipation en altitude :'
     193  write(lunout,*) '  dissip_fac_mid =', dissip_fac_mid
     194
     195! Second step if ok_strato:  from dissip_fac_mid to dissip_fac_up (in gcm.def)
     196!==========================
     197! Utilisation de la fonction tangente hyperbolique pour augmenter
     198! arbitrairement la dissipation et donc la stabilite du modele a
     199! partir d'une certaine altitude.
     200
     201!   Le facteur multiplicatif de basse atmosphere etant deja pris
     202!   en compte, il faut diviser le facteur multiplicatif de haute
     203!   atmosphere par celui-ci.
     204
     205  if (ok_strato) then
     206
     207    Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
     208    do l=1,llm
     209      zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1) &
     210                *(1-(0.5*(1+tanh(-6./dissip_deltaz*              &
     211               (-dissip_hdelta*log(presnivs(l)/Pup))  ))))  ))
     212    enddo
     213
     214    write(*,*) '  dissip_fac_up =', dissip_fac_up
     215    write(*,*) 'Transition mid /up:  Pupstart,delta =',           &
     216                   dissip_pupstart,'Pa', dissip_deltaz , '(km)'
     217
     218  endif
     219
     220
     221  write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale'
     222
     223  tetamin =  1.e+6
     224
     225  DO l=1,llm
     226     tetaudiv(l)   = zvert(l)/tetagdiv
     227     tetaurot(l)   = zvert(l)/tetagrot
     228     tetah(l)      = zvert(l)/tetatemp
     229
     230     IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l)
     231     IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l)
     232     IF( tetamin.GT. (1./   tetah(l)) ) tetamin = 1./    tetah(l)
     233  ENDDO
     234
     235  ! If dissip_period=0 calculate value for dissipation period, else keep value read from gcm.def
     236  IF (dissip_period == 0) THEN
     237     dissip_period = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod
     238     write(lunout,*)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ',tetamin,dtvr,iperiod,dissip_period
     239     dissip_period = MAX(iperiod,dissip_period)
     240  END IF
     241
     242  dtdiss  = dissip_period * dtvr
     243  write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr
     244
     245  DO l = 1,llm
     246     write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), &
     247          dtdiss*tetah(l)
     248  ENDDO
     249
     250END SUBROUTINE inidissip

  • trunk/LMDZ.COMMON/libf/dyn3d/integrd.F

    r7 r270  
    11!
    2 ! $Id: integrd.F 1446 2010-10-22 09:27:25Z emillour $
     2! $Id: integrd.F 1550 2011-07-05 09:44:55Z lguez $
    33!
    44      SUBROUTINE integrd
    … …  
    8989        IF( ps(ij).LT.0. ) THEN
    9090         PRINT*,' Au point ij = ',ij, ' , pression sol neg. ', ps(ij)
    91          STOP' dans integrd'
     91         print *, ' dans integrd'
     92         stop 1
    9293        ENDIF
    9394      ENDDO

  • trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F

    r130 r270  
    152152      logical ok_sync
    153153      parameter (ok_sync = .true.)
     154      logical physic
    154155
    155156      data callinigrads/.true./
    … …  
    222223
    223224      itau = 0
     225      physic=.true.
     226      if (iflag_phys==0.or.iflag_phys==2) physic=.false.
     227
    224228c      iday = day_ini+itau/day_step
    225229c      time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
    … …  
    322326      ! Purely Matsuno time stepping
    323327         IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
    324          IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
     328         IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward )
     329     s        apdiss = .TRUE.
    325330         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
    326      s          .and. iflag_phys.EQ.1                 ) apphys = .TRUE.
     331     s          .and. physic                        ) apphys = .TRUE.
    327332      ELSE
    328333      ! Leapfrog/Matsuno time stepping
    329334         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
    330          IF( MOD(itau+1,idissip)  .EQ. 0              ) apdiss = .TRUE.
    331          IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
     335         IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
     336     s        apdiss = .TRUE.
     337         IF( MOD(itau+1,iphysiq).EQ.0.AND.physic       ) apphys=.TRUE.
    332338      END IF
    333339
    … …  
    357363c   -------------------------------------------------------------
    358364
    359       IF( forward. OR . leapf )  THEN
    360 
     365!      IF( forward. OR . leapf )  THEN
     366      IF((.not.forward).OR. leapf )  THEN
     367        ! Ehouarn: gather mass fluxes during backward Matsuno or LF step
    361368         CALL caladvtrac(q,pbaru,pbarv,
    362369     *        p, masse, dq,  teta,
    … …  
    497504          ENDDO
    498505        ENDDO ! of DO l=1,llm
    499           call friction(ucov,vcov,dtvr)
    500506   
    501        if (planet_type.eq."giant") then
    502           ! Intrinsic heat flux
    503           ! Aymeric -- for giant planets
    504           if (ihf .gt. 1.e-6) then
    505           !print *, '**** INTRINSIC HEAT FLUX ****', ihf
    506             DO ij=1,ip1jmp1
    507               teta(ij,1) = teta(ij,1)
    508      &        + dtvr * aire(ij) * ihf / cpp / masse(ij,1)
    509             ENDDO
    510           !print *, '**** d teta '
    511           !print *, dtvr * aire(ijb:ije) * ihf / cpp / masse(ijb:ije,1)
    512           endif
    513        endif
     507        if (planet_type.eq."giant") then
     508          ! add an intrinsic heat flux at the base of the atmosphere
     509          teta(:,1)=teta(:,1)+dtvr*aire(:)*ihf/cpp/masse(:,1)
     510        endif
     511
     512        call friction(ucov,vcov,dtvr)
    514513
    515514   

  • trunk/LMDZ.COMMON/libf/dyn3d/limit_netcdf.F90

    r7 r270  
    11!
    2 ! $Id: limit_netcdf.F90 1441 2010-10-13 13:06:56Z emillour $
     2! $Id: limit_netcdf.F90 1508 2011-04-15 13:05:34Z jghattas $
    33!-------------------------------------------------------------------------------
    44!
    … …  
    4242  REAL, DIMENSION(iip1,jjp1), INTENT(IN) :: masque   ! land mask
    4343#ifndef CPP_EARTH
    44   WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
     44  CALL abort_gcm('limit_netcdf','Earth-specific routine, needs Earth physics',1)
    4545#else
    4646!-------------------------------------------------------------------------------
    … …  
    5252#include "indicesol.h"
    5353
    54 !--- For fractionary sub-cell use (old coding used soil index: 0,1,2,3) --------
    55   LOGICAL, PARAMETER :: fracterre=.TRUE.
    56 
    5754!--- INPUT NETCDF FILES NAMES --------------------------------------------------
    5855  CHARACTER(LEN=25) :: icefile, sstfile, dumstr
    5956  CHARACTER(LEN=25), PARAMETER :: famipsst='amipbc_sst_1x1.nc        ',        &
    6057                                  famipsic='amipbc_sic_1x1.nc        ',        &
    61                                   fclimsst='amipbc_sst_1x1_clim.nc   ',        &
    62                                   fclimsic='amipbc_sic_1x1_clim.nc   ',        &
    6358                                  fcpldsst='cpl_atm_sst.nc           ',        &
    6459                                  fcpldsic='cpl_atm_sic.nc           ',        &
     60                                  fhistsst='histmth_sst.nc           ',        &
     61                                  fhistsic='histmth_sic.nc           ',        &
    6562                                  frugo   ='Rugos.nc                 ',        &
    6663                                  falbe   ='Albedo.nc                '
    67 
     64  CHARACTER(LEN=10) :: varname
    6865!--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------
    6966  REAL,   DIMENSION(klon)                :: fi_ice, verif
    … …  
    8077  INTEGER :: id_FOCE, id_FSIC, id_FTER, id_FLIC
    8178  INTEGER :: NF90_FORMAT
    82   LOGICAL :: lCPL                    !--- T: IPCC-IPSL cpl model output files
    8379  INTEGER :: ndays                   !--- Depending on the output calendar
    8480
    … …  
    104100
    105101!--- RUGOSITY TREATMENT --------------------------------------------------------
    106   WRITE(lunout,*) 'Traitement de la rugosite'
    107   CALL get_2Dfield(frugo,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:))
     102  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la rugosite'
     103  varname='RUGOS'
     104  CALL get_2Dfield(frugo,varname,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:))
    108105
    109106!--- OCEAN TREATMENT -----------------------------------------------------------
    110   PRINT*, 'Traitement de la glace oceanique' ; icefile=''; lCPL=.FALSE.
     107  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la glace oceanique'
    111108
    112109! Input SIC file selection
    113   icefile='fake'
    114   IF(NF90_OPEN(famipsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(famipsic)
    115   IF(NF90_OPEN(fclimsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fclimsic)
    116   IF(NF90_OPEN(fcpldsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fcpldsic)
    117   SELECT CASE(icefile)
    118     CASE(famipsic); dumstr='Amip.'
    119     CASE(fclimsic); dumstr='Amip climatologique.'
    120     CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE.
    121     CASE('fake');   CALL abort_gcm('limit_netcdf','Fichier SIC non reconnu.',1)
    122   END SELECT
     110! Open file only to test if available
     111  IF ( NF90_OPEN(TRIM(famipsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
     112     icefile=TRIM(famipsic)
     113     varname='sicbcs'
     114  ELSE IF( NF90_OPEN(TRIM(fcpldsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
     115     icefile=TRIM(fcpldsic)
     116     varname='SIICECOV'
     117  ELSE IF ( NF90_OPEN(TRIM(fhistsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
     118     icefile=TRIM(fhistsic)
     119     varname='pourc_sic'
     120  ELSE
     121     WRITE(lunout,*) 'ERROR! No sea-ice input file was found.'
     122     WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsic),', ',trim(fcpldsic),', ',trim(fhistsic)
     123     CALL abort_gcm('limit_netcdf','No sea-ice file was found',1)
     124  END IF
    123125  ierr=NF90_CLOSE(nid)
    124   WRITE(lunout,*)'Pour la glace de mer a ete choisi un fichier '//TRIM(dumstr)
    125 
    126   CALL get_2Dfield(icefile,'SIC',interbar,ndays,phy_ice,flag=oldice,lCPL=lCPL)
     126  IF (prt_level>=0) WRITE(lunout,*)'Pour la glace de mer a ete choisi le fichier '//TRIM(icefile)
     127
     128  CALL get_2Dfield(icefile,varname, 'SIC',interbar,ndays,phy_ice,flag=oldice)
    127129
    128130  ALLOCATE(pctsrf_t(klon,nbsrf,ndays))
    129131  DO k=1,ndays
    130     fi_ice=phy_ice(:,k)
    131     WHERE(fi_ice>=1.0  ) fi_ice=1.0
    132     WHERE(fi_ice<EPSFRA) fi_ice=0.0
    133     IF(fracterre) THEN
    134       pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)       ! land soil
    135       pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic)       ! land ice
    136       IF(lCPL) THEN                               ! SIC=pICE*(1-LIC-TER)
    137         pctsrf_t(:,is_sic,k)=fi_ice*(1-pctsrf(:,is_lic)-pctsrf(:,is_ter))
    138       ELSE                                        ! SIC=pICE-LIC
     132     fi_ice=phy_ice(:,k)
     133     WHERE(fi_ice>=1.0  ) fi_ice=1.0
     134     WHERE(fi_ice<EPSFRA) fi_ice=0.0
     135     pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)       ! land soil
     136     pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic)       ! land ice
     137     IF (icefile==trim(fcpldsic)) THEN           ! SIC=pICE*(1-LIC-TER)
     138        pctsrf_t(:,is_sic,k)=fi_ice(:)*(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter))
     139     ELSE IF (icefile==trim(fhistsic)) THEN      ! SIC=pICE
     140        pctsrf_t(:,is_sic,k)=fi_ice(:)
     141     ELSE ! icefile==famipsic                    ! SIC=pICE-LIC
    139142        pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k)
    140       END IF
    141       WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
    142       WHERE(1.0-zmasq<EPSFRA)
     143     END IF
     144     WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
     145     WHERE(1.0-zmasq<EPSFRA)
    143146        pctsrf_t(:,is_sic,k)=0.0
    144147        pctsrf_t(:,is_oce,k)=0.0
    145       ELSEWHERE
     148     ELSEWHERE
    146149        WHERE(pctsrf_t(:,is_sic,k)>=1.0-zmasq)
    147           pctsrf_t(:,is_sic,k)=1.0-zmasq
    148           pctsrf_t(:,is_oce,k)=0.0
     150           pctsrf_t(:,is_sic,k)=1.0-zmasq
     151           pctsrf_t(:,is_oce,k)=0.0
    149152        ELSEWHERE
    150           pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)
    151           WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)
    152             pctsrf_t(:,is_oce,k)=0.0
    153             pctsrf_t(:,is_sic,k)=1.0-zmasq
    154           END WHERE
     153           pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)
     154           WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)
     155              pctsrf_t(:,is_oce,k)=0.0
     156              pctsrf_t(:,is_sic,k)=1.0-zmasq
     157           END WHERE
    155158        END WHERE
    156       END WHERE
    157       nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0)
    158       IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
    159       nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA)
    160       IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad
    161     ELSE
    162       pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)
    163       WHERE(NINT(pctsrf(:,is_ter))==1)
    164         pctsrf_t(:,is_sic,k)=0.
    165         pctsrf_t(:,is_oce,k)=0.                 
    166         WHERE(fi_ice>=1.e-5)
    167           pctsrf_t(:,is_lic,k)=fi_ice
    168           pctsrf_t(:,is_ter,k)=pctsrf_t(:,is_ter,k)-pctsrf_t(:,is_lic,k)
    169         ELSEWHERE
    170           pctsrf_t(:,is_lic,k)=0.0
    171         END WHERE
    172       ELSEWHERE
    173         pctsrf_t(:,is_lic,k) = 0.0
    174         WHERE(fi_ice>=1.e-5)
    175           pctsrf_t(:,is_ter,k)=0.0
    176           pctsrf_t(:,is_sic,k)=fi_ice
    177           pctsrf_t(:,is_oce,k)=1.0-pctsrf_t(:,is_sic,k)
    178         ELSEWHERE
    179           pctsrf_t(:,is_sic,k)=0.0
    180           pctsrf_t(:,is_oce,k)=1.0
    181         END WHERE
    182       END WHERE
    183       verif=sum(pctsrf_t(:,:,k),dim=2)
    184       nbad=COUNT(verif<1.0-1.e-5.OR.verif>1.0+1.e-5)
    185       IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
    186     END IF
     159     END WHERE
     160     nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0)
     161     IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
     162     nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA)
     163     IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad
    187164  END DO
    188165  DEALLOCATE(phy_ice)
    189166
    190167!--- SST TREATMENT -------------------------------------------------------------
    191   WRITE(lunout,*) 'Traitement de la sst' ; sstfile=''; lCPL=.FALSE.
     168  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la sst'
    192169
    193170! Input SST file selection
    194   sstfile='fake'
    195   IF(NF90_OPEN(famipsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(famipsst)
    196   IF(NF90_OPEN(fclimsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fclimsst)
    197   IF(NF90_OPEN(fcpldsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fcpldsst)
    198   SELECT CASE(icefile)
    199     CASE(famipsic); dumstr='Amip.'
    200     CASE(fclimsic); dumstr='Amip climatologique.'
    201     CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE.
    202     CASE('fake');   CALL abort_gcm('limit_netcdf','Fichier SST non reconnu',1)
    203   END SELECT
     171! Open file only to test if available
     172  IF ( NF90_OPEN(TRIM(famipsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
     173     sstfile=TRIM(famipsst)
     174     varname='tosbcs'
     175  ELSE IF ( NF90_OPEN(TRIM(fcpldsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
     176     sstfile=TRIM(fcpldsst)
     177     varname='SISUTESW'
     178  ELSE IF ( NF90_OPEN(TRIM(fhistsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
     179     sstfile=TRIM(fhistsst)
     180     varname='tsol_oce'
     181  ELSE
     182     WRITE(lunout,*) 'ERROR! No sst input file was found.'
     183     WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsst),trim(fcpldsst),trim(fhistsst)
     184     CALL abort_gcm('limit_netcdf','No sst file was found',1)
     185  END IF
    204186  ierr=NF90_CLOSE(nid)
    205   WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(dumstr)
    206 
    207   CALL get_2Dfield(trim(sstfile),'SST',interbar,ndays,phy_sst,flag=extrap)
     187  IF (prt_level>=0) WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(sstfile)
     188
     189  CALL get_2Dfield(sstfile,varname,'SST',interbar,ndays,phy_sst,flag=extrap)
    208190
    209191!--- ALBEDO TREATMENT ----------------------------------------------------------
    210   WRITE(lunout,*) 'Traitement de l albedo'
    211   CALL get_2Dfield(falbe,'ALB',interbar,ndays,phy_alb)
     192  IF (prt_level>1) WRITE(lunout,*) 'Traitement de l albedo'
     193  varname='ALBEDO'
     194  CALL get_2Dfield(falbe,varname,'ALB',interbar,ndays,phy_alb)
    212195
    213196!--- REFERENCE GROUND HEAT FLUX TREATMENT --------------------------------------
    … …  
    215198
    216199!--- OUTPUT FILE WRITING -------------------------------------------------------
    217   WRITE(lunout,*) 'Ecriture du fichier limit'
     200  IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : debut'
    218201
    219202  !--- File creation
    … …  
    264247  ierr=NF90_CLOSE(nid)
    265248
     249  IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : fin'
     250
    266251  DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug)
    267252
    … …  
    276261!-------------------------------------------------------------------------------
    277262!
    278 SUBROUTINE get_2Dfield(fnam, mode, ibar, ndays, champo, flag, mask, lCPL)
     263SUBROUTINE get_2Dfield(fnam, varname, mode, ibar, ndays, champo, flag, mask)
    279264!
    280265!-----------------------------------------------------------------------------
    … …  
    304289! Arguments:
    305290  CHARACTER(LEN=*),  INTENT(IN)     :: fnam     ! NetCDF file name
     291  CHARACTER(LEN=10), INTENT(IN)     :: varname  ! NetCDF variable name
    306292  CHARACTER(LEN=3),  INTENT(IN)     :: mode     ! RUG, SIC, SST or ALB
    307293  LOGICAL,           INTENT(IN)     :: ibar     ! interp on pressure levels
    308294  INTEGER,           INTENT(IN)     :: ndays    ! current year number of days
    309   REAL,    POINTER,  DIMENSION(:, :) :: champo   ! output field = f(t)
     295  REAL,    POINTER,  DIMENSION(:, :) :: champo  ! output field = f(t)
    310296  LOGICAL, OPTIONAL, INTENT(IN)     :: flag     ! extrapol. (SST) old ice (SIC)
    311297  REAL,    OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask
    312   LOGICAL, OPTIONAL, INTENT(IN)     :: lCPL     ! Coupled model flag (for ICE)
    313298!------------------------------------------------------------------------------
    314299! Local variables:
    … …  
    316301  INTEGER :: ncid, varid                  ! NetCDF identifiers
    317302  CHARACTER(LEN=30)               :: dnam       ! dimension name
    318   CHARACTER(LEN=80)               :: varname    ! NetCDF variable name
    319303!--- dimensions
    320304  INTEGER,           DIMENSION(4) :: dids       ! NetCDF dimensions identifiers
    … …  
    331315!--- input files
    332316  CHARACTER(LEN=20)                 :: cal_in   ! calendar
     317  CHARACTER(LEN=20)                 :: unit_sic ! attribute unit in sea-ice file
    333318  INTEGER                           :: ndays_in ! number of days
    334319!--- misc
    … …  
    337322  CHARACTER(LEN=25)                 :: title    ! for messages
    338323  LOGICAL                           :: extrp    ! flag for extrapolation
     324  LOGICAL                           :: oldice   ! flag for old way ice computation
    339325  REAL                              :: chmin, chmax
    340326  INTEGER ierr
    341327  integer n_extrap ! number of extrapolated points
    342328  logical skip
     329
    343330!------------------------------------------------------------------------------
    344331!---Variables depending on keyword 'mode' -------------------------------------
    345332  NULLIFY(champo)
     333
    346334  SELECT CASE(mode)
    347     CASE('RUG'); varname='RUGOS'; title='Rugosite'
    348     CASE('SIC'); varname='sicbcs'; title='Sea-ice'
    349     CASE('SST'); varname='tosbcs'; title='SST'
    350     CASE('ALB'); varname='ALBEDO'; title='Albedo'
     335  CASE('RUG'); title='Rugosite'
     336  CASE('SIC'); title='Sea-ice'
     337  CASE('SST'); title='SST'
     338  CASE('ALB'); title='Albedo'
    351339  END SELECT
     340 
     341
    352342  extrp=.FALSE.
     343  oldice=.FALSE.
    353344  IF ( PRESENT(flag) ) THEN
    354345    IF ( flag .AND. mode=='SST' ) extrp=.TRUE.
     346    IF ( flag .AND. mode=='SIC' ) oldice=.TRUE.
    355347  END IF
    356348
    357349!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
     350  IF (prt_level>5) WRITE(lunout,*) ' Now reading file : ',fnam
    358351  ierr=NF90_OPEN(fnam, NF90_NOWRITE, ncid);             CALL ncerr(ierr, fnam)
    359   ierr=NF90_INQ_VARID(ncid, varname, varid);            CALL ncerr(ierr, fnam)
     352  ierr=NF90_INQ_VARID(ncid, trim(varname), varid);            CALL ncerr(ierr, fnam)
    360353  ierr=NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids); CALL ncerr(ierr, fnam)
     354
     355!--- Read unit for sea-ice variable only
     356  IF (mode=='SIC') THEN
     357     ierr=NF90_GET_ATT(ncid, varid, 'units', unit_sic)
     358     IF(ierr/=NF90_NOERR) THEN
     359        IF (prt_level>5) WRITE(lunout,*) 'No unit was given in sea-ice file. Take percentage as default value'
     360        unit_sic='X'
     361     ELSE
     362        IF (prt_level>5) WRITE(lunout,*) ' Sea-ice cover has unit=',unit_sic
     363     END IF
     364  END IF
    361365
    362366!--- Longitude
    … …  
    365369  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
    366370  ierr=NF90_GET_VAR(ncid, varid, dlon_ini);              CALL ncerr(ierr, fnam)
    367   WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep
     371  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep
    368372
    369373!--- Latitude
    … …  
    372376  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
    373377  ierr=NF90_GET_VAR(ncid, varid, dlat_ini);              CALL ncerr(ierr, fnam)
    374   WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep
     378  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep
    375379
    376380!--- Time (variable is not needed - it is rebuilt - but calendar is)
    … …  
    385389      CASE('SIC', 'SST'); cal_in='gregorian'
    386390    END SELECT
    387     WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' &
     391    IF (prt_level>5) WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' &
    388392         // 'dans '//TRIM(fnam)//'. On choisit la valeur par defaut.'
    389393  END IF
    390   WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', &
     394  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', &
    391395       cal_in
    392396
     397 
    393398!--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION --------------------
    394399  !--- Determining input file number of days, depending on calendar
    … …  
    398403!--- If input records are not monthly, time sampling has to be constant !
    399404  timeyear=mid_months(anneeref, cal_in, lmdep)
    400   IF (lmdep /= 12) WRITE(lunout, '(a, i3, a)') 'Note : les fichiers de ' &
    401        // TRIM(mode) // ' ne comportent pas 12, mais ', lmdep, &
    402        ' enregistrements.'
     405  IF (lmdep /= 12) WRITE(lunout,*) 'Note : les fichiers de ', TRIM(mode), &
     406       ' ne comportent pas 12, mais ', lmdep, ' enregistrements.'
    403407
    404408!--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
    … …  
    406410  IF(extrp) ALLOCATE(work(imdep, jmdep))
    407411
    408   WRITE(lunout, *)
    409   WRITE(lunout, '(a, i3, a)')'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, &
    410        ' CHAMPS.'
     412  IF (prt_level>5) WRITE(lunout, *)
     413  IF (prt_level>5) WRITE(lunout,*)'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, ' CHAMPS.'
    411414  ierr=NF90_INQ_VARID(ncid, varname, varid);             CALL ncerr(ierr, fnam)
    412415  DO l=1, lmdep
    … …  
    419422         work)
    420423
    421     IF(ibar.AND..NOT.(mode=='SIC'.AND.flag)) THEN
    422       IF(l==1) THEN
    423         WRITE(lunout, *)                                                      &
    424   '-------------------------------------------------------------------------'
    425         WRITE(lunout, *)                                                     &
    426   'Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$'
    427         WRITE(lunout, *)                                                      &
    428   '-------------------------------------------------------------------------'
     424    IF(ibar .AND. .NOT.oldice) THEN
     425      IF(l==1 .AND. prt_level>5) THEN
     426        WRITE(lunout, *) '-------------------------------------------------------------------------'
     427        WRITE(lunout, *) 'Utilisation de l''interpolation barycentrique pour ',TRIM(title),' $$$'
     428        WRITE(lunout, *) '-------------------------------------------------------------------------'
    429429      END IF
    430430      IF(mode=='RUG') champ=LOG(champ)
    … …  
    453453
    454454!--- TIME INTERPOLATION ------------------------------------------------------
    455   WRITE(lunout, *)
    456   WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
    457   WRITE(lunout, "(2x, ' Vecteur temps en entree: ', 10f6.1)") timeyear
    458   WRITE(lunout, "(2x, ' Vecteur temps en sortie de 0 a ', i3)") ndays
     455  IF (prt_level>5) THEN
     456     WRITE(lunout, *)
     457     WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
     458     WRITE(lunout, *)' Vecteur temps en entree: ', timeyear
     459     WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays
     460  END IF
     461
    459462  ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays))
    460463  skip = .false.
    … …  
    471474  END DO
    472475  if (n_extrap /= 0) then
    473      print *, "get_2Dfield pchfe_95: n_extrap = ", n_extrap
     476     WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
    474477  end if
    475478  champan(iip1, :, :)=champan(1, :, :)
    … …  
    479482  DO j=1, jjp1
    480483    CALL minmax(iip1, champan(1, j, 10), chmin, chmax)
    481     WRITE(lunout, *)' '//TRIM(title)//' au temps 10 ', chmin, chmax, j
     484    IF (prt_level>5) WRITE(lunout, *)' ',TRIM(title),' au temps 10 ', chmin, chmax, j
    482485  END DO
    483486
    484487!--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
    485488  IF(mode=='SST') THEN
    486     WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38'
     489    IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38'
    487490    WHERE(champan<271.38) champan=271.38
    488491  END IF
    … …  
    490493!--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
    491494  IF(mode=='SIC') THEN
    492     WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
    493     IF(.NOT.lCPL) champan(:, :, :)=champan(:, :, :)/100.
     495    IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
     496
     497    IF (unit_sic=='1') THEN
     498       ! Nothing to be done for sea-ice field is already in fraction of 1
     499       ! This is the case for sea-ice in file cpl_atm_sic.nc
     500       IF (prt_level>5) WRITE(lunout,*) 'Sea-ice field already in fraction of 1'
     501    ELSE
     502       ! Convert sea ice from percentage to fraction of 1
     503       IF (prt_level>5) WRITE(lunout,*) 'Transformt sea-ice field from percentage to fraction of 1.'
     504       champan(:, :, :)=champan(:, :, :)/100.
     505    END IF
     506
    494507    champan(iip1, :, :)=champan(1, :, :)
    495508    WHERE(champan>1.0) champan=1.0
    496509    WHERE(champan<0.0) champan=0.0
    497   END IF
     510 END IF
    498511
    499512!--- DYNAMICAL TO PHYSICAL GRID ----------------------------------------------

  • trunk/LMDZ.COMMON/libf/dyn3d/vlsplt.F

    r109 r270  
    1 !
    2 ! $Header$
    3 !
    4 c
     1c
     2c $Id: vlsplt.F 1550 2011-07-05 09:44:55Z lguez $
    53c
    64
    … …  
    131129      IMPLICIT NONE
    132130c
    133 #include "dimensions.h"
    134 #include "paramet.h"
    135 #include "logic.h"
    136 #include "comvert.h"
    137 #include "comconst.h"
     131      include "dimensions.h"
     132      include "paramet.h"
     133      include "logic.h"
     134      include "comvert.h"
     135      include "comconst.h"
     136      include "iniprint.h"
    138137c
    139138c
    … …  
    353352
    354353      IF(n0.gt.0) THEN
    355       PRINT*,'Nombre de points pour lesquels on advect plus que le'
     354      if (prt_level > 2) PRINT *,
     355     $        'Nombre de points pour lesquels on advect plus que le'
    356356     &       ,'contenu de la maille : ',n0
    357357

  • trunk/LMDZ.COMMON/libf/dyn3d/write_paramLMDZ_dyn.h

    r1 r270  
    5151     .               zx_tmp_2d,iip1*jjp1,ndex2d)
    5252c
    53       zx_tmp_2d(1:iip1,1:jjp1)=REAL(idissip)
    54       CALL histwrite(nid_ctesGCM, "idissip", itau_w,
     53      zx_tmp_2d(1:iip1,1:jjp1)=REAL(dissip_period)
     54      CALL histwrite(nid_ctesGCM, "dissip_period", itau_w,
    5555     .               zx_tmp_2d,iip1*jjp1,ndex2d)
    5656c

  • trunk/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F90

    r269 r270  
    1 !
    2 ! $Id: advtrac_p.F 1446 2010-10-22 09:27:25Z emillour $
    3 !
    4 c
    5 c
    6       SUBROUTINE advtrac_p(pbaru,pbarv ,
    7      *                   p,  masse,q,iapptrac,teta,
    8      *                  flxw,
    9      *                  pk   )
    10 
    11 c     Auteur :  F. Hourdin
    12 c
    13 c     Modif. P. Le Van     (20/12/97)
    14 c            F. Codron     (10/99)
    15 c            D. Le Croller (07/2001)
    16 c            M.A Filiberti (04/2002)
    17 c
    18       USE parallel
    19       USE Write_Field_p
    20       USE Bands
    21       USE mod_hallo
    22       USE Vampir
    23       USE times
    24       USE infotrac
    25       USE control_mod
    26       IMPLICIT NONE
    27 c
    28 #include "dimensions.h"
    29 #include "paramet.h"
    30 #include "comconst.h"
    31 #include "comvert.h"
    32 #include "comdissip.h"
    33 #include "comgeom2.h"
    34 #include "logic.h"
    35 #include "temps.h"
    36 #include "ener.h"
    37 #include "description.h"
    38 
    39 c-------------------------------------------------------------------
    40 c     Arguments
    41 c-------------------------------------------------------------------
    42 c     Ajout PPM
    43 c--------------------------------------------------------
    44       REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
    45 c--------------------------------------------------------
    46       INTEGER iapptrac
    47       REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
    48       REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
    49       REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
    50       REAL pk(ip1jmp1,llm)
    51       REAL               :: flxw(ip1jmp1,llm)
    52 
    53 c-------------------------------------------------------------
    54 c     Variables locales
    55 c-------------------------------------------------------------
    56 
    57       REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
    58       REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
    59       REAL,SAVE::pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
    60       REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
    61       INTEGER iadvtr
    62       INTEGER ij,l,iq,iiq
    63       REAL zdpmin, zdpmax
    64       SAVE iadvtr, massem, pbaruc, pbarvc
    65       DATA iadvtr/0/
    66 c$OMP THREADPRIVATE(iadvtr)
    67 c----------------------------------------------------------
    68 c     Rajouts pour PPM
    69 c----------------------------------------------------------
    70       INTEGER indice,n
    71       REAL dtbon ! Pas de temps adaptatif pour que CFL<1
    72       REAL CFLmaxz,aaa,bbb ! CFL maximum
    73       REAL psppm(iim,jjp1) ! pression  au sol
    74       REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
    75       REAL qppm(iim*jjp1,llm,nqtot)
    76       REAL fluxwppm(iim,jjp1,llm)
    77       REAL apppm(llmp1), bpppm(llmp1)
    78       LOGICAL dum,fill
    79       DATA fill/.true./
    80       DATA dum/.true./
    81       REAL,SAVE :: finmasse(ip1jmp1,llm)
    82       integer ijb,ije,ijb_u,ijb_v,ije_u,ije_v,j
    83       type(Request) :: Request_vanleer
    84       REAL,SAVE :: p_tmp( ip1jmp1,llmp1 )
    85       REAL,SAVE :: teta_tmp(ip1jmp1,llm)
    86       REAL,SAVE :: pk_tmp(ip1jmp1,llm)
    87 
    88       ijb_u=ij_begin
    89       ije_u=ij_end
    90      
    91       ijb_v=ij_begin-iip1
    92       ije_v=ij_end
    93       if (pole_nord) ijb_v=ij_begin
    94       if (pole_sud)  ije_v=ij_end-iip1
    95 
    96       IF(iadvtr.EQ.0) THEN
    97 c         CALL initial0(ijp1llm,pbaruc)
    98 c         CALL initial0(ijmllm,pbarvc)
    99 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
    100         DO l=1,llm   
    101           pbaruc(ijb_u:ije_u,l)=0.
    102           pbarvc(ijb_v:ije_v,l)=0.
    103         ENDDO
    104 c$OMP END DO NOWAIT 
    105       ENDIF
    106 
    107 c   accumulation des flux de masse horizontaux
    108 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    109       DO l=1,llm
    110          DO ij = ijb_u,ije_u
    111             pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
    112          ENDDO
    113          DO ij = ijb_v,ije_v
    114             pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
    115          ENDDO
    116       ENDDO
    117 c$OMP END DO NOWAIT
    118 
    119 c   selection de la masse instantannee des mailles avant le transport.
    120       IF(iadvtr.EQ.0) THEN
    121 
    122 c         CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
    123           ijb=ij_begin
    124           ije=ij_end
    125 
    126 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    127        DO l=1,llm
    128           massem(ijb:ije,l)=masse(ijb:ije,l)
    129        ENDDO
    130 c$OMP END DO NOWAIT
    131 
    132 ccc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
    133 c
    134       ENDIF ! of IF(iadvtr.EQ.0)
    135 
    136       iadvtr   = iadvtr+1
    137 
    138 c$OMP MASTER
    139       iapptrac = iadvtr
    140 c$OMP END MASTER
    141 
    142 c   Test pour savoir si on advecte a ce pas de temps
    143 
    144       IF ( iadvtr.EQ.iapp_tracvl ) THEN
    145 c$OMP MASTER
    146         call suspend_timer(timer_caldyn)
    147 c$OMP END MASTER
    148      
    149       ijb=ij_begin
    150       ije=ij_end
    151      
    152 
    153 cc   ..  Modif P.Le Van  ( 20/12/97 )  ....
    154 cc
    155 
    156 c   traitement des flux de masse avant advection.
    157 c     1. calcul de w
    158 c     2. groupement des mailles pres du pole.
    159 
    160         CALL groupe_p( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
    161 
    162 c$OMP BARRIER
    163 
    164 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    165       DO l=1,llmp1
     1! $Id: advtrac_p.F90 1549 2011-07-05 08:41:12Z lguez $
     2
     3SUBROUTINE advtrac_p(pbaru,pbarv , p,  masse,q,iapptrac,teta, flxw, pk)
     4
     5  !     Auteur :  F. Hourdin
     6  !
     7  !     Modif. P. Le Van     (20/12/97)
     8  !            F. Codron     (10/99)
     9  !            D. Le Croller (07/2001)
     10  !            M.A Filiberti (04/2002)
     11  !
     12  USE parallel
     13  USE Write_Field_p
     14  USE Bands
     15  USE mod_hallo
     16  USE Vampir
     17  USE times
     18  USE infotrac
     19  USE control_mod
     20  IMPLICIT NONE
     21  !
     22  include "dimensions.h"
     23  include "paramet.h"
     24  include "comconst.h"
     25  include "comvert.h"
     26  include "comdissip.h"
     27  include "comgeom2.h"
     28  include "logic.h"
     29  include "temps.h"
     30  include "ener.h"
     31  include "description.h"
     32
     33  !-------------------------------------------------------------------
     34  !     Arguments
     35  !-------------------------------------------------------------------
     36  !     Ajout PPM
     37  !--------------------------------------------------------
     38  REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
     39  !--------------------------------------------------------
     40  INTEGER iapptrac
     41  REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
     42  REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
     43  REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
     44  REAL pk(ip1jmp1,llm)
     45  REAL               :: flxw(ip1jmp1,llm)
     46
     47  !-------------------------------------------------------------
     48  !     Variables locales
     49  !-------------------------------------------------------------
     50
     51  REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
     52  REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
     53  REAL,SAVE::pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
     54  REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
     55  INTEGER iadvtr
     56  INTEGER ij,l,iq,iiq
     57  REAL zdpmin, zdpmax
     58  SAVE iadvtr, massem, pbaruc, pbarvc
     59  DATA iadvtr/0/
     60  !$OMP THREADPRIVATE(iadvtr)
     61  !----------------------------------------------------------
     62  !     Rajouts pour PPM
     63  !----------------------------------------------------------
     64  INTEGER indice,n
     65  REAL dtbon ! Pas de temps adaptatif pour que CFL<1
     66  REAL CFLmaxz,aaa,bbb ! CFL maximum
     67  REAL psppm(iim,jjp1) ! pression  au sol
     68  REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
     69  REAL qppm(iim*jjp1,llm,nqtot)
     70  REAL fluxwppm(iim,jjp1,llm)
     71  REAL apppm(llmp1), bpppm(llmp1)
     72  LOGICAL dum,fill
     73  DATA fill/.true./
     74  DATA dum/.true./
     75  REAL,SAVE :: finmasse(ip1jmp1,llm)
     76  integer ijb,ije,ijb_u,ijb_v,ije_u,ije_v,j
     77  type(Request) :: Request_vanleer
     78  REAL,SAVE :: p_tmp( ip1jmp1,llmp1 )
     79  REAL,SAVE :: teta_tmp(ip1jmp1,llm)
     80  REAL,SAVE :: pk_tmp(ip1jmp1,llm)
     81
     82  ijb_u=ij_begin
     83  ije_u=ij_end
     84
     85  ijb_v=ij_begin-iip1
     86  ije_v=ij_end
     87  if (pole_nord) ijb_v=ij_begin
     88  if (pole_sud)  ije_v=ij_end-iip1
     89
     90  IF(iadvtr.EQ.0) THEN
     91     !         CALL initial0(ijp1llm,pbaruc)
     92     !         CALL initial0(ijmllm,pbarvc)
     93     !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
     94     DO l=1,llm   
     95        pbaruc(ijb_u:ije_u,l)=0.
     96        pbarvc(ijb_v:ije_v,l)=0.
     97     ENDDO
     98     !$OMP END DO NOWAIT 
     99  ENDIF
     100
     101  !   accumulation des flux de masse horizontaux
     102  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     103  DO l=1,llm
     104     DO ij = ijb_u,ije_u
     105        pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
     106     ENDDO
     107     DO ij = ijb_v,ije_v
     108        pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
     109     ENDDO
     110  ENDDO
     111  !$OMP END DO NOWAIT
     112
     113  !   selection de la masse instantannee des mailles avant le transport.
     114  IF(iadvtr.EQ.0) THEN
     115
     116     !         CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
     117     ijb=ij_begin
     118     ije=ij_end
     119
     120     !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     121     DO l=1,llm
     122        massem(ijb:ije,l)=masse(ijb:ije,l)
     123     ENDDO
     124     !$OMP END DO NOWAIT
     125
     126     !cc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
     127     !
     128  ENDIF ! of IF(iadvtr.EQ.0)
     129
     130  iadvtr   = iadvtr+1
     131
     132  !$OMP MASTER
     133  iapptrac = iadvtr
     134  !$OMP END MASTER
     135
     136  !   Test pour savoir si on advecte a ce pas de temps
     137
     138  IF ( iadvtr.EQ.iapp_tracvl ) THEN
     139     !$OMP MASTER
     140     call suspend_timer(timer_caldyn)
     141     !$OMP END MASTER
     142
     143     ijb=ij_begin
     144     ije=ij_end
     145
     146
     147     !c   ..  Modif P.Le Van  ( 20/12/97 )  ....
     148     !c
     149
     150     !   traitement des flux de masse avant advection.
     151     !     1. calcul de w
     152     !     2. groupement des mailles pres du pole.
     153
     154     CALL groupe_p( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
     155
     156     !$OMP BARRIER
     157
     158     !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     159     DO l=1,llmp1
    166160        p_tmp(ijb:ije,l)=p(ijb:ije,l)
    167       ENDDO
    168 c$OMP END DO NOWAIT
    169      
    170 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    171       DO l=1,llm
     161     ENDDO
     162     !$OMP END DO NOWAIT
     163
     164     !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     165     DO l=1,llm
    172166        pk_tmp(ijb:ije,l)=pk(ijb:ije,l)
    173167        teta_tmp(ijb:ije,l)=teta(ijb:ije,l)
    174       ENDDO
    175 c$OMP END DO NOWAIT
    176 
    177 c$OMP MASTER
    178       call VTb(VTHallo)
    179 c$OMP END MASTER
    180 
    181       call Register_SwapFieldHallo(pbarug,pbarug,ip1jmp1,llm,
    182      *                             jj_Nb_vanleer,0,0,Request_vanleer)
    183       call Register_SwapFieldHallo(pbarvg,pbarvg,ip1jm,llm,
    184      *                             jj_Nb_vanleer,1,0,Request_vanleer)
    185       call Register_SwapFieldHallo(massem,massem,ip1jmp1,llm,
    186      *                             jj_Nb_vanleer,0,0,Request_vanleer)
    187       call Register_SwapFieldHallo(wg,wg,ip1jmp1,llm,
    188      *                             jj_Nb_vanleer,0,0,Request_vanleer)
    189       call Register_SwapFieldHallo(teta_tmp,teta_tmp,ip1jmp1,llm,
    190      *                             jj_Nb_vanleer,1,1,Request_vanleer)
    191       call Register_SwapFieldHallo(p_tmp,p_tmp,ip1jmp1,llmp1,
    192      *                             jj_Nb_vanleer,1,1,Request_vanleer)
    193       call Register_SwapFieldHallo(pk_tmp,pk_tmp,ip1jmp1,llm,
    194      *                             jj_Nb_vanleer,1,1,Request_vanleer)
    195       do j=1,nqtot
    196         call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
    197      *                             jj_nb_vanleer,0,0,Request_vanleer)
    198       enddo
    199 
    200       call SendRequest(Request_vanleer)
    201 c$OMP BARRIER
    202       call WaitRequest(Request_vanleer)
    203 
    204 
    205 c$OMP BARRIER
    206 c$OMP MASTER     
    207       call SetDistrib(jj_nb_vanleer)
    208       call VTe(VTHallo)
    209       call VTb(VTadvection)
    210       call start_timer(timer_vanleer)
    211 c$OMP END MASTER
    212 c$OMP BARRIER
    213      
    214       ! ... Flux de masse diaganostiques traceurs
    215          ijb=ij_begin
    216          ije=ij_end
    217          flxw(ijb:ije,1:llm)=wg(ijb:ije,1:llm)/REAL(iapp_tracvl)
    218 
    219 c  test sur l'eventuelle creation de valeurs negatives de la masse
    220          ijb=ij_begin
    221          ije=ij_end
    222          if (pole_nord) ijb=ij_begin+iip1
    223          if (pole_sud) ije=ij_end-iip1
    224          
    225 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
    226          DO l=1,llm-1
    227             DO ij = ijb+1,ije
    228               zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l)
    229      s                  - pbarvg(ij-iip1,l) + pbarvg(ij,l)
    230      s                  +       wg(ij,l+1)  - wg(ij,l)
    231             ENDDO
    232            
    233 c            CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
    234 c ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
    235            
    236             do ij=ijb,ije-iip1+1,iip1
    237               zdp(ij)=zdp(ij+iip1-1)
    238             enddo
    239            
    240             DO ij = ijb,ije
    241                zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
    242             ENDDO
    243 
    244 
    245 c            CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
    246 c  ym ---> eventuellement a revoir
    247             CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
    248            
    249             IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
    250             PRINT*,'WARNING DP/P l=',l,'  MIN:',zdpmin,
    251      s        '   MAX:', zdpmax
    252             ENDIF
    253 
    254          ENDDO
    255 c$OMP END DO NOWAIT
    256 
    257 c-------------------------------------------------------------------
    258 c   Advection proprement dite (Modification Le Croller (07/2001)
    259 c-------------------------------------------------------------------
    260 
    261 c----------------------------------------------------
    262 c        Calcul des moyennes basées sur la masse
    263 c----------------------------------------------------
    264 
    265 cym      ----> Normalement, inutile pour les schémas classiques
    266 cym      ----> Revérifier lors de la parallélisation des autres schemas
    267    
    268 cym          call massbar_p(massem,massebx,masseby) 
    269 
    270           call vlspltgen_p( q,iadv, 2., massem, wg ,
    271      *                    pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
    272 
    273          
    274           GOTO 1234     
    275 c-----------------------------------------------------------
    276 c     Appel des sous programmes d'advection
    277 c-----------------------------------------------------------
    278       do iq=1,nqtot
    279 c        call clock(t_initial)
     168     ENDDO
     169     !$OMP END DO NOWAIT
     170
     171     !$OMP MASTER
     172     call VTb(VTHallo)
     173     !$OMP END MASTER
     174
     175     call Register_SwapFieldHallo(pbarug,pbarug,ip1jmp1,llm, &
     176          jj_Nb_vanleer,0,0,Request_vanleer)
     177     call Register_SwapFieldHallo(pbarvg,pbarvg,ip1jm,llm, &
     178          jj_Nb_vanleer,1,0,Request_vanleer)
     179     call Register_SwapFieldHallo(massem,massem,ip1jmp1,llm, &
     180          jj_Nb_vanleer,0,0,Request_vanleer)
     181     call Register_SwapFieldHallo(wg,wg,ip1jmp1,llm, &
     182          jj_Nb_vanleer,0,0,Request_vanleer)
     183     call Register_SwapFieldHallo(teta_tmp,teta_tmp,ip1jmp1,llm, &
     184          jj_Nb_vanleer,1,1,Request_vanleer)
     185     call Register_SwapFieldHallo(p_tmp,p_tmp,ip1jmp1,llmp1, &
     186          jj_Nb_vanleer,1,1,Request_vanleer)
     187     call Register_SwapFieldHallo(pk_tmp,pk_tmp,ip1jmp1,llm, &
     188          jj_Nb_vanleer,1,1,Request_vanleer)
     189     do j=1,nqtot
     190        call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm, &
     191             jj_nb_vanleer,0,0,Request_vanleer)
     192     enddo
     193
     194     call SendRequest(Request_vanleer)
     195     !$OMP BARRIER
     196     call WaitRequest(Request_vanleer)
     197
     198
     199     !$OMP BARRIER
     200     !$OMP MASTER     
     201     call SetDistrib(jj_nb_vanleer)
     202     call VTe(VTHallo)
     203     call VTb(VTadvection)
     204     call start_timer(timer_vanleer)
     205     !$OMP END MASTER
     206     !$OMP BARRIER
     207
     208     ! ... Flux de masse diaganostiques traceurs
     209     ijb=ij_begin
     210     ije=ij_end
     211     flxw(ijb:ije,1:llm)=wg(ijb:ije,1:llm)/REAL(iapp_tracvl)
     212
     213     !  test sur l'eventuelle creation de valeurs negatives de la masse
     214     ijb=ij_begin
     215     ije=ij_end
     216     if (pole_nord) ijb=ij_begin+iip1
     217     if (pole_sud) ije=ij_end-iip1
     218
     219     !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
     220     DO l=1,llm-1
     221        DO ij = ijb+1,ije
     222           zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l) &
     223                - pbarvg(ij-iip1,l) + pbarvg(ij,l) &
     224                +       wg(ij,l+1)  - wg(ij,l)
     225        ENDDO
     226
     227        !            CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
     228        ! ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
     229
     230        do ij=ijb,ije-iip1+1,iip1
     231           zdp(ij)=zdp(ij+iip1-1)
     232        enddo
     233
     234        DO ij = ijb,ije
     235           zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
     236        ENDDO
     237
     238
     239        !            CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
     240        !  ym ---> eventuellement a revoir
     241        CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
     242
     243        IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
     244           PRINT*,'WARNING DP/P l=',l,'  MIN:',zdpmin, &
     245                '   MAX:', zdpmax
     246        ENDIF
     247
     248     ENDDO
     249     !$OMP END DO NOWAIT
     250
     251     !-------------------------------------------------------------------
     252     !   Advection proprement dite (Modification Le Croller (07/2001)
     253     !-------------------------------------------------------------------
     254
     255     !----------------------------------------------------
     256     !        Calcul des moyennes basées sur la masse
     257     !----------------------------------------------------
     258
     259     !ym      ----> Normalement, inutile pour les schémas classiques
     260     !ym      ----> Revérifier lors de la parallélisation des autres schemas
     261
     262     !ym          call massbar_p(massem,massebx,masseby) 
     263
     264     call vlspltgen_p( q,iadv, 2., massem, wg , &
     265          pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
     266
     267
     268     GOTO 1234     
     269     !-----------------------------------------------------------
     270     !     Appel des sous programmes d'advection
     271     !-----------------------------------------------------------
     272     do iq=1,nqtot
     273        !        call clock(t_initial)
    280274        if(iadv(iq) == 0) cycle
    281 c   ----------------------------------------------------------------
    282 c   Schema de Van Leer I MUSCL
    283 c   ----------------------------------------------------------------
     275        !   ----------------------------------------------------------------
     276        !   Schema de Van Leer I MUSCL
     277        !   ----------------------------------------------------------------
    284278        if(iadv(iq).eq.10) THEN
    285      
    286             call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
    287 
    288 c   ----------------------------------------------------------------
    289 c   Schema "pseudo amont" + test sur humidite specifique
    290 C    pour la vapeur d'eau. F. Codron
    291 c   ----------------------------------------------------------------
     279
     280           call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
     281
     282           !   ----------------------------------------------------------------
     283           !   Schema "pseudo amont" + test sur humidite specifique
     284           !    pour la vapeur d'eau. F. Codron
     285           !   ----------------------------------------------------------------
    292286        else if(iadv(iq).eq.14) then
    293 c
    294 cym           stop 'advtrac : appel à vlspltqs :schema non parallelise'
    295            CALL vlspltqs_p( q(1,1,1), 2., massem, wg ,
    296      *                 pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
    297 c   ----------------------------------------------------------------
    298 c   Schema de Frederic Hourdin
    299 c   ----------------------------------------------------------------
     287           !
     288           !ym        stop 'advtrac : appel à vlspltqs :schema non parallelise'
     289           CALL vlspltqs_p( q(1,1,1), 2., massem, wg , &
     290                pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
     291           !   ----------------------------------------------------------------
     292           !   Schema de Frederic Hourdin
     293           !   ----------------------------------------------------------------
    300294        else if(iadv(iq).eq.12) then
    301           stop 'advtrac : schema non parallelise'
    302 c            Pas de temps adaptatif
     295           stop 'advtrac : schema non parallelise'
     296           !            Pas de temps adaptatif
    303297           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
    304298           if (n.GT.1) then
    305            write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
    306      s             dtvr,'n=',n
     299              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
     300                   dtvr,'n=',n
    307301           endif
    308302           do indice=1,n
    309             call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
     303              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
    310304           end do
    311305        else if(iadv(iq).eq.13) then
    312           stop 'advtrac : schema non parallelise'
    313 c            Pas de temps adaptatif
     306           stop 'advtrac : schema non parallelise'
     307           !            Pas de temps adaptatif
    314308           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
    315309           if (n.GT.1) then
    316            write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
    317      s             dtvr,'n=',n
     310              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
     311                   dtvr,'n=',n
    318312           endif
    319           do indice=1,n
    320             call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
    321           end do
    322 c   ----------------------------------------------------------------
    323 c   Schema de pente SLOPES
    324 c   ----------------------------------------------------------------
     313           do indice=1,n
     314              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
     315           end do
     316           !   ----------------------------------------------------------------
     317           !   Schema de pente SLOPES
     318           !   ----------------------------------------------------------------
    325319        else if (iadv(iq).eq.20) then
    326           stop 'advtrac : schema non parallelise'
    327 
    328             call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
    329 
    330 c   ----------------------------------------------------------------
    331 c   Schema de Prather
    332 c   ----------------------------------------------------------------
     320           stop 'advtrac : schema non parallelise'
     321
     322           call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
     323
     324           !   ----------------------------------------------------------------
     325           !   Schema de Prather
     326           !   ----------------------------------------------------------------
    333327        else if (iadv(iq).eq.30) then
    334           stop 'advtrac : schema non parallelise'
    335 c            Pas de temps adaptatif
     328           stop 'advtrac : schema non parallelise'
     329           !            Pas de temps adaptatif
    336330           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
    337331           if (n.GT.1) then
    338            write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
    339      s             dtvr,'n=',n
     332              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
     333                   dtvr,'n=',n
    340334           endif
    341            call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
    342      s                     n,dtbon)
    343 c   ----------------------------------------------------------------
    344 c   Schemas PPM Lin et Rood
    345 c   ----------------------------------------------------------------
    346        else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
    347      s                     iadv(iq).LE.18)) then
     335           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
     336                n,dtbon)
     337           !   ----------------------------------------------------------------
     338           !   Schemas PPM Lin et Rood
     339           !   ----------------------------------------------------------------
     340        else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
     341             iadv(iq).LE.18)) then
    348342
    349343           stop 'advtrac : schema non parallelise'
    350344
    351 c        Test sur le flux horizontal
    352 c        Pas de temps adaptatif
    353          call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
    354          if (n.GT.1) then
    355            write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
    356      s             dtvr,'n=',n
    357          endif
    358 c        Test sur le flux vertical
    359          CFLmaxz=0.
    360          do l=2,llm
    361            do ij=iip2,ip1jm
    362             aaa=wg(ij,l)*dtvr/massem(ij,l)
    363             CFLmaxz=max(CFLmaxz,aaa)
    364             bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
    365             CFLmaxz=max(CFLmaxz,bbb)
     345           !        Test sur le flux horizontal
     346           !        Pas de temps adaptatif
     347           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
     348           if (n.GT.1) then
     349              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
     350                   dtvr,'n=',n
     351           endif
     352           !        Test sur le flux vertical
     353           CFLmaxz=0.
     354           do l=2,llm
     355              do ij=iip2,ip1jm
     356                 aaa=wg(ij,l)*dtvr/massem(ij,l)
     357                 CFLmaxz=max(CFLmaxz,aaa)
     358                 bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
     359                 CFLmaxz=max(CFLmaxz,bbb)
     360              enddo
    366361           enddo
    367          enddo
    368          if (CFLmaxz.GE.1) then
    369             write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
    370          endif
    371 
    372 c-----------------------------------------------------------
    373 c        Ss-prg interface LMDZ.4->PPM3d
    374 c-----------------------------------------------------------
    375 
    376           call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
    377      s                 apppm,bpppm,massebx,masseby,pbarug,pbarvg,
    378      s                 unatppm,vnatppm,psppm)
    379 
    380           do indice=1,n
    381 c---------------------------------------------------------------------
    382 c                         VL (version PPM) horiz. et PPM vert.
    383 c---------------------------------------------------------------------
    384                 if (iadv(iq).eq.11) then
    385 c                  Ss-prg PPM3d de Lin
    386                   call ppm3d(1,qppm(1,1,iq),
    387      s                       psppm,psppm,
    388      s                       unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
    389      s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
    390      s                       fill,dum,220.)
    391 
    392 c----------------------------------------------------------------------
    393 c                           Monotonic PPM
    394 c----------------------------------------------------------------------
    395                else if (iadv(iq).eq.16) then
    396 c                  Ss-prg PPM3d de Lin
    397                   call ppm3d(1,qppm(1,1,iq),
    398      s                       psppm,psppm,
    399      s                       unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
    400      s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
    401      s                       fill,dum,220.)
    402 c---------------------------------------------------------------------
    403 
    404 c---------------------------------------------------------------------
    405 c                           Semi Monotonic PPM
    406 c---------------------------------------------------------------------
    407                else if (iadv(iq).eq.17) then
    408 c                  Ss-prg PPM3d de Lin
    409                   call ppm3d(1,qppm(1,1,iq),
    410      s                       psppm,psppm,
    411      s                       unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
    412      s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
    413      s                       fill,dum,220.)
    414 c---------------------------------------------------------------------
    415 
    416 c---------------------------------------------------------------------
    417 c                         Positive Definite PPM
    418 c---------------------------------------------------------------------
    419                 else if (iadv(iq).eq.18) then
    420 c                  Ss-prg PPM3d de Lin
    421                   call ppm3d(1,qppm(1,1,iq),
    422      s                       psppm,psppm,
    423      s                       unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
    424      s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
    425      s                       fill,dum,220.)
    426 c---------------------------------------------------------------------
    427                 endif
    428             enddo
    429 c-----------------------------------------------------------------
    430 c               Ss-prg interface PPM3d-LMDZ.4
    431 c-----------------------------------------------------------------
    432                   call interpost(q(1,1,iq),qppm(1,1,iq))
    433             endif
    434 c----------------------------------------------------------------------
    435 
    436 c-----------------------------------------------------------------
    437 c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
    438 c et Nord j=1
    439 c-----------------------------------------------------------------
    440 
    441 c                  call traceurpole(q(1,1,iq),massem)
    442 
    443 c calcul du temps cpu pour un schema donne
    444 
    445 c                  call clock(t_final)
    446 cym                  tps_cpu=t_final-t_initial
    447 cym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
    448 
    449        end DO
    450 
    451 1234  CONTINUE
    452 c$OMP BARRIER
    453 
    454       if (planet_type=="earth") then
     362           if (CFLmaxz.GE.1) then
     363              write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
     364           endif
     365
     366           !-----------------------------------------------------------
     367           !        Ss-prg interface LMDZ.4->PPM3d
     368           !-----------------------------------------------------------
     369
     370           call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
     371                apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
     372                unatppm,vnatppm,psppm)
     373
     374           do indice=1,n
     375              !----------------------------------------------------------------
     376              !                         VL (version PPM) horiz. et PPM vert.
     377              !----------------------------------------------------------------
     378              if (iadv(iq).eq.11) then
     379                 !                  Ss-prg PPM3d de Lin
     380                 call ppm3d(1,qppm(1,1,iq), &
     381                      psppm,psppm, &
     382                      unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
     383                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
     384                      fill,dum,220.)
     385
     386                 !-------------------------------------------------------------
     387                 !                           Monotonic PPM
     388                 !-------------------------------------------------------------
     389              else if (iadv(iq).eq.16) then
     390                 !                  Ss-prg PPM3d de Lin
     391                 call ppm3d(1,qppm(1,1,iq), &
     392                      psppm,psppm, &
     393                      unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
     394                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
     395                      fill,dum,220.)
     396                 !-------------------------------------------------------------
     397
     398                 !-------------------------------------------------------------
     399                 !                           Semi Monotonic PPM
     400                 !-------------------------------------------------------------
     401              else if (iadv(iq).eq.17) then
     402                 !                  Ss-prg PPM3d de Lin
     403                 call ppm3d(1,qppm(1,1,iq), &
     404                      psppm,psppm, &
     405                      unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
     406                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
     407                      fill,dum,220.)
     408                 !-------------------------------------------------------------
     409
     410                 !-------------------------------------------------------------
     411                 !                         Positive Definite PPM
     412                 !-------------------------------------------------------------
     413              else if (iadv(iq).eq.18) then
     414                 !                  Ss-prg PPM3d de Lin
     415                 call ppm3d(1,qppm(1,1,iq), &
     416                      psppm,psppm, &
     417                      unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
     418                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
     419                      fill,dum,220.)
     420                 !-------------------------------------------------------------
     421              endif
     422           enddo
     423           !-----------------------------------------------------------------
     424           !               Ss-prg interface PPM3d-LMDZ.4
     425           !-----------------------------------------------------------------
     426           call interpost(q(1,1,iq),qppm(1,1,iq))
     427        endif
     428        !----------------------------------------------------------------------
     429
     430        !-----------------------------------------------------------------
     431        ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
     432        ! et Nord j=1
     433        !-----------------------------------------------------------------
     434
     435        !                  call traceurpole(q(1,1,iq),massem)
     436
     437        ! calcul du temps cpu pour un schema donne
     438
     439        !                  call clock(t_final)
     440        !ym                  tps_cpu=t_final-t_initial
     441        !ym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
     442
     443     end DO
     444
     4451234 CONTINUE
     446     !$OMP BARRIER
     447
     448     if (planet_type=="earth") then
    455449
    456450        ijb=ij_begin
    457451        ije=ij_end
    458452
    459 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     453        !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    460454        DO l = 1, llm
    461          DO ij = ijb, ije
    462            finmasse(ij,l) =  p(ij,l) - p(ij,l+1)
    463          ENDDO
     455           DO ij = ijb, ije
     456              finmasse(ij,l) =  p(ij,l) - p(ij,l+1)
     457           ENDDO
    464458        ENDDO
    465 c$OMP END DO
     459        !$OMP END DO
    466460
    467461        CALL qminimum_p( q, 2, finmasse )
    468462
    469 c------------------------------------------------------------------
    470 c   on reinitialise a zero les flux de masse cumules
    471 c---------------------------------------------------
    472 c          iadvtr=0
    473 
    474 c$OMP MASTER
     463        !------------------------------------------------------------------
     464        !   on reinitialise a zero les flux de masse cumules
     465        !---------------------------------------------------
     466        !          iadvtr=0
     467
     468        !$OMP MASTER
    475469        call VTe(VTadvection)
    476470        call stop_timer(timer_vanleer)
    477471        call VTb(VThallo)
    478 c$OMP END MASTER
     472        !$OMP END MASTER
    479473
    480474        do j=1,nqtot
    481           call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
    482      *                             jj_nb_caldyn,0,0,Request_vanleer)
     475           call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm, &
     476                jj_nb_caldyn,0,0,Request_vanleer)
    483477        enddo
    484478
    485         call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
    486      *       jj_nb_caldyn,0,0,Request_vanleer)
     479        call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm, &
     480             jj_nb_caldyn,0,0,Request_vanleer)
    487481
    488482        call SendRequest(Request_vanleer)
    489 c$OMP BARRIER
     483        !$OMP BARRIER
    490484        call WaitRequest(Request_vanleer)     
    491485
    492 c$OMP BARRIER
    493 c$OMP MASTER
     486        !$OMP BARRIER
     487        !$OMP MASTER
    494488        call SetDistrib(jj_nb_caldyn)
    495489        call VTe(VThallo)
    496490        call resume_timer(timer_caldyn)
    497 c$OMP END MASTER
    498 c$OMP BARRIER   
    499           iadvtr=0
    500       endif ! of if (planet_type=="earth")
    501        ENDIF ! if iadvtr.EQ.iapp_tracvl
    502 
    503        RETURN
    504        END
    505 
     491 !$OMP END MASTER
     492 !$OMP BARRIER 
     493        iadvtr=0
     494     endif ! of if (planet_type=="earth")
     495  ENDIF ! if iadvtr.EQ.iapp_tracvl
     496
     497END SUBROUTINE advtrac_p

  • trunk/LMDZ.COMMON/libf/dyn3dpar/bilan_dyn_p.F

    r1 r270  
    389389         Q_cum(:,jjb:jje,:,:)=0.
    390390         flux_uQ_cum(:,jjb:jje,:,:)=0.
    391          flux_v_cum(:,jjb:jje,:)=0.
    392391         if (pole_sud) jje=jj_end-1
    393392         flux_v_cum(:,jjb:jje,:)=0.

  • trunk/LMDZ.COMMON/libf/dyn3dpar/ce0l.F90

    r101 r270  
    11!
    2 ! $Id: ce0l.F90 1482 2011-02-09 15:03:09Z jghattas $
     2! $Id: ce0l.F90 1511 2011-04-28 15:21:47Z jghattas $
    33!
    44!-------------------------------------------------------------------------------
    … …  
    1919  USE dimphy
    2020  USE comgeomphy
    21   USE mod_phys_lmdz_para
    22   USE mod_const_mpi
    2321  USE infotrac
    24   USE parallel, ONLY: finalize_parallel
    2522
    2623#ifdef CPP_IOIPSL
    … …  
    3027#endif
    3128  IMPLICIT NONE
    32 #include "iniprint.h"
    3329#ifndef CPP_EARTH
    3430  WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
    … …  
    4036#include "paramet.h"
    4137#include "indicesol.h"
     38#include "iniprint.h"
    4239#include "temps.h"
    4340#include "logic.h"
     41  INTEGER, PARAMETER            :: longcles=20
     42  REAL,    DIMENSION(longcles)  :: clesphy0
    4443  REAL,    DIMENSION(iip1,jjp1) :: masque
    4544  CHARACTER(LEN=15)             :: calnd
     45  REAL,    DIMENSION(iip1,jjp1) :: phis ! geopotentiel au sol
    4646!-------------------------------------------------------------------------------
    47   CALL conf_gcm( 99, .TRUE. )
    48 
    49   CALL init_mpi
     47  CALL conf_gcm( 99, .TRUE. , clesphy0 )
    5048
    5149  CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
    5250  WRITE(lunout,*)'---> klon=',klon
    53   IF (mpi_size>1 .OR. omp_size>1) THEN
    54        CALL abort_gcm('ce0l','In parallel mode,                         &
    55  &                 ce0l must be called only                             &
    56  &                 for 1 process and 1 task',1)
    57   ENDIF
    58 
    5951  CALL InitComgeomphy
    6052
    … …  
    8981  WRITE(lunout,'(//)')
    9082  WRITE(lunout,*) ' interbar = ',interbar
    91   CALL etat0_netcdf(interbar,masque,ok_etat0)
     83  CALL etat0_netcdf(interbar,masque,phis,ok_etat0)
    9284
    9385  IF(ok_limit) THEN
    … …  
    10092  END IF
    10193
    102 !$OMP MASTER
    103   CALL finalize_parallel
    104 !$OMP END MASTER
    105 
     94  IF (grilles_gcm_netcdf) THEN
     95     WRITE(lunout,'(//)')
     96     WRITE(lunout,*) '  ***************************  '
     97     WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
     98     WRITE(lunout,*) '  ***************************  '
     99     WRITE(lunout,'(//)')
     100     CALL grilles_gcm_netcdf_sub(masque,phis)
     101  END IF
    106102#endif
    107103! of #ifndef CPP_EARTH #else

  • trunk/LMDZ.COMMON/libf/dyn3dpar/comdissipn.h

    r1 r270  
    22! $Header$
    33!
    4 c-----------------------------------------------------------------------
    5 c INCLUDE comdissipn.h
     4!  Attention : ce fichier include est compatible format fixe/format libre
     5!                 veillez à n'utiliser que des ! pour les commentaires
     6!                 et à bien positionner les & des lignes de continuation
     7!                 (les placer en colonne 6 et en colonne 73)
     8!-----------------------------------------------------------------------
     9! INCLUDE comdissipn.h
    610
    711      REAL  tetaudiv, tetaurot, tetah, cdivu, crot, cdivh
    8 c
    9       COMMON/comdissipn/ tetaudiv(llm),tetaurot(llm),tetah(llm)   ,
    10      1                        cdivu,      crot,         cdivh
     12!
     13      COMMON/comdissipn/ tetaudiv(llm),tetaurot(llm),tetah(llm)   ,     &
     14     &                        cdivu,      crot,         cdivh
    1115
    12 c
    13 c    Les parametres de ce common proviennent des calculs effectues dans
    14 c             Inidissip  .
    15 c
    16 c-----------------------------------------------------------------------
     16!
     17!    Les parametres de ce common proviennent des calculs effectues dans
     18!             Inidissip  .
     19!
     20!-----------------------------------------------------------------------

  • trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F

    r119 r270  
    11!
    2 ! $Id: conf_gcm.F 1438 2010-10-08 10:19:34Z jghattas $
     2! $Id: conf_gcm.F 1418 2010-07-19 15:11:24Z jghattas $
    33!
    44c
    … …  
    66      SUBROUTINE conf_gcm( tapedef, etatinit )
    77c
     8      USE control_mod
    89#ifdef CPP_IOIPSL
    910      use IOIPSL
    … …  
    1213      use ioipsl_getincom
    1314#endif
    14       use misc_mod
    15       use mod_filtre_fft, ONLY : use_filtre_fft
    16       use mod_hallo, ONLY : use_mpi_alloc
    17       use parallel, ONLY : omp_chunk
    18       USE control_mod
    1915      IMPLICIT NONE
    2016c-----------------------------------------------------------------------
    … …  
    3733#include "serre.h"
    3834#include "comdissnew.h"
    39 !#include "clesphys.h"
    40 #include "iniprint.h"
    4135#include "temps.h"
    4236#include "comconst.h"
    4337
    4438! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
     39! #include "clesphys.h"
     40#include "iniprint.h"
    4541c
    4642c
    … …  
    5349      LOGICAL  fxyhypbb, ysinuss
    5450      INTEGER i
    55      
     51      LOGICAL use_filtre_fft
    5652c
    5753c  -------------------------------------------------------------------
    … …  
    8278c   initialisations:
    8379c   ----------------
    84       adjust=.false.
    85       call getin('adjust',adjust)
    86      
    87       itaumax=0
    88       call getin('itaumax',itaumax);
    89       if (itaumax<=0) itaumax=HUGE(itaumax)
    90      
     80
    9181!Config  Key  = lunout
    9282!Config  Desc = unite de fichier pour les impressions
    … …  
    190180
    191181!Config  Key  = nsplit_phys
    192 !Config  Desc = nombre d'iteration de la physique
    193 !Config  Def  = 240
    194 !Config  Help = nombre d'itration de la physique
    195 !
     182!Config  Desc = nombre de subdivisions par pas physique
     183!Config  Def  = 1
     184!Config  Help = nombre de subdivisions par pas physique
    196185       nsplit_phys = 1
    197186       CALL getin('nsplit_phys',nsplit_phys)
    … …  
    241230       CALL getin('output_grads_dyn',output_grads_dyn)
    242231
    243 !Config  Key  = idissip
     232!Config  Key  = dissip_period
    244233!Config  Desc = periode de la dissipation
    245 !Config  Def  = 10
     234!Config  Def  = 0
    246235!Config  Help = periode de la dissipation
    247 !Config         (en pas) ... a completer !
    248        idissip = 10
    249        CALL getin('idissip',idissip)
     236!Config  dissip_period=0 => la valeur sera calcule dans inidissip       
     237!Config  dissip_period>0 => on prend cette valeur
     238       dissip_period = 0
     239       CALL getin('dissip_period',dissip_period)
    250240
    251241ccc  ....   P. Le Van , modif le 29/04/97 .pour la dissipation  ...
    … …  
    340330       CALL getin('tau_top_bound',tau_top_bound)
    341331
    342 !
    343332!Config  Key  = coefdis
    344333!Config  Desc = coefficient pour gamdissip
    … …  
    407396       ip_ebil_dyn = 0
    408397       CALL getin('ip_ebil_dyn',ip_ebil_dyn)
    409 !
    410 
    411398
    412399ccc  ....   P. Le Van , ajout  le 7/03/95 .pour le zoom ...
    … …  
    613600       offline = .FALSE.
    614601       CALL getin('offline',offline)
    615        IF (offline .AND. adjust) THEN
    616           WRITE(lunout,*)
    617      &         'WARNING : option offline does not work with adjust=y :'
    618           WRITE(lunout,*) 'the files defstoke.nc, fluxstoke.nc ',
    619      &         'and fluxstokev.nc will not be created'
    620           WRITE(lunout,*)
    621      &         'only the file phystoke.nc will still be created '
    622        END IF
    623        
     602
    624603!Config  Key  = config_inca
    625604!Config  Desc = Choix de configuration de INCA
    … …  
    655634      ok_dyn_ave = .FALSE.
    656635      CALL getin('ok_dyn_ave',ok_dyn_ave)
     636
    657637
    658638      write(lunout,*)' #########################################'
    … …  
    669649      write(lunout,*)' day_step = ', day_step
    670650      write(lunout,*)' iperiod = ', iperiod
    671       write(lunout,*)' nsplit_phys = ', nsplit_phys
    672651      write(lunout,*)' iconser = ', iconser
    673652      write(lunout,*)' iecri = ', iecri
    674653      write(lunout,*)' periodav = ', periodav
    675654      write(lunout,*)' output_grads_dyn = ', output_grads_dyn
    676       write(lunout,*)' idissip = ', idissip
     655      write(lunout,*)' dissip_period = ', dissip_period
    677656      write(lunout,*)' lstardis = ', lstardis
    678657      write(lunout,*)' nitergdiv = ', nitergdiv
    … …  
    816795       offline = .FALSE.
    817796       CALL getin('offline',offline)
    818        IF (offline .AND. adjust) THEN
    819           WRITE(lunout,*)
    820      &         'WARNING : option offline does not work with adjust=y :'
    821           WRITE(lunout,*) 'the files defstoke.nc, fluxstoke.nc ',
    822      &         'and fluxstokev.nc will not be created'
    823           WRITE(lunout,*)
    824      &         'only the file phystoke.nc will still be created '
    825        END IF
    826797
    827798!Config  Key  = config_inca
    … …  
    836807
    837808!Config  Key  = ok_dynzon
    838 !Config  Desc = calcul et sortie des transports
     809!Config  Desc = sortie des transports zonaux dans la dynamique
    839810!Config  Def  = n
    840 !Config  Help = Permet de mettre en route le calcul des transports
     811!Config  Help =
    841812!Config         
    842       ok_dynzon = .FALSE.
    843       CALL getin('ok_dynzon',ok_dynzon)
     813       ok_dynzon = .FALSE.
     814       CALL getin('ok_dynzon',ok_dynzon)
    844815
    845816!Config  Key  = ok_dyn_ins
    … …  
    863834!Config  Def  = false
    864835!Config  Help = permet d'activer l'utilisation des FFT pour effectuer
    865 !Config         le filtrage aux poles.
     836!Config         le filtrage aux poles.
     837! Le filtre fft n'est pas implemente dans dyn3d
    866838      use_filtre_fft=.FALSE.
    867839      CALL getin('use_filtre_fft',use_filtre_fft)
    868840
    869       IF (use_filtre_fft .AND. grossismx /= 1.0) THEN
    870         write(lunout,*)'WARNING !!! '
    871         write(lunout,*)"Le zoom en longitude est incompatible",
    872      &                 " avec l'utilisation du filtre FFT ",
    873      &                 "---> filtre FFT désactivé "
    874        use_filtre_fft=.FALSE.
     841      IF (use_filtre_fft) THEN
     842        write(lunout,*)'STOP !!!'
     843        write(lunout,*)'use_filtre_fft n est pas implemente dans dyn3d'
     844        STOP
    875845      ENDIF
    876846     
    877  
    878      
    879 !Config  Key  = use_mpi_alloc
    880 !Config  Desc = Utilise un buffer MPI en m�moire globale
    881 !Config  Def  = false
    882 !Config  Help = permet d'activer l'utilisation d'un buffer MPI
    883 !Config         en m�moire globale a l'aide de la fonction MPI_ALLOC.
    884 !Config         Cela peut am�liorer la bande passante des transferts MPI
    885 !Config         d'un facteur 2 
    886       use_mpi_alloc=.FALSE.
    887       CALL getin('use_mpi_alloc',use_mpi_alloc)
    888 
    889 !Config  Key  = omp_chunk
    890 !Config  Desc = taille des blocs openmp
    891 !Config  Def  = 1
    892 !Config  Help = defini la taille des packets d'it�ration openmp
    893 !Config         distribu�e � chaque t�che lors de l'entr�e dans une
    894 !Config         boucle parall�lis�e
    895  
    896       omp_chunk=1
    897       CALL getin('omp_chunk',omp_chunk)
    898 
    899847!Config key = ok_strato
    900848!Config  Desc = activation de la version strato
    … …  
    902850!Config  Help = active la version stratosphérique de LMDZ de F. Lott
    903851
    904       ok_strato=.FALSE.
     852      ok_strato=.TRUE.
    905853      CALL getin('ok_strato',ok_strato)
    906854
    … …  
    928876      ok_etat0 = .TRUE.
    929877      CALL getin('ok_etat0',ok_etat0)
     878
     879!Config  Key  = grilles_gcm_netcdf
     880!Config  Desc = creation de fichier grilles_gcm.nc dans create_etat0_limit
     881!Config  Def  = n
     882      grilles_gcm_netcdf = .FALSE.
     883      CALL getin('grilles_gcm_netcdf',grilles_gcm_netcdf)
    930884
    931885      write(lunout,*)' #########################################'
    … …  
    943897      write(lunout,*)' periodav = ', periodav
    944898      write(lunout,*)' output_grads_dyn = ', output_grads_dyn
    945       write(lunout,*)' idissip = ', idissip
     899      write(lunout,*)' dissip_period = ', dissip_period
    946900      write(lunout,*)' lstardis = ', lstardis
    947901      write(lunout,*)' nitergdiv = ', nitergdiv
    … …  
    968922      write(lunout,*)' offline = ', offline
    969923      write(lunout,*)' config_inca = ', config_inca
    970       write(lunout,*)' ok_dynzon = ', ok_dynzon 
     924      write(lunout,*)' ok_dynzon = ', ok_dynzon
    971925      write(lunout,*)' ok_dyn_ins = ', ok_dyn_ins
    972926      write(lunout,*)' ok_dyn_ave = ', ok_dyn_ave
    973       write(lunout,*)' use_filtre_fft = ', use_filtre_fft
    974       write(lunout,*)' use_mpi_alloc = ', use_mpi_alloc
    975       write(lunout,*)' omp_chunk = ', omp_chunk
    976927      write(lunout,*)' ok_strato = ', ok_strato
    977928      write(lunout,*)' ok_gradsfile = ', ok_gradsfile
    978929      write(lunout,*)' ok_limit = ', ok_limit
    979930      write(lunout,*)' ok_etat0 = ', ok_etat0
     931      write(lunout,*)' grilles_gcm_netcdf = ', grilles_gcm_netcdf
    980932c
    981933      RETURN

  • trunk/LMDZ.COMMON/libf/dyn3dpar/control_mod.F90

    r97 r270  
    1212  REAL    :: periodav
    1313  INTEGER :: nday,day_step,iperiod,iapp_tracvl,nsplit_phys
    14   INTEGER :: iconser,iecri,idissip,iphysiq,iecrimoy
     14  INTEGER :: iconser,iecri,dissip_period,iphysiq,iecrimoy
    1515  INTEGER :: dayref,anneeref, raz_date, ip_ebil_dyn
    1616  LOGICAL :: offline

  • trunk/LMDZ.COMMON/libf/dyn3dpar/defrun.F

    r1 r270  
    66      SUBROUTINE defrun( tapedef, etatinit, clesphy0 )
    77c
     8! ========================== ATTENTION =============================
     9! COMMENTAIRE SL :
     10! NE SERT PLUS APPAREMMENT
     11! DONC PAS MIS A JOUR POUR L'UTILISATION AVEC LES PLANETES
     12! ==================================================================
     13
    814      USE control_mod
     15 
    916      IMPLICIT NONE
    1017c-----------------------------------------------------------------------
    … …  
    132139
    133140      READ (tapedef,9001) ch1,ch4
    134       READ (tapedef,*)    idissip
    135       WRITE(tapeout,9001) ch1,'idissip'
    136       WRITE(tapeout,*)    idissip
     141      READ (tapedef,*)    dissip_period
     142      WRITE(tapeout,9001) ch1,'dissip_period'
     143      WRITE(tapeout,*)    dissip_period
    137144
    138145ccc  ....   P. Le Van , modif le 29/04/97 .pour la dissipation  ...
    … …  
    196203
    197204
    198 ccc   .... P.Le Van, ajout le 03/01/96 pour l'ecriture phys ...
    199 c
    200205      READ (tapedef,9001) ch1,ch4
    201206      READ (tapedef,*)    cycle_diurne

  • trunk/LMDZ.COMMON/libf/dyn3dpar/exner_hyb_p.F

    r127 r270  
    125125      endif
    126126!$OMP END MASTER
    127 
     127!$OMP BARRIER
    128128        jjb=jj_begin
    129129        jje=jj_end
    … …  
    171171      endif
    172172c$OMP END MASTER
     173c$OMP BARRIER
    173174c
    174175c

  • trunk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu_p.F

    r127 r270  
    121121      endif
    122122!$OMP END MASTER
    123 
     123!$OMP BARRIER
    124124        jjb=jj_begin
    125125        jje=jj_end
    … …  
    169169      endif
    170170c$OMP END MASTER
     171c$OMP BARRIER
    171172c
    172173c

  • trunk/LMDZ.COMMON/libf/dyn3dpar/iniacademic.F90

    r127 r270  
    11!
    2 ! $Id: iniacademic.F90 1520 2011-05-23 11:37:09Z emillour $
     2! $Id: iniacademic.F90 1529 2011-05-26 15:17:33Z fairhead $
    33!
    44SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
    … …  
    115115  endif
    116116
    117   academic_case: if (iflag_phys == 2) then
     117  academic_case: if (iflag_phys >= 2) then
    118118     ! initializations
    119119
    … …  
    208208     IF (.NOT. read_start) THEN
    209209        ! surface pressure
    210         ps(:)=preff
     210        if (iflag_phys>2) then
     211           ps(:)=108080.  ! Earth aqua/terra planets
     212        else
     213           ps(:)=preff
     214        endif
    211215        ! ground geopotential
    212216        phis(:)=0.

  • trunk/LMDZ.COMMON/libf/dyn3dpar/iniconst.F

    r127 r270  
    4747c-----------------------------------------------------------------------
    4848
    49       dtdiss  = idissip * dtvr
    5049      dtphys  = iphysiq * dtvr
    5150      unsim   = 1./iim

  • trunk/LMDZ.COMMON/libf/dyn3dpar/inidissip.F90

    r269 r270  
    11!
    2 ! $Id: inidissip.F 1403 2010-07-01 09:02:53Z fairhead $
     2! $Id: inidissip.F90 1502 2011-03-21 16:07:54Z jghattas $
    33!
    4       SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh  ,
    5      *                       tetagdiv,tetagrot,tetatemp             )
    6 c=======================================================================
    7 c   initialisation de la dissipation horizontale
    8 c=======================================================================
    9 c-----------------------------------------------------------------------
    10 c   declarations:
    11 c   -------------
    12 
    13       USE control_mod
    14 
    15       IMPLICIT NONE
    16 #include "dimensions.h"
    17 #include "paramet.h"
    18 #include "comdissipn.h"
    19 #include "comconst.h"
    20 #include "comvert.h"
    21 #include "logic.h"
    22 
    23       LOGICAL lstardis
    24       INTEGER nitergdiv,nitergrot,niterh
    25       REAL    tetagdiv,tetagrot,tetatemp
    26       REAL fact,zvert(llm),zz
    27       REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)
    28       REAL ullm,vllm,umin,vmin,zhmin,zhmax
    29       REAL zllm,z1llm
    30 
    31       INTEGER l,ij,idum,ii
    32       REAL tetamin
    33       REAL Pup
    34 
    35       REAL ran1
    36 
    37 
    38 c-----------------------------------------------------------------------
    39 c
    40 c   calcul des valeurs propres des operateurs par methode iterrative:
    41 c   -----------------------------------------------------------------
    42 
    43       crot     = -1.
    44       cdivu    = -1.
    45       cdivh    = -1.
    46 
    47 c   calcul de la valeur propre de divgrad:
    48 c   --------------------------------------
    49       idum = 0
    50       DO l = 1, llm
    51        DO ij = 1, ip1jmp1
     4SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh  , &
     5     tetagdiv,tetagrot,tetatemp             )
     6  !=======================================================================
     7  !   initialisation de la dissipation horizontale
     8  !=======================================================================
     9  !-----------------------------------------------------------------------
     10  !   declarations:
     11  !   -------------
     12
     13  USE control_mod, only : dissip_period,iperiod
     14
     15  IMPLICIT NONE
     16  include "dimensions.h"
     17  include "paramet.h"
     18  include "comdissipn.h"
     19  include "comconst.h"
     20  include "comvert.h"
     21  include "logic.h"
     22  include "iniprint.h"
     23
     24  LOGICAL,INTENT(in) :: lstardis
     25  INTEGER,INTENT(in) :: nitergdiv,nitergrot,niterh
     26  REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp
     27
     28! Local variables:
     29  REAL fact,zvert(llm),zz
     30  REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm)
     31  REAL ullm,vllm,umin,vmin,zhmin,zhmax
     32  REAL zllm,z1llm
     33
     34  INTEGER l,ij,idum,ii
     35  REAL tetamin
     36  REAL Pup
     37  character (len=80) :: abort_message
     38
     39  REAL ran1
     40
     41
     42  !-----------------------------------------------------------------------
     43  !
     44  !   calcul des valeurs propres des operateurs par methode iterrative:
     45  !   -----------------------------------------------------------------
     46
     47  crot     = -1.
     48  cdivu    = -1.
     49  cdivh    = -1.
     50
     51  !   calcul de la valeur propre de divgrad:
     52  !   --------------------------------------
     53  idum = 0
     54  DO l = 1, llm
     55     DO ij = 1, ip1jmp1
    5256        deltap(ij,l) = 1.
    53        ENDDO
    54       ENDDO
    55 
    56       idum  = -1
    57       zh(1) = RAN1(idum)-.5
    58       idum  = 0
    59       DO ij = 2, ip1jmp1
    60         zh(ij) = RAN1(idum) -.5
    61       ENDDO
    62 
    63       CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)
    64 
    65       CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
    66 
    67       IF ( zhmin .GE. zhmax  )     THEN
    68          PRINT*,'  Inidissip  zh min max  ',zhmin,zhmax
    69          STOP'probleme generateur alleatoire dans inidissip'
    70       ENDIF
    71 
    72       zllm = ABS( zhmax )
    73       DO l = 1,50
    74          IF(lstardis) THEN
    75             CALL divgrad2(1,zh,deltap,niterh,zh)
    76          ELSE
    77             CALL divgrad (1,zh,niterh,zh)
    78          ENDIF
    79 
    80         CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
    81 
    82          zllm  = ABS( zhmax )
    83          z1llm = 1./zllm
    84          DO ij = 1,ip1jmp1
    85             zh(ij) = zh(ij)* z1llm
    86          ENDDO
    87       ENDDO
    88 
    89       IF(lstardis) THEN
    90          cdivh = 1./ zllm
    91       ELSE
    92          cdivh = zllm ** ( -1./niterh )
    93       ENDIF
    94 
    95 c   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)
    96 c   -----------------------------------------------------------------
    97       print*,'calcul des valeurs propres'
    98 
    99       DO  20  ii = 1, 2
    100 c
    101          DO ij = 1, ip1jmp1
    102            zu(ij)  = RAN1(idum) -.5
    103          ENDDO
    104          CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)
    105          DO ij = 1, ip1jm
    106             zv(ij) = RAN1(idum) -.5
    107          ENDDO
    108          CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)
    109 
    110          CALL minmax(iip1*jjp1,zu,umin,ullm )
    111          CALL minmax(iip1*jjm, zv,vmin,vllm )
    112 
    113          ullm = ABS ( ullm )
    114          vllm = ABS ( vllm )
    115 
    116          DO  5  l = 1, 50
    117             IF(ii.EQ.1) THEN
    118 ccccc             CALL covcont( 1,zu,zv,zu,zv )
    119                IF(lstardis) THEN
    120                   CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )
    121                ELSE
    122                   CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )
    123                ENDIF
    124             ELSE
    125                IF(lstardis) THEN
    126                   CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )
    127                ELSE
    128                   CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )
    129                ENDIF
    130             ENDIF
    131 
    132             CALL minmax(iip1*jjp1,zu,umin,ullm )
    133             CALL minmax(iip1*jjm, zv,vmin,vllm )
    134 
    135             ullm = ABS  ( ullm )
    136             vllm = ABS  ( vllm )
    137 
    138             zllm  = MAX( ullm,vllm )
    139             z1llm = 1./ zllm
    140             DO ij = 1, ip1jmp1
    141               zu(ij) = zu(ij)* z1llm
    142             ENDDO
    143             DO ij = 1, ip1jm
    144                zv(ij) = zv(ij)* z1llm
    145             ENDDO
    146  5       CONTINUE
    147 
    148          IF ( ii.EQ.1 ) THEN
    149             IF(lstardis) THEN
    150                cdivu  = 1./zllm
    151             ELSE
    152                cdivu  = zllm **( -1./nitergdiv )
    153             ENDIF
    154          ELSE
    155             IF(lstardis) THEN
    156                crot   = 1./ zllm
    157             ELSE
    158                crot   = zllm **( -1./nitergrot )
    159             ENDIF
    160          ENDIF
    161 
    162  20   CONTINUE
    163 
    164 c   petit test pour les operateurs non star:
    165 c   ----------------------------------------
    166 
    167 c     IF(.NOT.lstardis) THEN
    168          fact    = rad*24./REAL(jjm)
    169          fact    = fact*fact
    170          PRINT*,'coef u ', fact/cdivu, 1./cdivu
    171          PRINT*,'coef r ', fact/crot , 1./crot
    172          PRINT*,'coef h ', fact/cdivh, 1./cdivh
    173 c     ENDIF
    174 
    175 c-----------------------------------------------------------------------
    176 c   variation verticale du coefficient de dissipation:
    177 c   --------------------------------------------------
    178 
    179 c First step: going from 1 to dissip_fac_mid (in gcm.def)
    180 c============
    181       DO l=1,llm
    182          zz      = 1. - preff/presnivs(l)
    183          zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
    184       ENDDO
    185 
    186          write(*,*) 'Dissipation : '
    187          write(*,*) 'Multiplication de la dissipation en altitude :'
    188          write(*,*) '  dissip_fac_mid =', dissip_fac_mid
    189 
    190 c Second step if ok_strato:  from dissip_fac_mid to dissip_fac_up (in gcm.def)
    191 c==========================
    192 c Utilisation de la fonction tangente hyperbolique pour augmenter
    193 c arbitrairement la dissipation et donc la stabilite du modele a
    194 c partir d'une certaine altitude.
    195 
    196 c   Le facteur multiplicatif de basse atmosphere etant deja pris
    197 c   en compte, il faut diviser le facteur multiplicatif de haute
    198 c   atmosphere par celui-ci.
    199 
    200       if (ok_strato) then
    201 
    202        Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
    203        do l=1,llm
    204          zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1)
    205      &   *(1-(0.5*(1+tanh(-6./dissip_deltaz*
    206      &    (-dissip_hdelta*log(presnivs(l)/Pup))  ))))  ))
    207        enddo
    208 
    209          write(*,*) '  dissip_fac_up =', dissip_fac_up
    210          write(*,*) 'Transition mid /up:  Pupstart,delta =',
    211      &             dissip_pupstart,'Pa', dissip_deltaz , '(km)'
    212 
    213       endif
    214 
    215 
    216       PRINT*,'Constantes de temps de la diffusion horizontale'
    217 
    218       tetamin =  1.e+6
    219 
    220       DO l=1,llm
    221         tetaudiv(l)   = zvert(l)/tetagdiv
    222         tetaurot(l)   = zvert(l)/tetagrot
    223         tetah(l)      = zvert(l)/tetatemp
    224 
    225         IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l)
    226         IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l)
    227         IF( tetamin.GT. (1./   tetah(l)) ) tetamin = 1./    tetah(l)
    228       ENDDO
    229 
    230       PRINT *,' INIDI tetamin dtvr ',tetamin,dtvr,iperiod
    231       idissip = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod
    232       PRINT *,' INIDI tetamin idissip ',tetamin,idissip
    233       idissip = MAX(iperiod,idissip)
    234       dtdiss  = idissip * dtvr
    235       PRINT *,' INIDI idissip dtdiss ',idissip,dtdiss
    236 
    237       DO l = 1,llm
    238          PRINT*,zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l),
    239      *                   dtdiss*tetah(l)
    240       ENDDO
    241 
    242 c
    243       RETURN
    244       END
     57     ENDDO
     58  ENDDO
     59
     60  idum  = -1
     61  zh(1) = RAN1(idum)-.5
     62  idum  = 0
     63  DO ij = 2, ip1jmp1
     64     zh(ij) = RAN1(idum) -.5
     65  ENDDO
     66
     67  CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)
     68
     69  CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
     70
     71  IF ( zhmin .GE. zhmax  )     THEN
     72     write(lunout,*)'  Inidissip  zh min max  ',zhmin,zhmax
     73     abort_message='probleme generateur alleatoire dans inidissip'
     74     call abort_gcm('inidissip',abort_message,1)
     75  ENDIF
     76
     77  zllm = ABS( zhmax )
     78  DO l = 1,50
     79     IF(lstardis) THEN
     80        CALL divgrad2(1,zh,deltap,niterh,zh)
     81     ELSE
     82        CALL divgrad (1,zh,niterh,zh)
     83     ENDIF
     84
     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
     92  ENDDO
     93
     94  IF(lstardis) THEN
     95     cdivh = 1./ zllm
     96  ELSE
     97     cdivh = zllm ** ( -1./niterh )
     98  ENDIF
     99
     100  !   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)
     101  !   -----------------------------------------------------------------
     102  write(lunout,*)'inidissip: calcul des valeurs propres'
     103
     104  DO    ii = 1, 2
     105     !
     106     DO ij = 1, ip1jmp1
     107        zu(ij)  = RAN1(idum) -.5
     108     ENDDO
     109     CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)
     110     DO ij = 1, ip1jm
     111        zv(ij) = RAN1(idum) -.5
     112     ENDDO
     113     CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)
     114
     115     CALL minmax(iip1*jjp1,zu,umin,ullm )
     116     CALL minmax(iip1*jjm, zv,vmin,vllm )
     117
     118     ullm = ABS ( ullm )
     119     vllm = ABS ( vllm )
     120
     121     DO    l = 1, 50
     122        IF(ii.EQ.1) THEN
     123           !cccc             CALL covcont( 1,zu,zv,zu,zv )
     124           IF(lstardis) THEN
     125              CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv )
     126           ELSE
     127              CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv )
     128           ENDIF
     129        ELSE
     130           IF(lstardis) THEN
     131              CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv )
     132           ELSE
     133              CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv )
     134           ENDIF
     135        ENDIF
     136
     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
     151     end DO
     152
     153     IF ( ii.EQ.1 ) THEN
     154        IF(lstardis) THEN
     155           cdivu  = 1./zllm
     156        ELSE
     157           cdivu  = zllm **( -1./nitergdiv )
     158        ENDIF
     159     ELSE
     160        IF(lstardis) THEN
     161           crot   = 1./ zllm
     162        ELSE
     163           crot   = zllm **( -1./nitergrot )
     164        ENDIF
     165     ENDIF
     166
     167  end DO
     168
     169  !   petit test pour les operateurs non star:
     170  !   ----------------------------------------
     171
     172  !     IF(.NOT.lstardis) THEN
     173  fact    = rad*24./REAL(jjm)
     174  fact    = fact*fact
     175  write(lunout,*)'inidissip: coef u ', fact/cdivu, 1./cdivu
     176  write(lunout,*)'inidissip: coef r ', fact/crot , 1./crot
     177  write(lunout,*)'inidissip: coef h ', fact/cdivh, 1./cdivh
     178  !     ENDIF
     179
     180  !-----------------------------------------------------------------------
     181  !   variation verticale du coefficient de dissipation:
     182  !   --------------------------------------------------
     183
     184! First step: going from 1 to dissip_fac_mid (in gcm.def)
     185!============
     186  DO l=1,llm
     187     zz      = 1. - preff/presnivs(l)
     188     zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
     189  ENDDO
     190
     191  write(lunout,*) 'Dissipation : '
     192  write(lunout,*) 'Multiplication de la dissipation en altitude :'
     193  write(lunout,*) '  dissip_fac_mid =', dissip_fac_mid
     194
     195! Second step if ok_strato:  from dissip_fac_mid to dissip_fac_up (in gcm.def)
     196!==========================
     197! Utilisation de la fonction tangente hyperbolique pour augmenter
     198! arbitrairement la dissipation et donc la stabilite du modele a
     199! partir d'une certaine altitude.
     200
     201!   Le facteur multiplicatif de basse atmosphere etant deja pris
     202!   en compte, il faut diviser le facteur multiplicatif de haute
     203!   atmosphere par celui-ci.
     204
     205  if (ok_strato) then
     206
     207    Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
     208    do l=1,llm
     209      zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1) &
     210                *(1-(0.5*(1+tanh(-6./dissip_deltaz*              &
     211               (-dissip_hdelta*log(presnivs(l)/Pup))  ))))  ))
     212    enddo
     213
     214    write(*,*) '  dissip_fac_up =', dissip_fac_up
     215    write(*,*) 'Transition mid /up:  Pupstart,delta =',           &
     216                   dissip_pupstart,'Pa', dissip_deltaz , '(km)'
     217
     218  endif
     219
     220
     221  write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale'
     222
     223  tetamin =  1.e+6
     224
     225  DO l=1,llm
     226     tetaudiv(l)   = zvert(l)/tetagdiv
     227     tetaurot(l)   = zvert(l)/tetagrot
     228     tetah(l)      = zvert(l)/tetatemp
     229
     230     IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l)
     231     IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l)
     232     IF( tetamin.GT. (1./   tetah(l)) ) tetamin = 1./    tetah(l)
     233  ENDDO
     234
     235  ! If dissip_period=0 calculate value for dissipation period, else keep value read from gcm.def
     236  IF (dissip_period == 0) THEN
     237     dissip_period = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod
     238     write(lunout,*)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ',tetamin,dtvr,iperiod,dissip_period
     239     dissip_period = MAX(iperiod,dissip_period)
     240  END IF
     241
     242  dtdiss  = dissip_period * dtvr
     243  write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr
     244
     245  DO l = 1,llm
     246     write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), &
     247          dtdiss*tetah(l)
     248  ENDDO
     249
     250END SUBROUTINE inidissip

  • trunk/LMDZ.COMMON/libf/dyn3dpar/integrd_p.F

    r7 r270  
    11!
    2 ! $Id: integrd_p.F 1446 2010-10-22 09:27:25Z emillour $
     2! $Id: integrd_p.F 1550 2011-07-05 09:44:55Z lguez $
    33!
    44      SUBROUTINE integrd_p
    … …  
    128128         PRINT*,' Au point ij = ',stop_it, ' , pression sol neg. '
    129129     &         , ps(stop_it)
    130          STOP' dans integrd'
     130         print *, ' dans integrd'
     131         stop 1
    131132        ENDIF
    132133

  • trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F

    r130 r270  
    168168
    169169      character*80 dynhist_file, dynhistave_file
    170       character*20 modname
     170      character(len=*),parameter :: modname="leapfrog"
    171171      character*80 abort_message
    172172
    … …  
    226226      endif
    227227      itaufinp1 = itaufin +1
    228       modname="leapfrog_p"
    229228
    230229      itau = 0
    … …  
    392391      ! Purely Matsuno time stepping
    393392         IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
    394          IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
     393         IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward )
     394     s        apdiss = .TRUE.
    395395         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
    396396     s          .and. iflag_phys.EQ.1                 ) apphys = .TRUE.
    … …  
    398398      ! Leapfrog/Matsuno time stepping
    399399         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
    400          IF( MOD(itau+1,idissip)  .EQ. 0              ) apdiss = .TRUE.
     400         IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
     401     s        apdiss = .TRUE.
    401402         IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
    402403      END IF
    … …  
    650651c   -------------------------------------------------------------
    651652
    652       IF( forward. OR . leapf )  THEN
     653!      IF( forward. OR . leapf )  THEN
     654      IF((.not.forward).OR. leapf )  THEN
     655        ! Ehouarn: gather mass fluxes during backward Matsuno or LF step
    653656cc$OMP PARALLEL DEFAULT(SHARED)
    654657c

  • trunk/LMDZ.COMMON/libf/dyn3dpar/limit_netcdf.F90

    r7 r270  
    11!
    2 ! $Id: limit_netcdf.F90 1441 2010-10-13 13:06:56Z emillour $
     2! $Id: limit_netcdf.F90 1508 2011-04-15 13:05:34Z jghattas $
    33!-------------------------------------------------------------------------------
    44!
    … …  
    4242  REAL, DIMENSION(iip1,jjp1), INTENT(IN) :: masque   ! land mask
    4343#ifndef CPP_EARTH
    44   WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
     44  CALL abort_gcm('limit_netcdf','Earth-specific routine, needs Earth physics',1)
    4545#else
    4646!-------------------------------------------------------------------------------
    … …  
    5252#include "indicesol.h"
    5353
    54 !--- For fractionary sub-cell use (old coding used soil index: 0,1,2,3) --------
    55   LOGICAL, PARAMETER :: fracterre=.TRUE.
    56 
    5754!--- INPUT NETCDF FILES NAMES --------------------------------------------------
    5855  CHARACTER(LEN=25) :: icefile, sstfile, dumstr
    5956  CHARACTER(LEN=25), PARAMETER :: famipsst='amipbc_sst_1x1.nc        ',        &
    6057                                  famipsic='amipbc_sic_1x1.nc        ',        &
    61                                   fclimsst='amipbc_sst_1x1_clim.nc   ',        &
    62                                   fclimsic='amipbc_sic_1x1_clim.nc   ',        &
    6358                                  fcpldsst='cpl_atm_sst.nc           ',        &
    6459                                  fcpldsic='cpl_atm_sic.nc           ',        &
     60                                  fhistsst='histmth_sst.nc           ',        &
     61                                  fhistsic='histmth_sic.nc           ',        &
    6562                                  frugo   ='Rugos.nc                 ',        &
    6663                                  falbe   ='Albedo.nc                '
    67 
     64  CHARACTER(LEN=10) :: varname
    6865!--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------
    6966  REAL,   DIMENSION(klon)                :: fi_ice, verif
    … …  
    8077  INTEGER :: id_FOCE, id_FSIC, id_FTER, id_FLIC
    8178  INTEGER :: NF90_FORMAT
    82   LOGICAL :: lCPL                    !--- T: IPCC-IPSL cpl model output files
    8379  INTEGER :: ndays                   !--- Depending on the output calendar
    8480
    … …  
    104100
    105101!--- RUGOSITY TREATMENT --------------------------------------------------------
    106   WRITE(lunout,*) 'Traitement de la rugosite'
    107   CALL get_2Dfield(frugo,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:))
     102  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la rugosite'
     103  varname='RUGOS'
     104  CALL get_2Dfield(frugo,varname,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:))
    108105
    109106!--- OCEAN TREATMENT -----------------------------------------------------------
    110   PRINT*, 'Traitement de la glace oceanique' ; icefile=''; lCPL=.FALSE.
     107  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la glace oceanique'
    111108
    112109! Input SIC file selection
    113   icefile='fake'
    114   IF(NF90_OPEN(famipsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(famipsic)
    115   IF(NF90_OPEN(fclimsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fclimsic)
    116   IF(NF90_OPEN(fcpldsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fcpldsic)
    117   SELECT CASE(icefile)
    118     CASE(famipsic); dumstr='Amip.'
    119     CASE(fclimsic); dumstr='Amip climatologique.'
    120     CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE.
    121     CASE('fake');   CALL abort_gcm('limit_netcdf','Fichier SIC non reconnu.',1)
    122   END SELECT
     110! Open file only to test if available
     111  IF ( NF90_OPEN(TRIM(famipsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
     112     icefile=TRIM(famipsic)
     113     varname='sicbcs'
     114  ELSE IF( NF90_OPEN(TRIM(fcpldsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
     115     icefile=TRIM(fcpldsic)
     116     varname='SIICECOV'
     117  ELSE IF ( NF90_OPEN(TRIM(fhistsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
     118     icefile=TRIM(fhistsic)
     119     varname='pourc_sic'
     120  ELSE
     121     WRITE(lunout,*) 'ERROR! No sea-ice input file was found.'
     122     WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsic),', ',trim(fcpldsic),', ',trim(fhistsic)
     123     CALL abort_gcm('limit_netcdf','No sea-ice file was found',1)
     124  END IF
    123125  ierr=NF90_CLOSE(nid)
    124   WRITE(lunout,*)'Pour la glace de mer a ete choisi un fichier '//TRIM(dumstr)
    125 
    126   CALL get_2Dfield(icefile,'SIC',interbar,ndays,phy_ice,flag=oldice,lCPL=lCPL)
     126  IF (prt_level>=0) WRITE(lunout,*)'Pour la glace de mer a ete choisi le fichier '//TRIM(icefile)
     127
     128  CALL get_2Dfield(icefile,varname, 'SIC',interbar,ndays,phy_ice,flag=oldice)
    127129
    128130  ALLOCATE(pctsrf_t(klon,nbsrf,ndays))
    129131  DO k=1,ndays
    130     fi_ice=phy_ice(:,k)
    131     WHERE(fi_ice>=1.0  ) fi_ice=1.0
    132     WHERE(fi_ice<EPSFRA) fi_ice=0.0
    133     IF(fracterre) THEN
    134       pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)       ! land soil
    135       pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic)       ! land ice
    136       IF(lCPL) THEN                               ! SIC=pICE*(1-LIC-TER)
    137         pctsrf_t(:,is_sic,k)=fi_ice*(1-pctsrf(:,is_lic)-pctsrf(:,is_ter))
    138       ELSE                                        ! SIC=pICE-LIC
     132     fi_ice=phy_ice(:,k)
     133     WHERE(fi_ice>=1.0  ) fi_ice=1.0
     134     WHERE(fi_ice<EPSFRA) fi_ice=0.0
     135     pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)       ! land soil
     136     pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic)       ! land ice
     137     IF (icefile==trim(fcpldsic)) THEN           ! SIC=pICE*(1-LIC-TER)
     138        pctsrf_t(:,is_sic,k)=fi_ice(:)*(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter))
     139     ELSE IF (icefile==trim(fhistsic)) THEN      ! SIC=pICE
     140        pctsrf_t(:,is_sic,k)=fi_ice(:)
     141     ELSE ! icefile==famipsic                    ! SIC=pICE-LIC
    139142        pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k)
    140       END IF
    141       WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
    142       WHERE(1.0-zmasq<EPSFRA)
     143     END IF
     144     WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
     145     WHERE(1.0-zmasq<EPSFRA)
    143146        pctsrf_t(:,is_sic,k)=0.0
    144147        pctsrf_t(:,is_oce,k)=0.0
    145       ELSEWHERE
     148     ELSEWHERE
    146149        WHERE(pctsrf_t(:,is_sic,k)>=1.0-zmasq)
    147           pctsrf_t(:,is_sic,k)=1.0-zmasq
    148           pctsrf_t(:,is_oce,k)=0.0
     150           pctsrf_t(:,is_sic,k)=1.0-zmasq
     151           pctsrf_t(:,is_oce,k)=0.0
    149152        ELSEWHERE
    150           pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)
    151           WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)
    152             pctsrf_t(:,is_oce,k)=0.0
    153             pctsrf_t(:,is_sic,k)=1.0-zmasq
    154           END WHERE
     153           pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)
     154           WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)
     155              pctsrf_t(:,is_oce,k)=0.0
     156              pctsrf_t(:,is_sic,k)=1.0-zmasq
     157           END WHERE
    155158        END WHERE
    156       END WHERE
    157       nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0)
    158       IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
    159       nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA)
    160       IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad
    161     ELSE
    162       pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)
    163       WHERE(NINT(pctsrf(:,is_ter))==1)
    164         pctsrf_t(:,is_sic,k)=0.
    165         pctsrf_t(:,is_oce,k)=0.                 
    166         WHERE(fi_ice>=1.e-5)
    167           pctsrf_t(:,is_lic,k)=fi_ice
    168           pctsrf_t(:,is_ter,k)=pctsrf_t(:,is_ter,k)-pctsrf_t(:,is_lic,k)
    169         ELSEWHERE
    170           pctsrf_t(:,is_lic,k)=0.0
    171         END WHERE
    172       ELSEWHERE
    173         pctsrf_t(:,is_lic,k) = 0.0
    174         WHERE(fi_ice>=1.e-5)
    175           pctsrf_t(:,is_ter,k)=0.0
    176           pctsrf_t(:,is_sic,k)=fi_ice
    177           pctsrf_t(:,is_oce,k)=1.0-pctsrf_t(:,is_sic,k)
    178         ELSEWHERE
    179           pctsrf_t(:,is_sic,k)=0.0
    180           pctsrf_t(:,is_oce,k)=1.0
    181         END WHERE
    182       END WHERE
    183       verif=sum(pctsrf_t(:,:,k),dim=2)
    184       nbad=COUNT(verif<1.0-1.e-5.OR.verif>1.0+1.e-5)
    185       IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
    186     END IF
     159     END WHERE
     160     nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0)
     161     IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
     162     nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA)
     163     IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad
    187164  END DO
    188165  DEALLOCATE(phy_ice)
    189166
    190167!--- SST TREATMENT -------------------------------------------------------------
    191   WRITE(lunout,*) 'Traitement de la sst' ; sstfile=''; lCPL=.FALSE.
     168  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la sst'
    192169
    193170! Input SST file selection
    194   sstfile='fake'
    195   IF(NF90_OPEN(famipsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(famipsst)
    196   IF(NF90_OPEN(fclimsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fclimsst)
    197   IF(NF90_OPEN(fcpldsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fcpldsst)
    198   SELECT CASE(icefile)
    199     CASE(famipsic); dumstr='Amip.'
    200     CASE(fclimsic); dumstr='Amip climatologique.'
    201     CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE.
    202     CASE('fake');   CALL abort_gcm('limit_netcdf','Fichier SST non reconnu',1)
    203   END SELECT
     171! Open file only to test if available
     172  IF ( NF90_OPEN(TRIM(famipsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
     173     sstfile=TRIM(famipsst)
     174     varname='tosbcs'
     175  ELSE IF ( NF90_OPEN(TRIM(fcpldsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
     176     sstfile=TRIM(fcpldsst)
     177     varname='SISUTESW'
     178  ELSE IF ( NF90_OPEN(TRIM(fhistsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
     179     sstfile=TRIM(fhistsst)
     180     varname='tsol_oce'
     181  ELSE
     182     WRITE(lunout,*) 'ERROR! No sst input file was found.'
     183     WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsst),trim(fcpldsst),trim(fhistsst)
     184     CALL abort_gcm('limit_netcdf','No sst file was found',1)
     185  END IF
    204186  ierr=NF90_CLOSE(nid)
    205   WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(dumstr)
    206 
    207   CALL get_2Dfield(trim(sstfile),'SST',interbar,ndays,phy_sst,flag=extrap)
     187  IF (prt_level>=0) WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(sstfile)
     188
     189  CALL get_2Dfield(sstfile,varname,'SST',interbar,ndays,phy_sst,flag=extrap)
    208190
    209191!--- ALBEDO TREATMENT ----------------------------------------------------------
    210   WRITE(lunout,*) 'Traitement de l albedo'
    211   CALL get_2Dfield(falbe,'ALB',interbar,ndays,phy_alb)
     192  IF (prt_level>1) WRITE(lunout,*) 'Traitement de l albedo'
     193  varname='ALBEDO'
     194  CALL get_2Dfield(falbe,varname,'ALB',interbar,ndays,phy_alb)
    212195
    213196!--- REFERENCE GROUND HEAT FLUX TREATMENT --------------------------------------
    … …  
    215198
    216199!--- OUTPUT FILE WRITING -------------------------------------------------------
    217   WRITE(lunout,*) 'Ecriture du fichier limit'
     200  IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : debut'
    218201
    219202  !--- File creation
    … …  
    264247  ierr=NF90_CLOSE(nid)
    265248
     249  IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : fin'
     250
    266251  DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug)
    267252
    … …  
    276261!-------------------------------------------------------------------------------
    277262!
    278 SUBROUTINE get_2Dfield(fnam, mode, ibar, ndays, champo, flag, mask, lCPL)
     263SUBROUTINE get_2Dfield(fnam, varname, mode, ibar, ndays, champo, flag, mask)
    279264!
    280265!-----------------------------------------------------------------------------
    … …  
    304289! Arguments:
    305290  CHARACTER(LEN=*),  INTENT(IN)     :: fnam     ! NetCDF file name
     291  CHARACTER(LEN=10), INTENT(IN)     :: varname  ! NetCDF variable name
    306292  CHARACTER(LEN=3),  INTENT(IN)     :: mode     ! RUG, SIC, SST or ALB
    307293  LOGICAL,           INTENT(IN)     :: ibar     ! interp on pressure levels
    308294  INTEGER,           INTENT(IN)     :: ndays    ! current year number of days
    309   REAL,    POINTER,  DIMENSION(:, :) :: champo   ! output field = f(t)
     295  REAL,    POINTER,  DIMENSION(:, :) :: champo  ! output field = f(t)
    310296  LOGICAL, OPTIONAL, INTENT(IN)     :: flag     ! extrapol. (SST) old ice (SIC)
    311297  REAL,    OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask
    312   LOGICAL, OPTIONAL, INTENT(IN)     :: lCPL     ! Coupled model flag (for ICE)
    313298!------------------------------------------------------------------------------
    314299! Local variables:
    … …  
    316301  INTEGER :: ncid, varid                  ! NetCDF identifiers
    317302  CHARACTER(LEN=30)               :: dnam       ! dimension name
    318   CHARACTER(LEN=80)               :: varname    ! NetCDF variable name
    319303!--- dimensions
    320304  INTEGER,           DIMENSION(4) :: dids       ! NetCDF dimensions identifiers
    … …  
    331315!--- input files
    332316  CHARACTER(LEN=20)                 :: cal_in   ! calendar
     317  CHARACTER(LEN=20)                 :: unit_sic ! attribute unit in sea-ice file
    333318  INTEGER                           :: ndays_in ! number of days
    334319!--- misc
    … …  
    337322  CHARACTER(LEN=25)                 :: title    ! for messages
    338323  LOGICAL                           :: extrp    ! flag for extrapolation
     324  LOGICAL                           :: oldice   ! flag for old way ice computation
    339325  REAL                              :: chmin, chmax
    340326  INTEGER ierr
    341327  integer n_extrap ! number of extrapolated points
    342328  logical skip
     329
    343330!------------------------------------------------------------------------------
    344331!---Variables depending on keyword 'mode' -------------------------------------
    345332  NULLIFY(champo)
     333
    346334  SELECT CASE(mode)
    347     CASE('RUG'); varname='RUGOS'; title='Rugosite'
    348     CASE('SIC'); varname='sicbcs'; title='Sea-ice'
    349     CASE('SST'); varname='tosbcs'; title='SST'
    350     CASE('ALB'); varname='ALBEDO'; title='Albedo'
     335  CASE('RUG'); title='Rugosite'
     336  CASE('SIC'); title='Sea-ice'
     337  CASE('SST'); title='SST'
     338  CASE('ALB'); title='Albedo'
    351339  END SELECT
     340 
     341
    352342  extrp=.FALSE.
     343  oldice=.FALSE.
    353344  IF ( PRESENT(flag) ) THEN
    354345    IF ( flag .AND. mode=='SST' ) extrp=.TRUE.
     346    IF ( flag .AND. mode=='SIC' ) oldice=.TRUE.
    355347  END IF
    356348
    357349!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
     350  IF (prt_level>5) WRITE(lunout,*) ' Now reading file : ',fnam
    358351  ierr=NF90_OPEN(fnam, NF90_NOWRITE, ncid);             CALL ncerr(ierr, fnam)
    359   ierr=NF90_INQ_VARID(ncid, varname, varid);            CALL ncerr(ierr, fnam)
     352  ierr=NF90_INQ_VARID(ncid, trim(varname), varid);            CALL ncerr(ierr, fnam)
    360353  ierr=NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids); CALL ncerr(ierr, fnam)
     354
     355!--- Read unit for sea-ice variable only
     356  IF (mode=='SIC') THEN
     357     ierr=NF90_GET_ATT(ncid, varid, 'units', unit_sic)
     358     IF(ierr/=NF90_NOERR) THEN
     359        IF (prt_level>5) WRITE(lunout,*) 'No unit was given in sea-ice file. Take percentage as default value'
     360        unit_sic='X'
     361     ELSE
     362        IF (prt_level>5) WRITE(lunout,*) ' Sea-ice cover has unit=',unit_sic
     363     END IF
     364  END IF
    361365
    362366!--- Longitude
    … …  
    365369  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
    366370  ierr=NF90_GET_VAR(ncid, varid, dlon_ini);              CALL ncerr(ierr, fnam)
    367   WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep
     371  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep
    368372
    369373!--- Latitude
    … …  
    372376  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
    373377  ierr=NF90_GET_VAR(ncid, varid, dlat_ini);              CALL ncerr(ierr, fnam)
    374   WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep
     378  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep
    375379
    376380!--- Time (variable is not needed - it is rebuilt - but calendar is)
    … …  
    385389      CASE('SIC', 'SST'); cal_in='gregorian'
    386390    END SELECT
    387     WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' &
     391    IF (prt_level>5) WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' &
    388392         // 'dans '//TRIM(fnam)//'. On choisit la valeur par defaut.'
    389393  END IF
    390   WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', &
     394  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', &
    391395       cal_in
    392396
     397 
    393398!--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION --------------------
    394399  !--- Determining input file number of days, depending on calendar
    … …  
    398403!--- If input records are not monthly, time sampling has to be constant !
    399404  timeyear=mid_months(anneeref, cal_in, lmdep)
    400   IF (lmdep /= 12) WRITE(lunout, '(a, i3, a)') 'Note : les fichiers de ' &
    401        // TRIM(mode) // ' ne comportent pas 12, mais ', lmdep, &
    402        ' enregistrements.'
     405  IF (lmdep /= 12) WRITE(lunout,*) 'Note : les fichiers de ', TRIM(mode), &
     406       ' ne comportent pas 12, mais ', lmdep, ' enregistrements.'
    403407
    404408!--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
    … …  
    406410  IF(extrp) ALLOCATE(work(imdep, jmdep))
    407411
    408   WRITE(lunout, *)
    409   WRITE(lunout, '(a, i3, a)')'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, &
    410        ' CHAMPS.'
     412  IF (prt_level>5) WRITE(lunout, *)
     413  IF (prt_level>5) WRITE(lunout,*)'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, ' CHAMPS.'
    411414  ierr=NF90_INQ_VARID(ncid, varname, varid);             CALL ncerr(ierr, fnam)
    412415  DO l=1, lmdep
    … …  
    419422         work)
    420423
    421     IF(ibar.AND..NOT.(mode=='SIC'.AND.flag)) THEN
    422       IF(l==1) THEN
    423         WRITE(lunout, *)                                                      &
    424   '-------------------------------------------------------------------------'
    425         WRITE(lunout, *)                                                     &
    426   'Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$'
    427         WRITE(lunout, *)                                                      &
    428   '-------------------------------------------------------------------------'
     424    IF(ibar .AND. .NOT.oldice) THEN
     425      IF(l==1 .AND. prt_level>5) THEN
     426        WRITE(lunout, *) '-------------------------------------------------------------------------'
     427        WRITE(lunout, *) 'Utilisation de l''interpolation barycentrique pour ',TRIM(title),' $$$'
     428        WRITE(lunout, *) '-------------------------------------------------------------------------'
    429429      END IF
    430430      IF(mode=='RUG') champ=LOG(champ)
    … …  
    453453
    454454!--- TIME INTERPOLATION ------------------------------------------------------
    455   WRITE(lunout, *)
    456   WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
    457   WRITE(lunout, "(2x, ' Vecteur temps en entree: ', 10f6.1)") timeyear
    458   WRITE(lunout, "(2x, ' Vecteur temps en sortie de 0 a ', i3)") ndays
     455  IF (prt_level>5) THEN
     456     WRITE(lunout, *)
     457     WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
     458     WRITE(lunout, *)' Vecteur temps en entree: ', timeyear
     459     WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays
     460  END IF
     461
    459462  ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays))
    460463  skip = .false.
    … …  
    471474  END DO
    472475  if (n_extrap /= 0) then
    473      print *, "get_2Dfield pchfe_95: n_extrap = ", n_extrap
     476     WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
    474477  end if
    475478  champan(iip1, :, :)=champan(1, :, :)
    … …  
    479482  DO j=1, jjp1
    480483    CALL minmax(iip1, champan(1, j, 10), chmin, chmax)
    481     WRITE(lunout, *)' '//TRIM(title)//' au temps 10 ', chmin, chmax, j
     484    IF (prt_level>5) WRITE(lunout, *)' ',TRIM(title),' au temps 10 ', chmin, chmax, j
    482485  END DO
    483486
    484487!--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
    485488  IF(mode=='SST') THEN
    486     WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38'
     489    IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38'
    487490    WHERE(champan<271.38) champan=271.38
    488491  END IF
    … …  
    490493!--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
    491494  IF(mode=='SIC') THEN
    492     WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
    493     IF(.NOT.lCPL) champan(:, :, :)=champan(:, :, :)/100.
     495    IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
     496
     497    IF (unit_sic=='1') THEN
     498       ! Nothing to be done for sea-ice field is already in fraction of 1
     499       ! This is the case for sea-ice in file cpl_atm_sic.nc
     500       IF (prt_level>5) WRITE(lunout,*) 'Sea-ice field already in fraction of 1'
     501    ELSE
     502       ! Convert sea ice from percentage to fraction of 1
     503       IF (prt_level>5) WRITE(lunout,*) 'Transformt sea-ice field from percentage to fraction of 1.'
     504       champan(:, :, :)=champan(:, :, :)/100.
     505    END IF
     506
    494507    champan(iip1, :, :)=champan(1, :, :)
    495508    WHERE(champan>1.0) champan=1.0
    496509    WHERE(champan<0.0) champan=0.0
    497   END IF
     510 END IF
    498511
    499512!--- DYNAMICAL TO PHYSICAL GRID ----------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.