Changeset 66 for trunk/libf


Ignore:
Timestamp:
Feb 16, 2011, 4:57:45 PM (14 years ago)
Author:
emillour
Message:

EM: Mise a niveau par rapport a la version terrestre (LMDZ5V2.0-dev, rev 1487)

  • Mise a jour des scripts (terrestres) 'makegcm' et 'create_make_gcm'
  • Ajout du script 'makelmdz' (version "amelioree", en Bash, de makegcm)
  • Mise a jour des routines dans phylmd (sauf regr_lat_time_climoz_m.F)
  • disvert (dans dyn3d et dyn3dpar): passage au Fortran 90
  • parallel.F90 (dyn3dpar): correction bug
  • etat0_netcdf.F90 (dyn3d et dyn3dpar) : mise a jour mineure
  • ce0l.F90 (dyn3dpar) : correction bug
  • abort_gcm.F (dyn3dpar) : correction bug
  • ugeostr.F90 (dyn3d et dyn3dpar) : passage au Fortran 90
  • fluxstokenc_p.F (dyn3dpar) : correction bug
  • iniacademic.F90 (dyn3d et dyn3dpar) : passage au Fortran 90
  • friction_p.F (dyn3dpar) : correction bug
  • infotrac.F90 (dyn3d et dyn3dpar) : correction bug mineur sur lecture traceurs
  • caladvtrac.F (dyn3d) : modifications cosmetiques
Location:
trunk/libf
Files:
7 added
3 deleted
29 edited
6 moved

Legend:

Unmodified
Added
Removed

  • trunk/libf/dyn3d/caladvtrac.F

    r7 r66  
    88     *                   flxw, pk)
    99c
    10       USE infotrac
     10      USE infotrac, ONLY : nqtot
    1111      USE control_mod, ONLY : iapp_tracvl,planet_type
    1212 

  • trunk/libf/dyn3d/disvert.F90

    r64 r66  
    1 !
    2 ! $Id: disvert.F 1403 2010-07-01 09:02:53Z fairhead $
    3 !
    4       SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
     1! $Id: disvert.F90 1480 2011-01-31 21:29:58Z jghattas $
    52
    6 c    Auteur :  P. Le Van .
    7 c
    8       IMPLICIT NONE
     3SUBROUTINE disvert(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
    94
    10 #include "dimensions.h"
    11 #include "paramet.h"
    12 #include "iniprint.h"
    13 #include "logic.h"
    14 c
    15 c=======================================================================
    16 c
    17 c
    18 c    s = sigma ** kappa   :  coordonnee  verticale
    19 c    dsig(l)            : epaisseur de la couche l ds la coord.  s
    20 c    sig(l)             : sigma a l'interface des couches l et l-1
    21 c    ds(l)              : distance entre les couches l et l-1 en coord.s
    22 c
    23 c=======================================================================
    24 c
    25       REAL pa,preff
    26       REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
    27       REAL presnivs(llm)
    28 c
    29 c   declarations:
    30 c   -------------
    31 c
    32       REAL sig(llm+1),dsig(llm)
    33        real zzz(1:llm+1)
    34        real dzz(1:llm)
    35       real zk,zkm1,dzk1,dzk2,k0,k1
    36 c
    37       INTEGER l
    38       REAL snorm,dsigmin
    39       REAL alpha,beta,gama,delta,deltaz,h
    40       INTEGER np,ierr
    41       REAL pi,x
     5  ! Auteur : P. Le Van
    426
    43       REAL SSUM
    44 c
    45 c-----------------------------------------------------------------------
    46 c
    47       pi=2.*ASIN(1.)
     7  IMPLICIT NONE
    488
    49       OPEN(99,file='sigma.def',status='old',form='formatted',
    50      s   iostat=ierr)
     9  include "dimensions.h"
     10  include "paramet.h"
     11  include "iniprint.h"
     12  include "logic.h"
    5113
    52 c-----------------------------------------------------------------------
    53 c   cas 1 on lit les options dans sigma.def:
    54 c   ----------------------------------------
     14  ! s = sigma ** kappa : coordonnee verticale
     15  ! dsig(l) : epaisseur de la couche l ds la coord. s
     16  ! sig(l) : sigma a l'interface des couches l et l-1
     17  ! ds(l) : distance entre les couches l et l-1 en coord.s
    5518
    56       IF (ierr.eq.0) THEN
     19  REAL pa, preff
     20  REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1)
     21  REAL presnivs(llm)
    5722
    58       READ(99,*) h           ! hauteur d'echelle 8.
    59       READ(99,*) deltaz      ! epaiseur de la premiere couche 0.04
    60       READ(99,*) beta        ! facteur d'acroissement en haut 1.3
    61       READ(99,*) k0          ! nombre de couches dans la transition surf
    62       READ(99,*) k1          ! nombre de couches dans la transition haute
    63       CLOSE(99)
    64       alpha=deltaz/(llm*h)
    65       write(lunout,*)'h,alpha,k0,k1,beta'
     23  REAL sig(llm+1), dsig(llm)
     24  real zk, zkm1, dzk1, dzk2, k0, k1
    6625
    67 c     read(*,*) h,deltaz,beta,k0,k1 ! 8 0.04 4 20 1.2
     26  INTEGER l
     27  REAL dsigmin
     28  REAL alpha, beta, deltaz, h
     29  INTEGER iostat
     30  REAL pi, x
    6831
    69       alpha=deltaz/tanh(1./k0)*2.
    70       zkm1=0.
    71       sig(1)=1.
    72       do l=1,llm
    73         sig(l+1)=(cosh(l/k0))**(-alpha*k0/h)
    74      + *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
     32  !-----------------------------------------------------------------------
     33
     34  pi = 2 * ASIN(1.)
     35
     36  OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
     37
     38  IF (iostat == 0) THEN
     39     ! cas 1 on lit les options dans sigma.def:
     40     READ(99, *) h ! hauteur d'echelle 8.
     41     READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
     42     READ(99, *) beta ! facteur d'acroissement en haut 1.3
     43     READ(99, *) k0 ! nombre de couches dans la transition surf
     44     READ(99, *) k1 ! nombre de couches dans la transition haute
     45     CLOSE(99)
     46     alpha=deltaz/(llm*h)
     47     write(lunout, *)'disvert: h, alpha, k0, k1, beta'
     48
     49     alpha=deltaz/tanh(1./k0)*2.
     50     zkm1=0.
     51     sig(1)=1.
     52     do l=1, llm
     53        sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) &
     54             *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
    7555        zk=-h*log(sig(l+1))
    7656
    7757        dzk1=alpha*tanh(l/k0)
    7858        dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
    79         write(lunout,*)l,sig(l+1),zk,zk-zkm1,dzk1,dzk2
     59        write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2
    8060        zkm1=zk
    81       enddo
     61     enddo
    8262
    83       sig(llm+1)=0.
     63     sig(llm+1)=0.
    8464
    85 c
    86        DO 2  l = 1, llm
    87        dsig(l) = sig(l)-sig(l+1)
    88    2   CONTINUE
    89 c
     65     DO l = 1, llm
     66        dsig(l) = sig(l)-sig(l+1)
     67     end DO
     68  ELSE
     69     if (ok_strato) then
     70        if (llm==39) then
     71           dsigmin=0.3
     72        else if (llm==50) then
     73           dsigmin=1.
     74        else
     75           WRITE(LUNOUT,*)'disvert: ATTENTION discretisation z a ajuster'
     76           dsigmin=1.
     77        endif
     78        WRITE(LUNOUT,*) 'disvert: Discretisation verticale DSIGMIN=',dsigmin
     79     endif
    9080
    91       ELSE
    92 c-----------------------------------------------------------------------
    93 c   cas 2 ancienne discretisation (LMD5...):
    94 c   ----------------------------------------
     81     DO l = 1, llm
     82        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
    9583
    96       WRITE(LUNOUT,*)'WARNING!!! Ancienne discretisation verticale'
     84        IF (ok_strato) THEN
     85           dsig(l) =(dsigmin + 7. * SIN(x)**2) &
     86                *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
     87        ELSE
     88           dsig(l) = 1.0 + 7.0 * SIN(x)**2
     89        ENDIF
     90     ENDDO
     91     dsig = dsig / sum(dsig)
     92     sig(llm+1) = 0.
     93     DO l = llm, 1, -1
     94        sig(l) = sig(l+1) + dsig(l)
     95     ENDDO
     96  ENDIF
    9797
    98       if (ok_strato) then
    99          if (llm==39) then
    100             dsigmin=0.3
    101          else if (llm==50) then
    102             dsigmin=1.
    103          else
    104             WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster'
    105             dsigmin=1.
    106          endif
    107          WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin
    108       endif
     98  DO l=1, llm
     99     nivsigs(l) = REAL(l)
     100  ENDDO
    109101
    110       h=7.
    111       snorm  = 0.
    112       DO l = 1, llm
    113          x = 2.*asin(1.) * (REAL(l)-0.5) / REAL(llm+1)
     102  DO l=1, llmp1
     103     nivsig(l)= REAL(l)
     104  ENDDO
    114105
    115          IF (ok_strato) THEN
    116            dsig(l) =(dsigmin + 7.0 * SIN(x)**2)
    117      &            *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2       
    118          ELSE
    119            dsig(l) = 1.0 + 7.0 * SIN(x)**2
    120          ENDIF
     106  ! .... Calculs de ap(l) et de bp(l) ....
     107  ! ..... pa et preff sont lus sur les fichiers start par lectba .....
    121108
    122          snorm = snorm + dsig(l)
    123       ENDDO
    124       snorm = 1./snorm
    125       DO l = 1, llm
    126          dsig(l) = dsig(l)*snorm
    127       ENDDO
    128       sig(llm+1) = 0.
    129       DO l = llm, 1, -1
    130          sig(l) = sig(l+1) + dsig(l)
    131       ENDDO
     109  bp(llmp1) = 0.
    132110
    133       ENDIF
     111  DO l = 1, llm
     112     bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
     113     ap(l) = pa * ( sig(l) - bp(l) )
     114  ENDDO
    134115
     116  bp(1)=1.
     117  ap(1)=0.
    135118
    136       DO l=1,llm
    137         nivsigs(l) = REAL(l)
    138       ENDDO
     119  ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
    139120
    140       DO l=1,llmp1
    141         nivsig(l)= REAL(l)
    142       ENDDO
     121  write(lunout, *)'disvert: BP '
     122  write(lunout, *) bp
     123  write(lunout, *)'disvert: AP '
     124  write(lunout, *) ap
    143125
    144 c
    145 c    ....  Calculs  de ap(l) et de bp(l)  ....
    146 c    .........................................
    147 c
    148 c
    149 c   .....  pa et preff sont lus  sur les fichiers start par lectba  .....
    150 c
     126  write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
     127  write(lunout, *)'couches calcules pour une pression de surface =', preff
     128  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de'
     129  write(lunout, *)'8km'
     130  DO l = 1, llm
     131     dpres(l) = bp(l) - bp(l+1)
     132     presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
     133     write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
     134          log(preff/presnivs(l))*8. &
     135          , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ &
     136          max(ap(l+1)+bp(l+1)*preff, 1.e-10))
     137  ENDDO
    151138
    152       bp(llmp1) =   0.
     139  write(lunout, *) 'disvert: PRESNIVS '
     140  write(lunout, *) presnivs
    153141
    154       DO l = 1, llm
    155 cc
    156 ccc    ap(l) = 0.
    157 ccc    bp(l) = sig(l)
    158 
    159       bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
    160       ap(l) = pa * ( sig(l) - bp(l) )
    161 c
    162       ENDDO
    163 
    164       bp(1)=1.
    165       ap(1)=0.
    166 
    167       ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
    168 
    169       write(lunout,*)' BP '
    170       write(lunout,*)  bp
    171       write(lunout,*)' AP '
    172       write(lunout,*)  ap
    173 
    174       write(lunout,*)
    175      .'Niveaux de pressions approximatifs aux centres des'
    176       write(lunout,*)'couches calcules pour une pression de surface =',
    177      .                 preff
    178       write(lunout,*)
    179      .     'et altitudes equivalentes pour une hauteur d echelle de'
    180       write(lunout,*)'8km'
    181       DO l = 1, llm
    182        dpres(l) = bp(l) - bp(l+1)
    183        presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
    184        write(lunout,*)'PRESNIVS(',l,')=',presnivs(l),'    Z ~ ',
    185      .        log(preff/presnivs(l))*8.
    186      .  ,'   DZ ~ ',8.*log((ap(l)+bp(l)*preff)/
    187      .       max(ap(l+1)+bp(l+1)*preff,1.e-10))
    188       ENDDO
    189 
    190       write(lunout,*)' PRESNIVS '
    191       write(lunout,*)presnivs
    192 
    193       RETURN
    194       END
     142END SUBROUTINE disvert

  • trunk/libf/dyn3d/etat0_netcdf.F90

    r1 r66  
    11!
    2 ! $Id: etat0_netcdf.F90 1425 2010-09-02 13:45:23Z lguez $
     2! $Id: etat0_netcdf.F90 1486 2011-02-11 12:07:39Z fairhead $
    33!
    44!-------------------------------------------------------------------------------
    … …  
    9898  REAL    :: dummy
    9999  LOGICAL :: ok_newmicro, ok_journe, ok_mensuel, ok_instan, ok_hf
    100   LOGICAL :: ok_LES, ok_ade, ok_aie, aerosol_couple, new_aod
     100  LOGICAL :: ok_LES, ok_ade, ok_aie, aerosol_couple, new_aod, callstats
    101101  INTEGER :: iflag_radia, flag_aerosol
    102102  REAL    :: bl95_b0, bl95_b1, fact_cldcon, facttemps, ratqsbas, ratqshaut
    … …  
    130130!--- CONSTRUCT A GRID
    131131  CALL conf_phys(  ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES,     &
     132                   callstats,                                           &
    132133                   solarlong0,seuil_inversion,                          &
    133134                   fact_cldcon, facttemps,ok_newmicro,iflag_radia,      &

  • trunk/libf/dyn3d/infotrac.F90

    r37 r66  
    177177          ! Continue to read tracer.def
    178178          DO iq=1,nqtrue
    179              READ(90,999) hadv(iq),vadv(iq),tnom_0(iq)
     179             READ(90,*) hadv(iq),vadv(iq),tnom_0(iq)
    180180          END DO
    181181          CLOSE(90) 
    … …  
    234234          ! Continue to read tracer.def
    235235          DO iq=1,nqtrue
    236              READ(90,999) hadv(iq),vadv(iq),tnom_0(iq)
     236             READ(90,*) hadv(iq),vadv(iq),tnom_0(iq)
    237237          END DO
    238238          CLOSE(90) 
    … …  
    316316       tname(new_iq)= tnom_0(iq)
    317317       IF (iadv(new_iq)==0) THEN
    318           ttext(new_iq)=str1(1:len_trim(str1))
     318          ttext(new_iq)=trim(str1)
    319319       ELSE
    320           ttext(new_iq)=str1(1:len_trim(str1))//descrq(iadv(new_iq))
     320          ttext(new_iq)=trim(tnom_0(iq))//descrq(iadv(new_iq))
    321321       END IF
    322322
    … …  
    381381    DEALLOCATE(tracnam)
    382382
    383 999 FORMAT (i2,1x,i2,1x,a15)
    384 
    385383  END SUBROUTINE infotrac_init
    386384

  • trunk/libf/dyn3d/iniacademic.F90

    r64 r66  
    11!
    2 ! $Id: iniacademic.F 1446 2010-10-22 09:27:25Z emillour $
     2! $Id: iniacademic.F90 1474 2011-01-14 11:04:45Z lguez $
    33!
    4 c
    5 c
    6       SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
    7 
    8       USE filtreg_mod
    9       USE infotrac, ONLY : nqtot
    10       USE control_mod, ONLY: day_step,planet_type
     4SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
     5
     6  USE filtreg_mod
     7  USE infotrac, ONLY : nqtot
     8  USE control_mod, ONLY: day_step,planet_type
    119#ifdef CPP_IOIPSL
    12       USE IOIPSL
     10  USE IOIPSL
    1311#else
    14 ! if not using IOIPSL, we still need to use (a local version of) getin
    15       USE ioipsl_getincom
     12  ! if not using IOIPSL, we still need to use (a local version of) getin
     13  USE ioipsl_getincom
    1614#endif
    17       USE Write_Field
    18  
    19 
    20 c%W%    %G%
    21 c=======================================================================
    22 c
    23 c   Author:    Frederic Hourdin      original: 15/01/93
    24 c   -------
    25 c
    26 c   Subject:
    27 c   ------
    28 c
    29 c   Method:
    30 c   --------
    31 c
    32 c   Interface:
    33 c   ----------
    34 c
    35 c      Input:
    36 c      ------
    37 c
    38 c      Output:
    39 c      -------
    40 c
    41 c=======================================================================
    42       IMPLICIT NONE
    43 c-----------------------------------------------------------------------
    44 c   Declararations:
    45 c   ---------------
    46 
    47 #include "dimensions.h"
    48 #include "paramet.h"
    49 #include "comvert.h"
    50 #include "comconst.h"
    51 #include "comgeom.h"
    52 #include "academic.h"
    53 #include "ener.h"
    54 #include "temps.h"
    55 #include "iniprint.h"
    56 #include "logic.h"
    57 
    58 c   Arguments:
    59 c   ----------
    60 
    61       real time_0
    62 
    63 c   variables dynamiques
    64       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
    65       REAL teta(ip1jmp1,llm)                 ! temperature potentielle
    66       REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
    67       REAL ps(ip1jmp1)                       ! pression  au sol
    68       REAL masse(ip1jmp1,llm)                ! masse d'air
    69       REAL phis(ip1jmp1)                     ! geopotentiel au sol
    70 
    71 c   Local:
    72 c   ------
    73 
    74       REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
    75       REAL pks(ip1jmp1)                      ! exner au  sol
    76       REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
    77       REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
    78       REAL phi(ip1jmp1,llm)                  ! geopotentiel
    79       REAL ddsin,tetastrat,zsig,tetapv,w_pv  ! variables auxiliaires
    80       real tetajl(jjp1,llm)
    81       INTEGER i,j,l,lsup,ij
    82 
    83       REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
    84       REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
    85       LOGICAL ok_geost             ! Initialisation vent geost. ou nul
    86       LOGICAL ok_pv                ! Polar Vortex
    87       REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex
    88      
    89       real zz,ran1
    90       integer idum
    91 
    92       REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
    93 
    94 c-----------------------------------------------------------------------
    95 ! 1. Initializations for Earth-like case
    96 ! --------------------------------------
    97 c
    98 
    99         print *, 'This is iniacademic'
    100 
    101         ! initialize planet radius, rotation rate,...
    102         call conf_planete
    103 
    104         time_0=0.
    105         day_ref=1
    106         annee_ref=0
    107 
    108         im         = iim
    109         jm         = jjm
    110         day_ini    = 1
    111         dtvr    = daysec/REAL(day_step)
    112         zdtvr=dtvr
    113         etot0      = 0.
    114         ptot0      = 0.
    115         ztot0      = 0.
    116         stot0      = 0.
    117         ang0       = 0.
    118 
    119         if (llm.eq.1) then
    120           ! specific initializations for the shallow water case
    121           kappa=1
     15  USE Write_Field
     16
     17  !   Author:    Frederic Hourdin      original: 15/01/93
     18  ! The forcing defined here is from Held and Suarez, 1994, Bulletin
     19  ! of the American Meteorological Society, 75, 1825.
     20
     21  IMPLICIT NONE
     22
     23  !   Declararations:
     24  !   ---------------
     25
     26  include "dimensions.h"
     27  include "paramet.h"
     28  include "comvert.h"
     29  include "comconst.h"
     30  include "comgeom.h"
     31  include "academic.h"
     32  include "ener.h"
     33  include "temps.h"
     34  include "iniprint.h"
     35  include "logic.h"
     36
     37  !   Arguments:
     38  !   ----------
     39
     40  real time_0
     41
     42  !   variables dynamiques
     43  REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
     44  REAL teta(ip1jmp1,llm)                 ! temperature potentielle
     45  REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
     46  REAL ps(ip1jmp1)                       ! pression  au sol
     47  REAL masse(ip1jmp1,llm)                ! masse d'air
     48  REAL phis(ip1jmp1)                     ! geopotentiel au sol
     49
     50  !   Local:
     51  !   ------
     52
     53  REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
     54  REAL pks(ip1jmp1)                      ! exner au  sol
     55  REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
     56  REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
     57  REAL phi(ip1jmp1,llm)                  ! geopotentiel
     58  REAL ddsin,zsig,tetapv,w_pv  ! variables auxiliaires
     59  real tetastrat ! potential temperature in the stratosphere, in K
     60  real tetajl(jjp1,llm)
     61  INTEGER i,j,l,lsup,ij
     62
     63  REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
     64  REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
     65  LOGICAL ok_geost             ! Initialisation vent geost. ou nul
     66  LOGICAL ok_pv                ! Polar Vortex
     67  REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex
     68
     69  real zz,ran1
     70  integer idum
     71
     72  REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
     73
     74  !-----------------------------------------------------------------------
     75  ! 1. Initializations for Earth-like case
     76  ! --------------------------------------
     77  !
     78 
     79  write(lunout,*) 'Iniacademic'
     80 
     81  ! initialize planet radius, rotation rate,...
     82  call conf_planete
     83
     84  time_0=0.
     85  day_ref=1
     86  annee_ref=0
     87
     88  im         = iim
     89  jm         = jjm
     90  day_ini    = 1
     91  dtvr    = daysec/REAL(day_step)
     92  zdtvr=dtvr
     93  etot0      = 0.
     94  ptot0      = 0.
     95  ztot0      = 0.
     96  stot0      = 0.
     97  ang0       = 0.
     98
     99  if (llm == 1) then
     100     ! specific initializations for the shallow water case
     101     kappa=1
     102  endif
     103
     104  CALL iniconst
     105  CALL inigeom
     106  CALL inifilr
     107
     108  if (llm == 1) then
     109     ! initialize fields for the shallow water case, if required
     110     if (.not.read_start) then
     111        phis(:)=0.
     112        q(:,:,:)=0
     113        CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
     114     endif
     115  endif
     116
     117  academic_case: if (iflag_phys == 2) then
     118     ! initializations
     119
     120     ! 1. local parameters
     121     ! by convention, winter is in the southern hemisphere
     122     ! Geostrophic wind or no wind?
     123     ok_geost=.TRUE.
     124     CALL getin('ok_geost',ok_geost)
     125     ! Constants for Newtonian relaxation and friction
     126     k_f=1.                !friction
     127     CALL getin('k_j',k_f)
     128     k_f=1./(daysec*k_f)
     129     k_c_s=4.  !cooling surface
     130     CALL getin('k_c_s',k_c_s)
     131     k_c_s=1./(daysec*k_c_s)
     132     k_c_a=40. !cooling free atm
     133     CALL getin('k_c_a',k_c_a)
     134     k_c_a=1./(daysec*k_c_a)
     135     ! Constants for Teta equilibrium profile
     136     teta0=315.     ! mean Teta (S.H. 315K)
     137     CALL getin('teta0',teta0)
     138     write(lunout,*) 'Iniacademic - teta0 ',teta0
     139     write(lunout,*) 'Iniacademic - rad ',rad
     140     ttp=200.       ! Tropopause temperature (S.H. 200K)
     141     CALL getin('ttp',ttp)
     142     eps=0.         ! Deviation to N-S symmetry(~0-20K)
     143     CALL getin('eps',eps)
     144     delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
     145     CALL getin('delt_y',delt_y)
     146     delt_z=10.     ! Vertical Gradient (S.H. 10K)
     147     CALL getin('delt_z',delt_z)
     148     ! Polar vortex
     149     ok_pv=.false.
     150     CALL getin('ok_pv',ok_pv)
     151     phi_pv=-50.            ! Latitude of edge of vortex
     152     CALL getin('phi_pv',phi_pv)
     153     phi_pv=phi_pv*pi/180.
     154     dphi_pv=5.             ! Width of the edge
     155     CALL getin('dphi_pv',dphi_pv)
     156     dphi_pv=dphi_pv*pi/180.
     157     gam_pv=4.              ! -dT/dz vortex (in K/km)
     158     CALL getin('gam_pv',gam_pv)
     159
     160     ! 2. Initialize fields towards which to relax
     161     ! Friction
     162     knewt_g=k_c_a
     163     DO l=1,llm
     164        zsig=presnivs(l)/preff
     165        knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
     166        kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
     167     ENDDO
     168     DO j=1,jjp1
     169        clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
     170     ENDDO
     171
     172     ! Potential temperature
     173     DO l=1,llm
     174        zsig=presnivs(l)/preff
     175        tetastrat=ttp*zsig**(-kappa)
     176        tetapv=tetastrat
     177        IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
     178           tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
     179        ENDIF
     180        DO j=1,jjp1
     181           ! Troposphere
     182           ddsin=sin(rlatu(j))
     183           tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin &
     184                -delt_z*(1.-ddsin*ddsin)*log(zsig)
     185           !! Aymeric -- tests particuliers
     186           if (planet_type=="giant") then
     187             tetajl(j,l)=teta0+(delt_y*                   &
     188                ((sin(rlatu(j)*3.14159*eps+0.0001))**2)   &
     189                / ((rlatu(j)*3.14159*eps+0.0001)**2))     &
     190                -delt_z*log(zsig)
     191           endif
     192           ! Profil stratospherique isotherme (+vortex)
     193           w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
     194           tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
     195           tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 
     196        ENDDO
     197     ENDDO
     198
     199     !          CALL writefield('theta_eq',tetajl)
     200
     201     do l=1,llm
     202        do j=1,jjp1
     203           do i=1,iip1
     204              ij=(j-1)*iip1+i
     205              tetarappel(ij,l)=tetajl(j,l)
     206           enddo
     207        enddo
     208     enddo
     209     write(lunout,*) 'Iniacademic - check',tetajl(:,int(llm/2)),rlatu(:)
     210
     211     ! 3. Initialize fields (if necessary)
     212     IF (.NOT. read_start) THEN
     213        ! surface pressure
     214        ps(:)=preff
     215        ! ground geopotential
     216        phis(:)=0.
     217
     218        CALL pression ( ip1jmp1, ap, bp, ps, p       )
     219        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     220        CALL massdair(p,masse)
     221
     222        ! bulk initialization of temperature
     223        teta(:,:)=tetarappel(:,:)
     224
     225        ! geopotential
     226        CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
     227
     228        ! winds
     229        if (ok_geost) then
     230           call ugeostr(phi,ucov)
     231        else
     232           ucov(:,:)=0.
    122233        endif
    123        
    124         CALL iniconst
    125         CALL inigeom
    126         CALL inifilr
    127 
    128         if (llm.eq.1) then
    129           ! initialize fields for the shallow water case, if required
    130           if (.not.read_start) then
    131             phis(:)=0.
    132             q(:,:,:)=0
    133             CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
    134           endif
    135         endif
    136 
    137         if (iflag_phys.eq.2) then
    138           ! initializations for the academic case
    139          
    140 !         if (planet_type=="earth") then
    141 
    142           ! 1. local parameters
    143           ! by convention, winter is in the southern hemisphere
    144           ! Geostrophic wind or no wind?
    145           ok_geost=.TRUE.
    146           CALL getin('ok_geost',ok_geost)
    147           ! Constants for Newtonian relaxation and friction
    148           k_f=1.                !friction
    149           CALL getin('k_j',k_f)
    150           k_f=1./(daysec*k_f)
    151           k_c_s=4.  !cooling surface
    152           CALL getin('k_c_s',k_c_s)
    153           k_c_s=1./(daysec*k_c_s)
    154           k_c_a=40. !cooling free atm
    155           CALL getin('k_c_a',k_c_a)
    156           k_c_a=1./(daysec*k_c_a)
    157           ! Constants for Teta equilibrium profile
    158           teta0=315.     ! mean Teta (S.H. 315K)
    159           CALL getin('teta0',teta0)
    160           print *, 'iniacademic - teta0 ', teta0
    161           print *, 'iniacademic - rad ', rad
    162           ttp=200.       ! Tropopause temperature (S.H. 200K)
    163           CALL getin('ttp',ttp)
    164           eps=0.         ! Deviation to N-S symmetry(~0-20K)
    165           CALL getin('eps',eps)
    166           delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
    167           CALL getin('delt_y',delt_y)
    168           delt_z=10.     ! Vertical Gradient (S.H. 10K)
    169           CALL getin('delt_z',delt_z)
    170           ! Polar vortex
    171           ok_pv=.false.
    172           CALL getin('ok_pv',ok_pv)
    173           phi_pv=-50.            ! Latitude of edge of vortex
    174           CALL getin('phi_pv',phi_pv)
    175           phi_pv=phi_pv*pi/180.
    176           dphi_pv=5.             ! Width of the edge
    177           CALL getin('dphi_pv',dphi_pv)
    178           dphi_pv=dphi_pv*pi/180.
    179           gam_pv=4.              ! -dT/dz vortex (in K/km)
    180           CALL getin('gam_pv',gam_pv)
    181          
    182           ! 2. Initialize fields towards which to relax
    183           ! Friction
    184           knewt_g=k_c_a
    185           DO l=1,llm
    186            zsig=presnivs(l)/preff
    187            knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
    188            kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
    189           ENDDO
    190           DO j=1,jjp1
    191             clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
    192           ENDDO
    193          
    194           ! Potential temperature
    195           DO l=1,llm
    196            zsig=presnivs(l)/preff
    197            tetastrat=ttp*zsig**(-kappa)
    198            tetapv=tetastrat
    199            IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
    200                tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
    201            ENDIF
    202            DO j=1,jjp1
    203              ! Troposphere
    204              ddsin=sin(rlatu(j))
    205              tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin
    206      &           -delt_z*(1.-ddsin*ddsin)*log(zsig)
    207              !! Aymeric -- tests particuliers
    208              if (planet_type=="giant") then
    209              tetajl(j,l)=teta0+(delt_y*
    210      &          ((sin(rlatu(j)*3.14159*eps+0.0001))**2)
    211      &          / ((rlatu(j)*3.14159*eps+0.0001)**2))
    212      &          -delt_z*log(zsig)
    213 !!!             ddsin=sin(2.5*3.14159*rlatu(j))
    214 !!!             tetajl(j,l)=teta0-delt_y*ddsin*ddsin
    215 !!!!     &           -delt_z*(1.-ddsin*ddsin)*log(zsig)
    216              endif
    217              ! Profil stratospherique isotherme (+vortex)
    218              w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
    219              tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
    220              tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 
    221            ENDDO
    222           ENDDO ! of DO l=1,llm
    223  
    224 !          CALL writefield('theta_eq',tetajl)
    225 
    226           do l=1,llm
    227             do j=1,jjp1
    228               do i=1,iip1
    229                  ij=(j-1)*iip1+i
    230                  tetarappel(ij,l)=tetajl(j,l)
    231               enddo
    232             enddo
    233           enddo
    234           PRINT *, 'iniacademic - check',tetajl(:,int(llm/2)),rlatu(:)
    235 
    236 
    237 !         else
    238 !          write(lunout,*)"iniacademic: planet types other than earth",
    239 !     &                   " not implemented (yet)."
    240 !          stop
    241 !         endif ! of if (planet_type=="earth")
    242 
    243           ! 3. Initialize fields (if necessary)
    244           IF (.NOT. read_start) THEN
    245             ! surface pressure
    246             ps(:)=preff
    247             ! ground geopotential
    248             phis(:)=0.
    249            
    250             CALL pression ( ip1jmp1, ap, bp, ps, p       )
    251             CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
    252             CALL massdair(p,masse)
    253 
    254             ! bulk initialization of temperature
    255             teta(:,:)=tetarappel(:,:)
    256            
    257             ! geopotential
    258             CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
    259            
    260             ! winds
    261             if (ok_geost) then
    262               call ugeostr(phi,ucov)
    263             else
    264               ucov(:,:)=0.
    265             endif
    266             vcov(:,:)=0.
    267            
    268             ! bulk initialization of tracers
    269             if (planet_type=="earth") then
    270               ! Earth: first two tracers will be water
    271               do i=1,nqtot
    272                 if (i.eq.1) q(:,:,i)=1.e-10
    273                 if (i.eq.2) q(:,:,i)=1.e-15
    274                 if (i.gt.2) q(:,:,i)=0.
    275               enddo
    276             else
    277               q(:,:,:)=0
    278             endif ! of if (planet_type=="earth")
    279 
    280             ! add random perturbation to temperature
    281             idum  = -1
    282             zz = ran1(idum)
    283             idum  = 0
    284             do l=1,llm
    285               do ij=iip2,ip1jm
    286                 teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
    287               enddo
    288             enddo
    289 
    290             ! maintain periodicity in longitude
    291             do l=1,llm
    292               do ij=1,ip1jmp1,iip1
    293                 teta(ij+iim,l)=teta(ij,l)
    294               enddo
    295             enddo
    296 
    297           ENDIF ! of IF (.NOT. read_start)
    298         endif ! of if (iflag_phys.eq.2)
    299        
    300       END
    301 c-----------------------------------------------------------------------
     234        vcov(:,:)=0.
     235
     236        ! bulk initialization of tracers
     237        if (planet_type=="earth") then
     238           ! Earth: first two tracers will be water
     239           do i=1,nqtot
     240              if (i == 1) q(:,:,i)=1.e-10
     241              if (i == 2) q(:,:,i)=1.e-15
     242              if (i.gt.2) q(:,:,i)=0.
     243           enddo
     244        else
     245           q(:,:,:)=0
     246        endif ! of if (planet_type=="earth")
     247
     248        ! add random perturbation to temperature
     249        idum  = -1
     250        zz = ran1(idum)
     251        idum  = 0
     252        do l=1,llm
     253           do ij=iip2,ip1jm
     254              teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
     255           enddo
     256        enddo
     257
     258        ! maintain periodicity in longitude
     259        do l=1,llm
     260           do ij=1,ip1jmp1,iip1
     261              teta(ij+iim,l)=teta(ij,l)
     262           enddo
     263        enddo
     264
     265     ENDIF ! of IF (.NOT. read_start)
     266  endif academic_case
     267
     268END SUBROUTINE iniacademic

  • trunk/libf/dyn3d/ugeostr.F90

    r64 r66  
    11!
    2 ! $Id: ugeostr.F 1403 2010-07-01 09:02:53Z fairhead $
     2! $Id: ugeostr.F90 1474 2011-01-14 11:04:45Z lguez $
    33!
    4       subroutine ugeostr(phi,ucov)
     4subroutine ugeostr(phi,ucov)
    55
     6  ! Calcul du vent covariant geostrophique a partir du champ de
     7  ! geopotentiel.
     8  ! We actually compute: (1 - cos^8 \phi) u_g
     9  ! to have a wind going smoothly to 0 at the equator.
     10  ! We assume that the surface pressure is uniform so that model
     11  ! levels are pressure levels.
    612
    7 c  Calcul du vent covariant geostrophique a partir du champs de
    8 c  geopotentiel en supposant que le vent au sol est nul.
     13  implicit none
    914
    10       implicit none
     15  include "dimensions.h"
     16  include "paramet.h"
     17  include "comconst.h"
     18  include "comgeom2.h"
    1119
    12 #include "dimensions.h"
    13 #include "paramet.h"
    14 #include "comconst.h"
    15 #include "comgeom2.h"
     20  real ucov(iip1,jjp1,llm),phi(iip1,jjp1,llm)
     21  real um(jjm,llm),fact,u(iip1,jjm,llm)
     22  integer i,j,l
    1623
    17       real ucov(iip1,jjp1,llm),phi(iip1,jjp1,llm)
    18       real um(jjm,llm),fact,u(iip1,jjm,llm)
    19       integer i,j,l
     24  real zlat
    2025
    21       real zlat
     26  um(:,:)=0 ! initialize um()
    2227
    23       um(:,:)=0 ! initialize um()
     28  DO j=1,jjm
    2429
    25       DO j=1,jjm
     30     if (abs(sin(rlatv(j))).lt.1.e-4) then
     31        zlat=1.e-4
     32     else
     33        zlat=rlatv(j)
     34     endif
     35     fact=cos(zlat)
     36     fact=fact*fact
     37     fact=fact*fact
     38     fact=fact*fact
     39     fact=(1.-fact)/ &
     40          (2.*omeg*sin(zlat)*(rlatu(j+1)-rlatu(j)))
     41     fact=-fact/rad
     42     DO l=1,llm
     43        DO i=1,iim
     44           u(i,j,l)=fact*(phi(i,j+1,l)-phi(i,j,l))
     45           um(j,l)=um(j,l)+u(i,j,l)/REAL(iim)
     46        ENDDO
     47     ENDDO
     48  ENDDO
     49  call dump2d(jjm,llm,um,'Vent-u geostrophique')
    2650
    27          if (abs(sin(rlatv(j))).lt.1.e-4) then
    28              zlat=1.e-4
    29          else
    30              zlat=rlatv(j)
    31          endif
    32          fact=cos(zlat)
    33          fact=fact*fact
    34          fact=fact*fact
    35          fact=fact*fact
    36          fact=(1.-fact)/
    37      s    (2.*omeg*sin(zlat)*(rlatu(j+1)-rlatu(j)))
    38          fact=-fact/rad
    39          DO l=1,llm
    40             DO i=1,iim
    41                u(i,j,l)=fact*(phi(i,j+1,l)-phi(i,j,l))
    42                um(j,l)=um(j,l)+u(i,j,l)/REAL(iim)
    43             ENDDO
    44          ENDDO
    45       ENDDO
    46       call dump2d(jjm,llm,um,'Vent-u geostrophique')
     51  !   calcul des champ de vent:
    4752
    48 c
    49 c-----------------------------------------------------------------------
    50 c   calcul des champ de vent:
    51 c   -------------------------
     53  DO l=1,llm
     54     DO i=1,iip1
     55        ucov(i,1,l)=0.
     56        ucov(i,jjp1,l)=0.
     57     end DO
     58     DO  j=2,jjm
     59        DO  i=1,iim
     60           ucov(i,j,l) = 0.5*(u(i,j,l)+u(i,j-1,l))*cu(i,j)
     61        end DO
     62        ucov(iip1,j,l)=ucov(1,j,l)
     63     end DO
     64  end DO
    5265
    53       DO 301 l=1,llm
    54          DO 302 i=1,iip1
    55             ucov(i,1,l)=0.
    56             ucov(i,jjp1,l)=0.
    57 302      CONTINUE
    58          DO 304 j=2,jjm
    59             DO 305 i=1,iim
    60                ucov(i,j,l) = 0.5*(u(i,j,l)+u(i,j-1,l))*cu(i,j)
    61 305         CONTINUE
    62             ucov(iip1,j,l)=ucov(1,j,l)
    63 304      CONTINUE
    64 301   CONTINUE
     66  print *, 301
    6567
    66       print*,301
    67 
    68       return
    69       end
     68end subroutine ugeostr

  • trunk/libf/dyn3dpar/abort_gcm.F

    r1 r66  
    11!
    2 ! $Id: abort_gcm.F 1425 2010-09-02 13:45:23Z lguez $
     2! $Id: abort_gcm.F 1475 2011-01-21 14:41:03Z emillour $
    33!
    44c
    … …  
    4545      if (ierr .eq. 0) then
    4646        write(lunout,*) 'Everything is cool'
    47         stop
    4847      else
    4948        write(lunout,*) 'Houston, we have a problem ', ierr

  • trunk/libf/dyn3dpar/ce0l.F90

    r1 r66  
    11!
    2 ! $Id: ce0l.F90 1425 2010-09-02 13:45:23Z lguez $
     2! $Id: ce0l.F90 1482 2011-02-09 15:03:09Z jghattas $
    33!
    44!-------------------------------------------------------------------------------
    … …  
    2222  USE mod_const_mpi
    2323  USE infotrac
     24  USE parallel, ONLY: finalize_parallel
    2425
    2526#ifdef CPP_IOIPSL
    … …  
    5556       CALL abort_gcm('ce0l','In parallel mode,                         &
    5657 &                 ce0l must be called only                             &
    57  &                 for 1 process and 1 task')
     58 &                 for 1 process and 1 task',1)
    5859  ENDIF
    5960
    … …  
    101102  END IF
    102103
     104!$OMP MASTER
     105  CALL finalize_parallel
     106!$OMP END MASTER
     107
    103108#endif
    104109! of #ifndef CPP_EARTH #else

  • trunk/libf/dyn3dpar/disvert.F90

    r64 r66  
    1 !
    2 ! $Id: disvert.F 1403 2010-07-01 09:02:53Z fairhead $
    3 !
    4       SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
     1! $Id: disvert.F90 1480 2011-01-31 21:29:58Z jghattas $
    52
    6 c    Auteur :  P. Le Van .
    7 c
    8       IMPLICIT NONE
     3SUBROUTINE disvert(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
    94
    10 #include "dimensions.h"
    11 #include "paramet.h"
    12 #include "iniprint.h"
    13 #include "logic.h"
    14 c
    15 c=======================================================================
    16 c
    17 c
    18 c    s = sigma ** kappa   :  coordonnee  verticale
    19 c    dsig(l)            : epaisseur de la couche l ds la coord.  s
    20 c    sig(l)             : sigma a l'interface des couches l et l-1
    21 c    ds(l)              : distance entre les couches l et l-1 en coord.s
    22 c
    23 c=======================================================================
    24 c
    25       REAL pa,preff
    26       REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
    27       REAL presnivs(llm)
    28 c
    29 c   declarations:
    30 c   -------------
    31 c
    32       REAL sig(llm+1),dsig(llm)
    33        real zzz(1:llm+1)
    34        real dzz(1:llm)
    35       real zk,zkm1,dzk1,dzk2,k0,k1
    36 c
    37       INTEGER l
    38       REAL snorm,dsigmin
    39       REAL alpha,beta,gama,delta,deltaz,h
    40       INTEGER np,ierr
    41       REAL pi,x
     5  ! Auteur : P. Le Van
    426
    43       REAL SSUM
    44 c
    45 c-----------------------------------------------------------------------
    46 c
    47       pi=2.*ASIN(1.)
     7  IMPLICIT NONE
    488
    49       OPEN(99,file='sigma.def',status='old',form='formatted',
    50      s   iostat=ierr)
     9  include "dimensions.h"
     10  include "paramet.h"
     11  include "iniprint.h"
     12  include "logic.h"
    5113
    52 c-----------------------------------------------------------------------
    53 c   cas 1 on lit les options dans sigma.def:
    54 c   ----------------------------------------
     14  ! s = sigma ** kappa : coordonnee verticale
     15  ! dsig(l) : epaisseur de la couche l ds la coord. s
     16  ! sig(l) : sigma a l'interface des couches l et l-1
     17  ! ds(l) : distance entre les couches l et l-1 en coord.s
    5518
    56       IF (ierr.eq.0) THEN
     19  REAL pa, preff
     20  REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1)
     21  REAL presnivs(llm)
    5722
    58       READ(99,*) h           ! hauteur d'echelle 8.
    59       READ(99,*) deltaz      ! epaiseur de la premiere couche 0.04
    60       READ(99,*) beta        ! facteur d'acroissement en haut 1.3
    61       READ(99,*) k0          ! nombre de couches dans la transition surf
    62       READ(99,*) k1          ! nombre de couches dans la transition haute
    63       CLOSE(99)
    64       alpha=deltaz/(llm*h)
    65       write(lunout,*)'h,alpha,k0,k1,beta'
     23  REAL sig(llm+1), dsig(llm)
     24  real zk, zkm1, dzk1, dzk2, k0, k1
    6625
    67 c     read(*,*) h,deltaz,beta,k0,k1 ! 8 0.04 4 20 1.2
     26  INTEGER l
     27  REAL dsigmin
     28  REAL alpha, beta, deltaz, h
     29  INTEGER iostat
     30  REAL pi, x
    6831
    69       alpha=deltaz/tanh(1./k0)*2.
    70       zkm1=0.
    71       sig(1)=1.
    72       do l=1,llm
    73         sig(l+1)=(cosh(l/k0))**(-alpha*k0/h)
    74      + *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
     32  !-----------------------------------------------------------------------
     33
     34  pi = 2 * ASIN(1.)
     35
     36  OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
     37
     38  IF (iostat == 0) THEN
     39     ! cas 1 on lit les options dans sigma.def:
     40     READ(99, *) h ! hauteur d'echelle 8.
     41     READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
     42     READ(99, *) beta ! facteur d'acroissement en haut 1.3
     43     READ(99, *) k0 ! nombre de couches dans la transition surf
     44     READ(99, *) k1 ! nombre de couches dans la transition haute
     45     CLOSE(99)
     46     alpha=deltaz/(llm*h)
     47     write(lunout, *)'disvert: h, alpha, k0, k1, beta'
     48
     49     alpha=deltaz/tanh(1./k0)*2.
     50     zkm1=0.
     51     sig(1)=1.
     52     do l=1, llm
     53        sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) &
     54             *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
    7555        zk=-h*log(sig(l+1))
    7656
    7757        dzk1=alpha*tanh(l/k0)
    7858        dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
    79         write(lunout,*)l,sig(l+1),zk,zk-zkm1,dzk1,dzk2
     59        write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2
    8060        zkm1=zk
    81       enddo
     61     enddo
    8262
    83       sig(llm+1)=0.
     63     sig(llm+1)=0.
    8464
    85 c
    86        DO 2  l = 1, llm
    87        dsig(l) = sig(l)-sig(l+1)
    88    2   CONTINUE
    89 c
     65     DO l = 1, llm
     66        dsig(l) = sig(l)-sig(l+1)
     67     end DO
     68  ELSE
     69     if (ok_strato) then
     70        if (llm==39) then
     71           dsigmin=0.3
     72        else if (llm==50) then
     73           dsigmin=1.
     74        else
     75           WRITE(LUNOUT,*)'disvert: ATTENTION discretisation z a ajuster'
     76           dsigmin=1.
     77        endif
     78        WRITE(LUNOUT,*) 'disvert: Discretisation verticale DSIGMIN=',dsigmin
     79     endif
    9080
    91       ELSE
    92 c-----------------------------------------------------------------------
    93 c   cas 2 ancienne discretisation (LMD5...):
    94 c   ----------------------------------------
     81     DO l = 1, llm
     82        x = 2*asin(1.) * (l - 0.5) / (llm + 1)
    9583
    96       WRITE(LUNOUT,*)'WARNING!!! Ancienne discretisation verticale'
     84        IF (ok_strato) THEN
     85           dsig(l) =(dsigmin + 7. * SIN(x)**2) &
     86                *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2
     87        ELSE
     88           dsig(l) = 1.0 + 7.0 * SIN(x)**2
     89        ENDIF
     90     ENDDO
     91     dsig = dsig / sum(dsig)
     92     sig(llm+1) = 0.
     93     DO l = llm, 1, -1
     94        sig(l) = sig(l+1) + dsig(l)
     95     ENDDO
     96  ENDIF
    9797
    98       if (ok_strato) then
    99          if (llm==39) then
    100             dsigmin=0.3
    101          else if (llm==50) then
    102             dsigmin=1.
    103          else
    104             WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster'
    105             dsigmin=1.
    106          endif
    107          WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin
    108       endif
     98  DO l=1, llm
     99     nivsigs(l) = REAL(l)
     100  ENDDO
    109101
    110       h=7.
    111       snorm  = 0.
    112       DO l = 1, llm
    113          x = 2.*asin(1.) * (REAL(l)-0.5) / REAL(llm+1)
     102  DO l=1, llmp1
     103     nivsig(l)= REAL(l)
     104  ENDDO
    114105
    115          IF (ok_strato) THEN
    116            dsig(l) =(dsigmin + 7.0 * SIN(x)**2)
    117      &            *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2       
    118          ELSE
    119            dsig(l) = 1.0 + 7.0 * SIN(x)**2
    120          ENDIF
     106  ! .... Calculs de ap(l) et de bp(l) ....
     107  ! ..... pa et preff sont lus sur les fichiers start par lectba .....
    121108
    122          snorm = snorm + dsig(l)
    123       ENDDO
    124       snorm = 1./snorm
    125       DO l = 1, llm
    126          dsig(l) = dsig(l)*snorm
    127       ENDDO
    128       sig(llm+1) = 0.
    129       DO l = llm, 1, -1
    130          sig(l) = sig(l+1) + dsig(l)
    131       ENDDO
     109  bp(llmp1) = 0.
    132110
    133       ENDIF
     111  DO l = 1, llm
     112     bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
     113     ap(l) = pa * ( sig(l) - bp(l) )
     114  ENDDO
    134115
     116  bp(1)=1.
     117  ap(1)=0.
    135118
    136       DO l=1,llm
    137         nivsigs(l) = REAL(l)
    138       ENDDO
     119  ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
    139120
    140       DO l=1,llmp1
    141         nivsig(l)= REAL(l)
    142       ENDDO
     121  write(lunout, *)'disvert: BP '
     122  write(lunout, *) bp
     123  write(lunout, *)'disvert: AP '
     124  write(lunout, *) ap
    143125
    144 c
    145 c    ....  Calculs  de ap(l) et de bp(l)  ....
    146 c    .........................................
    147 c
    148 c
    149 c   .....  pa et preff sont lus  sur les fichiers start par lectba  .....
    150 c
     126  write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
     127  write(lunout, *)'couches calcules pour une pression de surface =', preff
     128  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de'
     129  write(lunout, *)'8km'
     130  DO l = 1, llm
     131     dpres(l) = bp(l) - bp(l+1)
     132     presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
     133     write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
     134          log(preff/presnivs(l))*8. &
     135          , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ &
     136          max(ap(l+1)+bp(l+1)*preff, 1.e-10))
     137  ENDDO
    151138
    152       bp(llmp1) =   0.
     139  write(lunout, *) 'disvert: PRESNIVS '
     140  write(lunout, *) presnivs
    153141
    154       DO l = 1, llm
    155 cc
    156 ccc    ap(l) = 0.
    157 ccc    bp(l) = sig(l)
    158 
    159       bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
    160       ap(l) = pa * ( sig(l) - bp(l) )
    161 c
    162       ENDDO
    163 
    164       bp(1)=1.
    165       ap(1)=0.
    166 
    167       ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
    168 
    169       write(lunout,*)' BP '
    170       write(lunout,*)  bp
    171       write(lunout,*)' AP '
    172       write(lunout,*)  ap
    173 
    174       write(lunout,*)
    175      .'Niveaux de pressions approximatifs aux centres des'
    176       write(lunout,*)'couches calcules pour une pression de surface =',
    177      .                 preff
    178       write(lunout,*)
    179      .     'et altitudes equivalentes pour une hauteur d echelle de'
    180       write(lunout,*)'8km'
    181       DO l = 1, llm
    182        dpres(l) = bp(l) - bp(l+1)
    183        presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
    184        write(lunout,*)'PRESNIVS(',l,')=',presnivs(l),'    Z ~ ',
    185      .        log(preff/presnivs(l))*8.
    186      .  ,'   DZ ~ ',8.*log((ap(l)+bp(l)*preff)/
    187      .       max(ap(l+1)+bp(l+1)*preff,1.e-10))
    188       ENDDO
    189 
    190       write(lunout,*)' PRESNIVS '
    191       write(lunout,*)presnivs
    192 
    193       RETURN
    194       END
     142END SUBROUTINE disvert

  • trunk/libf/dyn3dpar/etat0_netcdf.F90

    r1 r66  
    11!
    2 ! $Id: etat0_netcdf.F90 1425 2010-09-02 13:45:23Z lguez $
     2! $Id: etat0_netcdf.F90 1486 2011-02-11 12:07:39Z fairhead $
    33!
    44!-------------------------------------------------------------------------------
    … …  
    9898  REAL    :: dummy
    9999  LOGICAL :: ok_newmicro, ok_journe, ok_mensuel, ok_instan, ok_hf
    100   LOGICAL :: ok_LES, ok_ade, ok_aie, aerosol_couple, new_aod
     100  LOGICAL :: ok_LES, ok_ade, ok_aie, aerosol_couple, new_aod, callstats
    101101  INTEGER :: iflag_radia, flag_aerosol
    102102  REAL    :: bl95_b0, bl95_b1, fact_cldcon, facttemps, ratqsbas, ratqshaut
    … …  
    130130!--- CONSTRUCT A GRID
    131131  CALL conf_phys(  ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES,     &
     132                   callstats,                                           &
    132133                   solarlong0,seuil_inversion,                          &
    133134                   fact_cldcon, facttemps,ok_newmicro,iflag_radia,      &

  • trunk/libf/dyn3dpar/fluxstokenc_p.F

    r1 r66  
    11!
    2 ! $Id: fluxstokenc_p.F 1403 2010-07-01 09:02:53Z fairhead $
     2! $Id: fluxstokenc_p.F 1454 2010-11-18 12:01:24Z fairhead $
    33!
    44      SUBROUTINE fluxstokenc_p(pbaru,pbarv,masse,teta,phi,phis,
    … …  
    5757
    5858c AC initialisations
    59 cym      pbarug(:,:)   = 0.
     59      pbarug(:,:)   = 0.
    6060cym      pbarvg(:,:,:) = 0.
    6161cym      wg(:,:)       = 0.

  • trunk/libf/dyn3dpar/friction_p.F

    r1 r66  
    3434
    3535! arguments:
    36       REAL,INTENT(out) :: ucov( iip1,jjp1,llm )
    37       REAL,INTENT(out) :: vcov( iip1,jjm,llm )
     36      REAL,INTENT(inout) :: ucov( iip1,jjp1,llm )
     37      REAL,INTENT(inout) :: vcov( iip1,jjm,llm )
    3838      REAL,INTENT(in) :: pdt ! time step
    3939

  • trunk/libf/dyn3dpar/infotrac.F90

    r37 r66  
    177177          ! Continue to read tracer.def
    178178          DO iq=1,nqtrue
    179              READ(90,999) hadv(iq),vadv(iq),tnom_0(iq)
     179             READ(90,*) hadv(iq),vadv(iq),tnom_0(iq)
    180180          END DO
    181181          CLOSE(90) 
    … …  
    234234          ! Continue to read tracer.def
    235235          DO iq=1,nqtrue
    236              READ(90,999) hadv(iq),vadv(iq),tnom_0(iq)
     236             READ(90,*) hadv(iq),vadv(iq),tnom_0(iq)
    237237          END DO
    238238          CLOSE(90) 
    … …  
    316316       tname(new_iq)= tnom_0(iq)
    317317       IF (iadv(new_iq)==0) THEN
    318           ttext(new_iq)=str1(1:len_trim(str1))
     318          ttext(new_iq)=trim(str1)
    319319       ELSE
    320           ttext(new_iq)=str1(1:len_trim(str1))//descrq(iadv(new_iq))
     320          ttext(new_iq)=trim(tnom_0(iq))//descrq(iadv(new_iq))
    321321       END IF
    322322
    … …  
    381381    DEALLOCATE(tracnam)
    382382
    383 999 FORMAT (i2,1x,i2,1x,a15)
    384 
    385383  END SUBROUTINE infotrac_init
    386384

  • trunk/libf/dyn3dpar/iniacademic.F90

    r64 r66  
    11!
    2 ! $Id: iniacademic.F 1446 2010-10-22 09:27:25Z emillour $
     2! $Id: iniacademic.F90 1474 2011-01-14 11:04:45Z lguez $
    33!
    4 c
    5 c
    6       SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
    7 
    8       USE filtreg_mod
    9       USE infotrac, ONLY : nqtot
    10       USE control_mod, ONLY: day_step,planet_type
     4SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
     5
     6  USE filtreg_mod
     7  USE infotrac, ONLY : nqtot
     8  USE control_mod, ONLY: day_step,planet_type
    119#ifdef CPP_IOIPSL
    12       USE IOIPSL
     10  USE IOIPSL
    1311#else
    14 ! if not using IOIPSL, we still need to use (a local version of) getin
    15       USE ioipsl_getincom
     12  ! if not using IOIPSL, we still need to use (a local version of) getin
     13  USE ioipsl_getincom
    1614#endif
    17       USE Write_Field
    18  
    19 
    20 c%W%    %G%
    21 c=======================================================================
    22 c
    23 c   Author:    Frederic Hourdin      original: 15/01/93
    24 c   -------
    25 c
    26 c   Subject:
    27 c   ------
    28 c
    29 c   Method:
    30 c   --------
    31 c
    32 c   Interface:
    33 c   ----------
    34 c
    35 c      Input:
    36 c      ------
    37 c
    38 c      Output:
    39 c      -------
    40 c
    41 c=======================================================================
    42       IMPLICIT NONE
    43 c-----------------------------------------------------------------------
    44 c   Declararations:
    45 c   ---------------
    46 
    47 #include "dimensions.h"
    48 #include "paramet.h"
    49 #include "comvert.h"
    50 #include "comconst.h"
    51 #include "comgeom.h"
    52 #include "academic.h"
    53 #include "ener.h"
    54 #include "temps.h"
    55 #include "iniprint.h"
    56 #include "logic.h"
    57 
    58 c   Arguments:
    59 c   ----------
    60 
    61       real time_0
    62 
    63 c   variables dynamiques
    64       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
    65       REAL teta(ip1jmp1,llm)                 ! temperature potentielle
    66       REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
    67       REAL ps(ip1jmp1)                       ! pression  au sol
    68       REAL masse(ip1jmp1,llm)                ! masse d'air
    69       REAL phis(ip1jmp1)                     ! geopotentiel au sol
    70 
    71 c   Local:
    72 c   ------
    73 
    74       REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
    75       REAL pks(ip1jmp1)                      ! exner au  sol
    76       REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
    77       REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
    78       REAL phi(ip1jmp1,llm)                  ! geopotentiel
    79       REAL ddsin,tetastrat,zsig,tetapv,w_pv  ! variables auxiliaires
    80       real tetajl(jjp1,llm)
    81       INTEGER i,j,l,lsup,ij
    82 
    83       REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
    84       REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
    85       LOGICAL ok_geost             ! Initialisation vent geost. ou nul
    86       LOGICAL ok_pv                ! Polar Vortex
    87       REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex
    88      
    89       real zz,ran1
    90       integer idum
    91 
    92       REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
    93 
    94 c-----------------------------------------------------------------------
    95 ! 1. Initializations for Earth-like case
    96 ! --------------------------------------
    97 c
    98 
    99         print *, 'This is iniacademic'
    100 
    101         ! initialize planet radius, rotation rate,...
    102         call conf_planete
    103 
    104         time_0=0.
    105         day_ref=1
    106         annee_ref=0
    107 
    108         im         = iim
    109         jm         = jjm
    110         day_ini    = 1
    111         dtvr    = daysec/REAL(day_step)
    112         zdtvr=dtvr
    113         etot0      = 0.
    114         ptot0      = 0.
    115         ztot0      = 0.
    116         stot0      = 0.
    117         ang0       = 0.
    118 
    119         if (llm.eq.1) then
    120           ! specific initializations for the shallow water case
    121           kappa=1
     15  USE Write_Field
     16
     17  !   Author:    Frederic Hourdin      original: 15/01/93
     18  ! The forcing defined here is from Held and Suarez, 1994, Bulletin
     19  ! of the American Meteorological Society, 75, 1825.
     20
     21  IMPLICIT NONE
     22
     23  !   Declararations:
     24  !   ---------------
     25
     26  include "dimensions.h"
     27  include "paramet.h"
     28  include "comvert.h"
     29  include "comconst.h"
     30  include "comgeom.h"
     31  include "academic.h"
     32  include "ener.h"
     33  include "temps.h"
     34  include "iniprint.h"
     35  include "logic.h"
     36
     37  !   Arguments:
     38  !   ----------
     39
     40  real time_0
     41
     42  !   variables dynamiques
     43  REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
     44  REAL teta(ip1jmp1,llm)                 ! temperature potentielle
     45  REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
     46  REAL ps(ip1jmp1)                       ! pression  au sol
     47  REAL masse(ip1jmp1,llm)                ! masse d'air
     48  REAL phis(ip1jmp1)                     ! geopotentiel au sol
     49
     50  !   Local:
     51  !   ------
     52
     53  REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
     54  REAL pks(ip1jmp1)                      ! exner au  sol
     55  REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
     56  REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
     57  REAL phi(ip1jmp1,llm)                  ! geopotentiel
     58  REAL ddsin,zsig,tetapv,w_pv  ! variables auxiliaires
     59  real tetastrat ! potential temperature in the stratosphere, in K
     60  real tetajl(jjp1,llm)
     61  INTEGER i,j,l,lsup,ij
     62
     63  REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
     64  REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
     65  LOGICAL ok_geost             ! Initialisation vent geost. ou nul
     66  LOGICAL ok_pv                ! Polar Vortex
     67  REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex
     68
     69  real zz,ran1
     70  integer idum
     71
     72  REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
     73
     74  !-----------------------------------------------------------------------
     75  ! 1. Initializations for Earth-like case
     76  ! --------------------------------------
     77  !
     78 
     79  write(lunout,*) 'Iniacademic'
     80 
     81  ! initialize planet radius, rotation rate,...
     82  call conf_planete
     83
     84  time_0=0.
     85  day_ref=1
     86  annee_ref=0
     87
     88  im         = iim
     89  jm         = jjm
     90  day_ini    = 1
     91  dtvr    = daysec/REAL(day_step)
     92  zdtvr=dtvr
     93  etot0      = 0.
     94  ptot0      = 0.
     95  ztot0      = 0.
     96  stot0      = 0.
     97  ang0       = 0.
     98
     99  if (llm == 1) then
     100     ! specific initializations for the shallow water case
     101     kappa=1
     102  endif
     103
     104  CALL iniconst
     105  CALL inigeom
     106  CALL inifilr
     107
     108  if (llm == 1) then
     109     ! initialize fields for the shallow water case, if required
     110     if (.not.read_start) then
     111        phis(:)=0.
     112        q(:,:,:)=0
     113        CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
     114     endif
     115  endif
     116
     117  academic_case: if (iflag_phys == 2) then
     118     ! initializations
     119
     120     ! 1. local parameters
     121     ! by convention, winter is in the southern hemisphere
     122     ! Geostrophic wind or no wind?
     123     ok_geost=.TRUE.
     124     CALL getin('ok_geost',ok_geost)
     125     ! Constants for Newtonian relaxation and friction
     126     k_f=1.                !friction
     127     CALL getin('k_j',k_f)
     128     k_f=1./(daysec*k_f)
     129     k_c_s=4.  !cooling surface
     130     CALL getin('k_c_s',k_c_s)
     131     k_c_s=1./(daysec*k_c_s)
     132     k_c_a=40. !cooling free atm
     133     CALL getin('k_c_a',k_c_a)
     134     k_c_a=1./(daysec*k_c_a)
     135     ! Constants for Teta equilibrium profile
     136     teta0=315.     ! mean Teta (S.H. 315K)
     137     CALL getin('teta0',teta0)
     138     write(lunout,*) 'Iniacademic - teta0 ',teta0
     139     write(lunout,*) 'Iniacademic - rad ',rad
     140     ttp=200.       ! Tropopause temperature (S.H. 200K)
     141     CALL getin('ttp',ttp)
     142     eps=0.         ! Deviation to N-S symmetry(~0-20K)
     143     CALL getin('eps',eps)
     144     delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
     145     CALL getin('delt_y',delt_y)
     146     delt_z=10.     ! Vertical Gradient (S.H. 10K)
     147     CALL getin('delt_z',delt_z)
     148     ! Polar vortex
     149     ok_pv=.false.
     150     CALL getin('ok_pv',ok_pv)
     151     phi_pv=-50.            ! Latitude of edge of vortex
     152     CALL getin('phi_pv',phi_pv)
     153     phi_pv=phi_pv*pi/180.
     154     dphi_pv=5.             ! Width of the edge
     155     CALL getin('dphi_pv',dphi_pv)
     156     dphi_pv=dphi_pv*pi/180.
     157     gam_pv=4.              ! -dT/dz vortex (in K/km)
     158     CALL getin('gam_pv',gam_pv)
     159
     160     ! 2. Initialize fields towards which to relax
     161     ! Friction
     162     knewt_g=k_c_a
     163     DO l=1,llm
     164        zsig=presnivs(l)/preff
     165        knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
     166        kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
     167     ENDDO
     168     DO j=1,jjp1
     169        clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
     170     ENDDO
     171
     172     ! Potential temperature
     173     DO l=1,llm
     174        zsig=presnivs(l)/preff
     175        tetastrat=ttp*zsig**(-kappa)
     176        tetapv=tetastrat
     177        IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
     178           tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
     179        ENDIF
     180        DO j=1,jjp1
     181           ! Troposphere
     182           ddsin=sin(rlatu(j))
     183           tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin &
     184                -delt_z*(1.-ddsin*ddsin)*log(zsig)
     185           !! Aymeric -- tests particuliers
     186           if (planet_type=="giant") then
     187             tetajl(j,l)=teta0+(delt_y*                   &
     188                ((sin(rlatu(j)*3.14159*eps+0.0001))**2)   &
     189                / ((rlatu(j)*3.14159*eps+0.0001)**2))     &
     190                -delt_z*log(zsig)
     191           endif
     192           ! Profil stratospherique isotherme (+vortex)
     193           w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
     194           tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
     195           tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 
     196        ENDDO
     197     ENDDO
     198
     199     !          CALL writefield('theta_eq',tetajl)
     200
     201     do l=1,llm
     202        do j=1,jjp1
     203           do i=1,iip1
     204              ij=(j-1)*iip1+i
     205              tetarappel(ij,l)=tetajl(j,l)
     206           enddo
     207        enddo
     208     enddo
     209     write(lunout,*) 'Iniacademic - check',tetajl(:,int(llm/2)),rlatu(:)
     210
     211     ! 3. Initialize fields (if necessary)
     212     IF (.NOT. read_start) THEN
     213        ! surface pressure
     214        ps(:)=preff
     215        ! ground geopotential
     216        phis(:)=0.
     217
     218        CALL pression ( ip1jmp1, ap, bp, ps, p       )
     219        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     220        CALL massdair(p,masse)
     221
     222        ! bulk initialization of temperature
     223        teta(:,:)=tetarappel(:,:)
     224
     225        ! geopotential
     226        CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
     227
     228        ! winds
     229        if (ok_geost) then
     230           call ugeostr(phi,ucov)
     231        else
     232           ucov(:,:)=0.
    122233        endif
    123        
    124         CALL iniconst
    125         CALL inigeom
    126         CALL inifilr
    127 
    128         if (llm.eq.1) then
    129           ! initialize fields for the shallow water case, if required
    130           if (.not.read_start) then
    131             phis(:)=0.
    132             q(:,:,:)=0
    133             CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
    134           endif
    135         endif
    136 
    137         if (iflag_phys.eq.2) then
    138           ! initializations for the academic case
    139          
    140 !         if (planet_type=="earth") then
    141 
    142           ! 1. local parameters
    143           ! by convention, winter is in the southern hemisphere
    144           ! Geostrophic wind or no wind?
    145           ok_geost=.TRUE.
    146           CALL getin('ok_geost',ok_geost)
    147           ! Constants for Newtonian relaxation and friction
    148           k_f=1.                !friction
    149           CALL getin('k_j',k_f)
    150           k_f=1./(daysec*k_f)
    151           k_c_s=4.  !cooling surface
    152           CALL getin('k_c_s',k_c_s)
    153           k_c_s=1./(daysec*k_c_s)
    154           k_c_a=40. !cooling free atm
    155           CALL getin('k_c_a',k_c_a)
    156           k_c_a=1./(daysec*k_c_a)
    157           ! Constants for Teta equilibrium profile
    158           teta0=315.     ! mean Teta (S.H. 315K)
    159           CALL getin('teta0',teta0)
    160           print *, 'iniacademic - teta0 ', teta0
    161           print *, 'iniacademic - rad ', rad
    162           ttp=200.       ! Tropopause temperature (S.H. 200K)
    163           CALL getin('ttp',ttp)
    164           eps=0.         ! Deviation to N-S symmetry(~0-20K)
    165           CALL getin('eps',eps)
    166           delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
    167           CALL getin('delt_y',delt_y)
    168           delt_z=10.     ! Vertical Gradient (S.H. 10K)
    169           CALL getin('delt_z',delt_z)
    170           ! Polar vortex
    171           ok_pv=.false.
    172           CALL getin('ok_pv',ok_pv)
    173           phi_pv=-50.            ! Latitude of edge of vortex
    174           CALL getin('phi_pv',phi_pv)
    175           phi_pv=phi_pv*pi/180.
    176           dphi_pv=5.             ! Width of the edge
    177           CALL getin('dphi_pv',dphi_pv)
    178           dphi_pv=dphi_pv*pi/180.
    179           gam_pv=4.              ! -dT/dz vortex (in K/km)
    180           CALL getin('gam_pv',gam_pv)
    181          
    182           ! 2. Initialize fields towards which to relax
    183           ! Friction
    184           knewt_g=k_c_a
    185           DO l=1,llm
    186            zsig=presnivs(l)/preff
    187            knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
    188            kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
    189           ENDDO
    190           DO j=1,jjp1
    191             clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
    192           ENDDO
    193          
    194           ! Potential temperature
    195           DO l=1,llm
    196            zsig=presnivs(l)/preff
    197            tetastrat=ttp*zsig**(-kappa)
    198            tetapv=tetastrat
    199            IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
    200                tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
    201            ENDIF
    202            DO j=1,jjp1
    203              ! Troposphere
    204              ddsin=sin(rlatu(j))
    205              tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin
    206      &           -delt_z*(1.-ddsin*ddsin)*log(zsig)
    207              !! Aymeric -- tests particuliers
    208              if (planet_type=="giant") then
    209              tetajl(j,l)=teta0+(delt_y*
    210      &          ((sin(rlatu(j)*3.14159*eps+0.0001))**2)
    211      &          / ((rlatu(j)*3.14159*eps+0.0001)**2))
    212      &          -delt_z*log(zsig)
    213 !!!             ddsin=sin(2.5*3.14159*rlatu(j))
    214 !!!             tetajl(j,l)=teta0-delt_y*ddsin*ddsin
    215 !!!!     &           -delt_z*(1.-ddsin*ddsin)*log(zsig)
    216              endif
    217              ! Profil stratospherique isotherme (+vortex)
    218              w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
    219              tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
    220              tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 
    221            ENDDO
    222           ENDDO ! of DO l=1,llm
    223  
    224 !          CALL writefield('theta_eq',tetajl)
    225 
    226           do l=1,llm
    227             do j=1,jjp1
    228               do i=1,iip1
    229                  ij=(j-1)*iip1+i
    230                  tetarappel(ij,l)=tetajl(j,l)
    231               enddo
    232             enddo
    233           enddo
    234           PRINT *, 'iniacademic - check',tetajl(:,int(llm/2)),rlatu(:)
    235 
    236 
    237 !         else
    238 !          write(lunout,*)"iniacademic: planet types other than earth",
    239 !     &                   " not implemented (yet)."
    240 !          stop
    241 !         endif ! of if (planet_type=="earth")
    242 
    243           ! 3. Initialize fields (if necessary)
    244           IF (.NOT. read_start) THEN
    245             ! surface pressure
    246             ps(:)=preff
    247             ! ground geopotential
    248             phis(:)=0.
    249            
    250             CALL pression ( ip1jmp1, ap, bp, ps, p       )
    251             CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
    252             CALL massdair(p,masse)
    253 
    254             ! bulk initialization of temperature
    255             teta(:,:)=tetarappel(:,:)
    256            
    257             ! geopotential
    258             CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
    259            
    260             ! winds
    261             if (ok_geost) then
    262               call ugeostr(phi,ucov)
    263             else
    264               ucov(:,:)=0.
    265             endif
    266             vcov(:,:)=0.
    267            
    268             ! bulk initialization of tracers
    269             if (planet_type=="earth") then
    270               ! Earth: first two tracers will be water
    271               do i=1,nqtot
    272                 if (i.eq.1) q(:,:,i)=1.e-10
    273                 if (i.eq.2) q(:,:,i)=1.e-15
    274                 if (i.gt.2) q(:,:,i)=0.
    275               enddo
    276             else
    277               q(:,:,:)=0
    278             endif ! of if (planet_type=="earth")
    279 
    280             ! add random perturbation to temperature
    281             idum  = -1
    282             zz = ran1(idum)
    283             idum  = 0
    284             do l=1,llm
    285               do ij=iip2,ip1jm
    286                 teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
    287               enddo
    288             enddo
    289 
    290             ! maintain periodicity in longitude
    291             do l=1,llm
    292               do ij=1,ip1jmp1,iip1
    293                 teta(ij+iim,l)=teta(ij,l)
    294               enddo
    295             enddo
    296 
    297           ENDIF ! of IF (.NOT. read_start)
    298         endif ! of if (iflag_phys.eq.2)
    299        
    300       END
    301 c-----------------------------------------------------------------------
     234        vcov(:,:)=0.
     235
     236        ! bulk initialization of tracers
     237        if (planet_type=="earth") then
     238           ! Earth: first two tracers will be water
     239           do i=1,nqtot
     240              if (i == 1) q(:,:,i)=1.e-10
     241              if (i == 2) q(:,:,i)=1.e-15
     242              if (i.gt.2) q(:,:,i)=0.
     243           enddo
     244        else
     245           q(:,:,:)=0
     246        endif ! of if (planet_type=="earth")
     247
     248        ! add random perturbation to temperature
     249        idum  = -1
     250        zz = ran1(idum)
     251        idum  = 0
     252        do l=1,llm
     253           do ij=iip2,ip1jm
     254              teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
     255           enddo
     256        enddo
     257
     258        ! maintain periodicity in longitude
     259        do l=1,llm
     260           do ij=1,ip1jmp1,iip1
     261              teta(ij+iim,l)=teta(ij,l)
     262           enddo
     263        enddo
     264
     265     ENDIF ! of IF (.NOT. read_start)
     266  endif academic_case
     267
     268END SUBROUTINE iniacademic

  • trunk/libf/dyn3dpar/parallel.F90

    r1 r66  
    11!
    2 ! $Id: parallel.F90 1279 2009-12-10 09:02:56Z fairhead $
     2! $Id: parallel.F90 1487 2011-02-11 15:07:54Z jghattas $
    33!
    44  module parallel
    55  USE mod_const_mpi
    66   
    7     LOGICAL,SAVE :: using_mpi
     7    LOGICAL,SAVE :: using_mpi=.TRUE.
    88    LOGICAL,SAVE :: using_omp
    99   
    … …  
    208208      integer :: ierr
    209209      integer :: i
    210       deallocate(jj_begin_para)
    211       deallocate(jj_end_para)
    212       deallocate(jj_nb_para)
     210
     211      if (allocated(jj_begin_para)) deallocate(jj_begin_para)
     212      if (allocated(jj_end_para))   deallocate(jj_end_para)
     213      if (allocated(jj_nb_para))    deallocate(jj_nb_para)
    213214
    214215      if (type_ocean == 'couple') then
    … …  
    549550       
    550551   
    551     /* 
    552   Subroutine verif_hallo(Field,ij,ll,up,down)
    553     implicit none
    554 #include "dimensions.h"
    555 #include "paramet.h"   
    556     include 'mpif.h'
    557    
    558       INTEGER :: ij,ll
    559       REAL, dimension(ij,ll) :: Field
    560       INTEGER :: up,down
    561      
    562       REAL,dimension(ij,ll): NewField
    563      
    564       NewField=0
    565      
    566       ijb=ij_begin
    567       ije=ij_end
    568       if (pole_nord)
    569       NewField(ij_be       
    570 */
     552!  Subroutine verif_hallo(Field,ij,ll,up,down)
     553!    implicit none
     554!#include "dimensions.h"
     555!#include "paramet.h"   
     556!    include 'mpif.h'
     557!   
     558!      INTEGER :: ij,ll
     559!      REAL, dimension(ij,ll) :: Field
     560!      INTEGER :: up,down
     561!     
     562!      REAL,dimension(ij,ll): NewField
     563!     
     564!      NewField=0
     565!     
     566!      ijb=ij_begin
     567!      ije=ij_end
     568!      if (pole_nord)
     569!      NewField(ij_be       
     570
    571571  end module parallel

  • trunk/libf/dyn3dpar/ugeostr.F90

    r64 r66  
    11!
    2 ! $Id: ugeostr.F 1403 2010-07-01 09:02:53Z fairhead $
     2! $Id: ugeostr.F90 1474 2011-01-14 11:04:45Z lguez $
    33!
    4       subroutine ugeostr(phi,ucov)
     4subroutine ugeostr(phi,ucov)
    55
     6  ! Calcul du vent covariant geostrophique a partir du champ de
     7  ! geopotentiel.
     8  ! We actually compute: (1 - cos^8 \phi) u_g
     9  ! to have a wind going smoothly to 0 at the equator.
     10  ! We assume that the surface pressure is uniform so that model
     11  ! levels are pressure levels.
    612
    7 c  Calcul du vent covariant geostrophique a partir du champs de
    8 c  geopotentiel en supposant que le vent au sol est nul.
     13  implicit none
    914
    10       implicit none
     15  include "dimensions.h"
     16  include "paramet.h"
     17  include "comconst.h"
     18  include "comgeom2.h"
    1119
    12 #include "dimensions.h"
    13 #include "paramet.h"
    14 #include "comconst.h"
    15 #include "comgeom2.h"
     20  real ucov(iip1,jjp1,llm),phi(iip1,jjp1,llm)
     21  real um(jjm,llm),fact,u(iip1,jjm,llm)
     22  integer i,j,l
    1623
    17       real ucov(iip1,jjp1,llm),phi(iip1,jjp1,llm)
    18       real um(jjm,llm),fact,u(iip1,jjm,llm)
    19       integer i,j,l
     24  real zlat
    2025
    21       real zlat
     26  um(:,:)=0 ! initialize um()
    2227
    23       um(:,:)=0 ! initialize um()
     28  DO j=1,jjm
    2429
    25       DO j=1,jjm
     30     if (abs(sin(rlatv(j))).lt.1.e-4) then
     31        zlat=1.e-4
     32     else
     33        zlat=rlatv(j)
     34     endif
     35     fact=cos(zlat)
     36     fact=fact*fact
     37     fact=fact*fact
     38     fact=fact*fact
     39     fact=(1.-fact)/ &
     40          (2.*omeg*sin(zlat)*(rlatu(j+1)-rlatu(j)))
     41     fact=-fact/rad
     42     DO l=1,llm
     43        DO i=1,iim
     44           u(i,j,l)=fact*(phi(i,j+1,l)-phi(i,j,l))
     45           um(j,l)=um(j,l)+u(i,j,l)/REAL(iim)
     46        ENDDO
     47     ENDDO
     48  ENDDO
     49  call dump2d(jjm,llm,um,'Vent-u geostrophique')
    2650
    27          if (abs(sin(rlatv(j))).lt.1.e-4) then
    28              zlat=1.e-4
    29          else
    30              zlat=rlatv(j)
    31          endif
    32          fact=cos(zlat)
    33          fact=fact*fact
    34          fact=fact*fact
    35          fact=fact*fact
    36          fact=(1.-fact)/
    37      s    (2.*omeg*sin(zlat)*(rlatu(j+1)-rlatu(j)))
    38          fact=-fact/rad
    39          DO l=1,llm
    40             DO i=1,iim
    41                u(i,j,l)=fact*(phi(i,j+1,l)-phi(i,j,l))
    42                um(j,l)=um(j,l)+u(i,j,l)/REAL(iim)
    43             ENDDO
    44          ENDDO
    45       ENDDO
    46       call dump2d(jjm,llm,um,'Vent-u geostrophique')
     51  !   calcul des champ de vent:
    4752
    48 c
    49 c-----------------------------------------------------------------------
    50 c   calcul des champ de vent:
    51 c   -------------------------
     53  DO l=1,llm
     54     DO i=1,iip1
     55        ucov(i,1,l)=0.
     56        ucov(i,jjp1,l)=0.
     57     end DO
     58     DO  j=2,jjm
     59        DO  i=1,iim
     60           ucov(i,j,l) = 0.5*(u(i,j,l)+u(i,j-1,l))*cu(i,j)
     61        end DO
     62        ucov(iip1,j,l)=ucov(1,j,l)
     63     end DO
     64  end DO
    5265
    53       DO 301 l=1,llm
    54          DO 302 i=1,iip1
    55             ucov(i,1,l)=0.
    56             ucov(i,jjp1,l)=0.
    57 302      CONTINUE
    58          DO 304 j=2,jjm
    59             DO 305 i=1,iim
    60                ucov(i,j,l) = 0.5*(u(i,j,l)+u(i,j-1,l))*cu(i,j)
    61 305         CONTINUE
    62             ucov(iip1,j,l)=ucov(1,j,l)
    63 304      CONTINUE
    64 301   CONTINUE
     66  print *, 301
    6567
    66       print*,301
    67 
    68       return
    69       end
     68end subroutine ugeostr

  • trunk/libf/phylmd/aeropt_5wv.F90

    r1 r66  
    11!
    2 ! $Id: aeropt_5wv.F90 1418 2010-07-19 15:11:24Z jghattas $
     2! $Id: aeropt_5wv.F90 1470 2010-12-21 15:51:32Z idelkadi $
    33!
    44
    … …  
    756756    ENDIF
    757757
    758     used_tau(spsol)=.TRUE.
     758!Bug 21 12 10 AI
     759!    used_tau(spsol)=.TRUE.
     760    IF (soluble) then
     761      used_tau(spsol)=.TRUE.
     762       ELSE
     763      used_tau(naero_soluble+spinsol)=.TRUE.
     764    ENDIF
     765
    759766    DO la=1,las
    760767

  • trunk/libf/phylmd/calcul_divers.h

    r1 r66  
    22c $Header$
    33c
    4 c
    5 c initialisations diverses au "debut" du mois
    6 c
    7       IF(MOD(itap,INT(ecrit_mth/dtime)).EQ.1) THEN
     4
     5c     Initialisations diverses au "debut" du mois
     6      IF(debut) THEN
     7         nday_rain(:)=0.
     8
     9c        surface terre
     10         paire_ter(:)=0.
    811         DO i=1, klon
    9           nday_rain(i)=0.
    10          ENDDO !i
    11 c
    12 c surface terre
    13        DO i=1, klon
    14          IF(pctsrf(i,is_ter).GT.0.) THEN
    15             paire_ter(i)=airephy(i)*pctsrf(i,is_ter)
    16          ENDIF
    17        ENDDO
    18 c
    19       ENDIF !MOD(itap,INT(ecrit_mth)).EQ.1
    20 c
    21       IF(MOD(itap,INT(ecrit_day/dtime)).EQ.0) THEN
    22 c
    23 cIM calcul total_rain, nday_rain
    24 c
    25        DO i = 1, klon
    26         total_rain(i)=rain_fall(i)+snow_fall(i) 
    27         IF(total_rain(i).GT.0.) nday_rain(i)=nday_rain(i)+1.
    28        ENDDO
    29       ENDIF !itap.EQ.ecrit_mth
     12            IF(pctsrf(i,is_ter).GT.0.) THEN
     13               paire_ter(i)=airephy(i)*pctsrf(i,is_ter)
     14            ENDIF
     15         ENDDO
     16      ENDIF
     17
     18cIM   Calcul une fois par jour : total_rain, nday_rain
     19      IF(MOD(itap,INT(un_jour/dtime)).EQ.0) THEN
     20         DO i = 1, klon
     21            total_rain(i)=rain_fall(i)+snow_fall(i) 
     22            IF(total_rain(i).GT.0.) nday_rain(i)=nday_rain(i)+1.
     23         ENDDO
     24      ENDIF

  • trunk/libf/phylmd/carbon_cycle_mod.F90

    r1 r66  
    11MODULE carbon_cycle_mod
    2 
     2! Controle module for the carbon CO2 tracers :
     3!   - Identification
     4!   - Get concentrations comming from coupled model or read from file to tracers
     5!   - Calculate new RCO2 for radiation scheme
     6!   - Calculate new carbon flux for sending to coupled models (PISCES and ORCHIDEE)
     7!
    38! Author : Josefine GHATTAS, Patricia CADULE
    49
    … …  
    1318  LOGICAL, PUBLIC :: carbon_cycle_cpl       ! Coupling of CO2 fluxes between LMDZ/ORCHIDEE and LMDZ/OCEAN(PISCES)
    1419!$OMP THREADPRIVATE(carbon_cycle_cpl)
     20
    1521  LOGICAL :: carbon_cycle_emis_comp=.FALSE. ! Calculation of emission compatible
     22!$OMP THREADPRIVATE(carbon_cycle_emis_comp)
     23
     24  LOGICAL :: RCO2_inter  ! RCO2 interactive : if true calculate new value RCO2 for the radiation scheme
     25!$OMP THREADPRIVATE(RCO2_inter)
    1626
    1727! Scalare values when no transport, from physiq.def
    … …  
    2131!$OMP THREADPRIVATE(emis_land_s)
    2232
    23   INTEGER :: ntr_co2                ! Number of tracers concerning the carbon cycle
    24   INTEGER :: id_fco2_tot            ! Tracer index
    25   INTEGER :: id_fco2_ocn            !  - " -
    26   INTEGER :: id_fco2_land           !  - " -
    27   INTEGER :: id_fco2_land_use       !  - " -
    28   INTEGER :: id_fco2_fos_fuel       !  - " -
    29 !$OMP THREADPRIVATE(ntr_co2, id_fco2_tot, id_fco2_ocn, id_fco2_land, id_fco2_land_use, id_fco2_fos_fuel) 
    30 
    31   REAL, DIMENSION(:), ALLOCATABLE :: fos_fuel        ! CO2 fossil fuel emission from file [gC/m2/d]
    32 !$OMP THREADPRIVATE(fos_fuel)
    33   REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ocn_day ! flux CO2 from ocean for 1 day (cumulated) [gC/m2/d]
     33  REAL :: airetot     ! Total area of the earth surface
     34!$OMP THREADPRIVATE(airetot)
     35
     36  INTEGER :: ntr_co2  ! Number of tracers concerning the carbon cycle
     37!$OMP THREADPRIVATE(ntr_co2)
     38
     39! fco2_ocn_day : flux CO2 from ocean for 1 day (cumulated) [gC/m2/d]. Allocation and initalization done in cpl_mod
     40  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ocn_day
    3441!$OMP THREADPRIVATE(fco2_ocn_day)
     42
    3543  REAL, DIMENSION(:), ALLOCATABLE :: fco2_land_day   ! flux CO2 from land for 1 day (cumulated)  [gC/m2/d]
    3644!$OMP THREADPRIVATE(fco2_land_day)
    … …  
    3846!$OMP THREADPRIVATE(fco2_lu_day)
    3947
    40 ! Following 2 fields will be initialized in surf_land_orchidee at each time step
     48  REAL, DIMENSION(:,:), ALLOCATABLE :: dtr_add       ! Tracer concentration to be injected
     49!$OMP THREADPRIVATE(dtr_add)
     50
     51! Following 2 fields will be allocated and initialized in surf_land_orchidee
    4152  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_inst  ! flux CO2 from land at one time step
    4253!$OMP THREADPRIVATE(fco2_land_inst)
    … …  
    4556
    4657! Calculated co2 field to be send to the ocean via the coupler and to ORCHIDEE
    47   REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: co2_send
     58  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: co2_send ! Field allocated in phyetat0
    4859!$OMP THREADPRIVATE(co2_send)
     60
     61
     62  TYPE, PUBLIC ::   co2_trac_type
     63     CHARACTER(len = 8) :: name       ! Tracer name in tracer.def
     64     INTEGER            :: id         ! Index in total tracer list, tr_seri
     65     CHARACTER(len=30)  :: file       ! File name
     66     LOGICAL            :: cpl        ! True if this tracers is coupled from ORCHIDEE or PISCES.
     67                                      ! False if read from file.
     68     INTEGER            :: updatefreq ! Frequence to inject in second
     69     INTEGER            :: readstep   ! Actual time step to read in file
     70     LOGICAL            :: updatenow  ! True if this tracer should be updated this time step
     71  END TYPE co2_trac_type
     72  INTEGER,PARAMETER :: maxco2trac=5  ! Maximum number of different CO2 fluxes
     73  TYPE(co2_trac_type), DIMENSION(maxco2trac) :: co2trac
    4974
    5075CONTAINS
    5176 
    52   SUBROUTINE carbon_cycle_init(tr_seri, aerosol, radio)
     77  SUBROUTINE carbon_cycle_init(tr_seri, pdtphys, aerosol, radio)
     78! This subroutine is called from traclmdz_init, only at first timestep.
     79! - Read controle parameters from .def input file
     80! - Search for carbon tracers and set default values
     81! - Allocate variables
     82! - Test for compatibility
     83
    5384    USE dimphy
     85    USE comgeomphy
     86    USE mod_phys_lmdz_transfert_para
    5487    USE infotrac
    5588    USE IOIPSL
    5689    USE surface_data, ONLY : ok_veget, type_ocean
     90    USE phys_cal_mod, ONLY : mth_len
    5791
    5892    IMPLICIT NONE
    … …  
    6296! Input argument
    6397    REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr_seri ! Concentration Traceur [U/KgA] 
     98    REAL,INTENT(IN)                           :: pdtphys ! length of time step in physiq (sec)
    6499
    65100! InOutput arguments
    … …  
    68103
    69104! Local variables
    70     INTEGER               :: ierr, it, iiq
    71     REAL, DIMENSION(klon) :: tr_seri_sum
    72 
    73 
    74 ! 0) Test for compatibility
    75     IF (carbon_cycle_cpl .AND. type_ocean/='couple') &
    76          CALL abort_gcm('carbon_cycle_init', 'Coupling with ocean model is needed for carbon_cycle_cpl',1)
    77     IF (carbon_cycle_cpl .AND..NOT. ok_veget) &
    78          CALL abort_gcm('carbon_cycle_init', 'Coupling with surface land model ORCHDIEE is needed for carbon_cycle_cpl',1)
    79 
    80 
    81 ! 1) Check if transport of one tracer flux CO2 or 4 separated tracers
    82     IF (carbon_cycle_tr) THEN
    83        id_fco2_tot=0
    84        id_fco2_ocn=0
    85        id_fco2_land=0
    86        id_fco2_land_use=0
    87        id_fco2_fos_fuel=0
    88        
    89        ! Search in tracer list
    90        DO it=1,nbtr
    91           iiq=niadv(it+2)
    92           IF (tname(iiq) == "fCO2" ) THEN
    93              id_fco2_tot=it
    94           ELSE IF (tname(iiq) == "fCO2_ocn" ) THEN
    95              id_fco2_ocn=it
    96           ELSE IF (tname(iiq) == "fCO2_land" ) THEN
    97              id_fco2_land=it
    98           ELSE IF (tname(iiq) == "fCO2_land_use" ) THEN
    99              id_fco2_land_use=it
    100           ELSE IF (tname(iiq) == "fCO2_fos_fuel" ) THEN
    101              id_fco2_fos_fuel=it
    102           END IF
    103        END DO
    104 
    105        ! Count tracers found
    106        IF (id_fco2_tot /= 0 .AND. &
    107             id_fco2_ocn==0 .AND. id_fco2_land==0 .AND. id_fco2_land_use==0 .AND. id_fco2_fos_fuel==0) THEN
    108          
    109           ! transport  1 tracer flux CO2
    110           ntr_co2 = 1
    111          
    112        ELSE IF (id_fco2_tot==0 .AND. &
    113             id_fco2_ocn /=0 .AND. id_fco2_land/=0 .AND. id_fco2_land_use/=0 .AND. id_fco2_fos_fuel/=0) THEN
    114          
    115           ! transport 4 tracers seperatively
    116           ntr_co2 = 4
    117          
    118        ELSE
    119           CALL abort_gcm('carbon_cycle_init', 'error in coherence between traceur.def and gcm.def',1)
    120        END IF
    121 
    122        ! Definition of control varaiables for the tracers
    123        DO it=1,nbtr
    124           IF (it==id_fco2_tot .OR. it==id_fco2_ocn .OR. it==id_fco2_land .OR. &
    125                it==id_fco2_land_use .OR. it==id_fco2_fos_fuel) THEN
    126              aerosol(it) = .FALSE.
    127              radio(it)   = .FALSE.
    128           END IF
    129        END DO
    130 
    131     ELSE
    132        ! No transport of CO2
    133        ntr_co2 = 0
    134     END IF ! carbon_cycle_tr
    135 
    136 
    137 ! 2) Allocate variable for CO2 fossil fuel emission
    138     IF (carbon_cycle_tr) THEN
    139        ! Allocate 2D variable
    140        ALLOCATE(fos_fuel(klon), stat=ierr)
    141        IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 1',1)
    142     ELSE
    143        ! No transport : read value from .def
     105    INTEGER               :: ierr, it, iiq, itc
     106    INTEGER               :: teststop
     107
     108
     109
     110! 1) Read controle parameters from .def input file
     111! ------------------------------------------------
     112    ! Read fosil fuel value if no transport
     113    IF (.NOT. carbon_cycle_tr) THEN
    144114       fos_fuel_s = 0.
    145115       CALL getin ('carbon_cycle_fos_fuel',fos_fuel_s)
    … …  
    148118
    149119
    150 ! 3) Allocate and initialize fluxes
    151     IF (carbon_cycle_cpl) THEN
    152        IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 2',1)
    153        ALLOCATE(fco2_land_day(klon), stat=ierr)
    154        IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 3',1)
    155        ALLOCATE(fco2_lu_day(klon), stat=ierr)
    156        IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 4',1)
    157 
    158        fco2_land_day(:) = 0.  ! JG : Doit prend valeur de restart
    159        fco2_lu_day(:)   = 0.  ! JG : Doit prend valeur de restart
    160 
    161        ! fco2_ocn_day is allocated in cpl_mod
    162        ! fco2_land_inst and fco2_lu_inst are allocated in surf_land_orchidee
    163        
    164        ALLOCATE(co2_send(klon), stat=ierr)
    165        IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 7',1)
    166        
    167        ! Calculate using restart tracer values
    168        IF (carbon_cycle_tr) THEN
    169           IF (ntr_co2==1) THEN
    170              co2_send(:) = tr_seri(:,1,id_fco2_tot) + co2_ppm0
    171           ELSE ! ntr_co2==4
    172              ! Calculate the delta CO2 flux
    173              tr_seri_sum(:) = tr_seri(:,1,id_fco2_fos_fuel) + tr_seri(:,1,id_fco2_land_use) + &
    174                   tr_seri(:,1,id_fco2_land) + tr_seri(:,1,id_fco2_ocn)
    175              co2_send(:) = tr_seri_sum(:) + co2_ppm0
    176           END IF
    177        ELSE
    178           ! Send a scalare value in 2D variable to ocean and land model (PISCES and ORCHIDEE)
    179           co2_send(:) = co2_ppm
    180        END IF
    181 
    182 
    183     ELSE
    184        IF (carbon_cycle_tr) THEN
    185           ! No coupling of CO2 fields :
    186           ! corresponding fields will instead be read from files
    187           ALLOCATE(fco2_ocn_day(klon), stat=ierr)
    188           IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 8',1)
    189           ALLOCATE(fco2_land_day(klon), stat=ierr)
    190           IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 9',1)
    191           ALLOCATE(fco2_lu_day(klon), stat=ierr)
    192           IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 10',1)       
    193        END IF
    194     END IF
    195 
    196 ! 4) Read parmeter for calculation of emission compatible
     120    ! Read parmeter for calculation compatible emission
    197121    IF (.NOT. carbon_cycle_tr) THEN
    198122       carbon_cycle_emis_comp=.FALSE.
    199123       CALL getin('carbon_cycle_emis_comp',carbon_cycle_emis_comp)
    200124       WRITE(lunout,*) 'carbon_cycle_emis_comp = ',carbon_cycle_emis_comp
    201     END IF
     125       IF (carbon_cycle_emis_comp) THEN
     126          CALL abort_gcm('carbon_cycle_init', 'carbon_cycle_emis_comp option not yet implemented!!',1)
     127       END IF
     128    END IF
     129
     130    ! Read parameter for interactive calculation of the CO2 value for the radiation scheme
     131    RCO2_inter=.FALSE.
     132    CALL getin('RCO2_inter',RCO2_inter)
     133    WRITE(lunout,*) 'RCO2_inter = ', RCO2_inter
     134    IF (RCO2_inter) THEN
     135       WRITE(lunout,*) 'RCO2 will be recalculated once a day'
     136       WRITE(lunout,*) 'RCO2 initial = ', RCO2
     137    END IF
     138
     139
     140! 2) Search for carbon tracers and set default values
     141! ---------------------------------------------------
     142    itc=0
     143    DO it=1,nbtr
     144       iiq=niadv(it+2)
     145       
     146       SELECT CASE(tname(iiq))
     147       CASE("fCO2_ocn")
     148          itc = itc + 1
     149          co2trac(itc)%name='fCO2_ocn'
     150          co2trac(itc)%id=it
     151          co2trac(itc)%file='fl_co2_ocean.nc'
     152          IF (carbon_cycle_cpl .AND. type_ocean=='couple') THEN
     153             co2trac(itc)%cpl=.TRUE.
     154             co2trac(itc)%updatefreq = 86400 ! Once a day as the coupling with OASIS/PISCES
     155          ELSE
     156             co2trac(itc)%cpl=.FALSE.
     157             co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
     158          END IF
     159       CASE("fCO2_land")
     160          itc = itc + 1
     161          co2trac(itc)%name='fCO2_land'
     162          co2trac(itc)%id=it
     163          co2trac(itc)%file='fl_co2_land.nc'
     164          IF (carbon_cycle_cpl .AND. ok_veget) THEN
     165             co2trac(itc)%cpl=.TRUE.
     166             co2trac(itc)%updatefreq = INT(pdtphys) ! Each timestep as the coupling with ORCHIDEE
     167          ELSE
     168             co2trac(itc)%cpl=.FALSE.
     169!             co2trac(itc)%updatefreq = 10800   ! 10800sec = 3H
     170             co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
     171          END IF
     172       CASE("fCO2_land_use")
     173          itc = itc + 1
     174          co2trac(itc)%name='fCO2_land_use'
     175          co2trac(itc)%id=it
     176          co2trac(itc)%file='fl_co2_land_use.nc'
     177          IF (carbon_cycle_cpl .AND. ok_veget) THEN
     178             co2trac(it)%cpl=.TRUE.
     179             co2trac(itc)%updatefreq = INT(pdtphys) ! Each timestep as the coupling with ORCHIDEE
     180          ELSE
     181             co2trac(itc)%cpl=.FALSE.
     182             co2trac(itc)%updatefreq = 10800   ! 10800sec = 3H
     183          END IF
     184       CASE("fCO2_fos_fuel")
     185          itc = itc + 1
     186          co2trac(itc)%name='fCO2_fos_fuel'
     187          co2trac(itc)%id=it
     188          co2trac(itc)%file='fossil_fuel.nc'
     189          co2trac(itc)%cpl=.FALSE.       ! This tracer always read from file
     190!         co2trac(itc)%updatefreq = 86400  ! 86400sec = 24H Cadule case
     191          co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
     192       CASE("fCO2_bbg")
     193          itc = itc + 1
     194          co2trac(itc)%name='fCO2_bbg'
     195          co2trac(itc)%id=it
     196          co2trac(itc)%file='fl_co2_bbg.nc'
     197          co2trac(itc)%cpl=.FALSE.       ! This tracer always read from file
     198          co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
     199       CASE("fCO2")
     200          ! fCO2 : One tracer transporting the total CO2 flux
     201          itc = itc + 1
     202          co2trac(itc)%name='fCO2'
     203          co2trac(itc)%id=it
     204          co2trac(itc)%file='fl_co2.nc'
     205          IF (carbon_cycle_cpl) THEN
     206             co2trac(itc)%cpl=.TRUE.
     207          ELSE
     208             co2trac(itc)%cpl=.FALSE.
     209          END IF
     210          co2trac(itc)%updatefreq = 86400
     211          ! DOES THIS WORK ???? Problematic due to implementation of the coupled fluxes...
     212          CALL abort_gcm('carbon_cycle_init','transport of total CO2 has to be implemented and tested',1)
     213       END SELECT
     214    END DO
     215
     216    ! Total number of carbon CO2 tracers
     217    ntr_co2 = itc
     218   
     219    ! Definition of control varaiables for the tracers
     220    DO it=1,ntr_co2
     221       aerosol(co2trac(it)%id) = .FALSE.
     222       radio(co2trac(it)%id)   = .FALSE.
     223    END DO
     224   
     225    ! Vector indicating which timestep to read for each tracer
     226    ! Always start read in the beginning of the file
     227    co2trac(:)%readstep = 0
     228   
     229
     230! 3) Allocate variables
     231! ---------------------
     232    ! Allocate vector for storing fluxes to inject
     233    ALLOCATE(dtr_add(klon,maxco2trac), stat=ierr)
     234    IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 11',1)       
     235   
     236    ! Allocate variables for cumulating fluxes from ORCHIDEE
     237    IF (RCO2_inter) THEN
     238       IF (.NOT. carbon_cycle_tr .AND. carbon_cycle_cpl) THEN
     239          ALLOCATE(fco2_land_day(klon), stat=ierr)
     240          IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 2',1)
     241          fco2_land_day(1:klon) = 0.
     242         
     243          ALLOCATE(fco2_lu_day(klon), stat=ierr)
     244          IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 3',1)
     245          fco2_lu_day(1:klon)   = 0.
     246       END IF
     247    END IF
     248
     249
     250! 4) Test for compatibility
     251! -------------------------
     252!    IF (carbon_cycle_cpl .AND. type_ocean/='couple') THEN
     253!       WRITE(lunout,*) 'Coupling with ocean model is needed for carbon_cycle_cpl'
     254!       CALL abort_gcm('carbon_cycle_init', 'coupled ocean is needed for carbon_cycle_cpl',1)
     255!    END IF
     256!
     257!    IF (carbon_cycle_cpl .AND..NOT. ok_veget) THEN
     258!       WRITE(lunout,*) 'Coupling with surface land model ORCHDIEE is needed for carbon_cycle_cpl'
     259!       CALL abort_gcm('carbon_cycle_init', 'ok_veget is needed for carbon_cycle_cpl',1)
     260!    END IF
     261
     262    ! Compiler test : following should never happen
     263    teststop=0
     264    DO it=1,teststop
     265       CALL abort_gcm('carbon_cycle_init', 'Entering loop from 1 to 0',1)
     266    END DO
     267
     268    IF (ntr_co2==0) THEN
     269       ! No carbon tracers found in tracer.def. It is not possible to do carbon cycle
     270       WRITE(lunout,*) 'No carbon tracers found in tracer.def. Not ok with carbon_cycle_tr and/or carbon_cycle_cp'
     271       CALL abort_gcm('carbon_cycle_init', 'No carbon tracers found in tracer.def',1)
     272    END IF
     273   
     274! 5) Calculate total area of the earth surface
     275! --------------------------------------------
     276    CALL reduce_sum(SUM(airephy),airetot)
     277    CALL bcast(airetot)
    202278
    203279  END SUBROUTINE carbon_cycle_init
    204280
     281
     282  SUBROUTINE carbon_cycle(nstep, pdtphys, pctsrf, tr_seri, source)
     283! Subroutine for injection of co2 in the tracers
    205284!
    206 !
    207 !
    208 
    209   SUBROUTINE carbon_cycle(nstep, pdtphys, pctsrf, tr_seri)
    210    
     285! - Find out if it is time to update
     286! - Get tracer from coupled model or from file
     287! - Calculate new RCO2 value for the radiation scheme
     288! - Calculate CO2 flux to send to ocean and land models (PISCES and ORCHIDEE)
     289
    211290    USE infotrac
    212291    USE dimphy
    213     USE mod_phys_lmdz_transfert_para, ONLY : reduce_sum
     292    USE mod_phys_lmdz_transfert_para
    214293    USE phys_cal_mod, ONLY : mth_cur, mth_len
    215294    USE phys_cal_mod, ONLY : day_cur
    … …  
    220299    INCLUDE "clesphys.h"
    221300    INCLUDE "indicesol.h"
     301    INCLUDE "iniprint.h"
     302    INCLUDE "YOMCST.h"
    222303
    223304! In/Output arguments
    224305    INTEGER,INTENT(IN) :: nstep      ! time step in physiq
    225306    REAL,INTENT(IN)    :: pdtphys    ! length of time step in physiq (sec)
    226     REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
    227     REAL, DIMENSION(klon,klev,nbtr), INTENT(INOUT)  :: tr_seri
     307    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf            ! Surface fraction
     308    REAL, DIMENSION(klon,klev,nbtr), INTENT(INOUT)  :: tr_seri ! All tracers
     309    REAL, DIMENSION(klon,nbtr), INTENT(INOUT)       :: source  ! Source for all tracers
    228310
    229311! Local variables
     312    INTEGER :: it
    230313    LOGICAL :: newmonth ! indicates if a new month just started
    231314    LOGICAL :: newday   ! indicates if a new day just started
    … …  
    233316
    234317    REAL, PARAMETER :: fact=1.E-15/2.12  ! transformation factor from gC/m2/day => ppm/m2/day
    235     REAL, DIMENSION(klon) :: fco2_tmp, tr_seri_sum
     318    REAL, DIMENSION(klon) :: fco2_tmp
    236319    REAL :: sumtmp
    237     REAL :: airetot     ! Total area the earth
    238320    REAL :: delta_co2_ppm
    239321   
    240 ! -) Calculate logicals indicating if it is a new month, new day or the last time step in a day (end day)
     322
     323! 1) Calculate logicals indicating if it is a new month, new day or the last time step in a day (end day)
     324! -------------------------------------------------------------------------------------------------------
    241325
    242326    newday = .FALSE.; endday = .FALSE.; newmonth = .FALSE.
    … …  
    245329    IF (MOD(nstep,INT(86400./pdtphys))==0) endday=.TRUE.
    246330    IF (newday .AND. day_cur==1) newmonth=.TRUE.
    247    
    248 ! -) Read new maps if new month started
    249     IF (newmonth .AND. carbon_cycle_tr) THEN
    250        CALL read_map2D('fossil_fuel.nc','fos_fuel', mth_cur, .FALSE., fos_fuel)
    251        
    252        ! division by month lenght to get dayly value
    253        fos_fuel(:) = fos_fuel(:)/mth_len
    254        
    255        IF (.NOT. carbon_cycle_cpl) THEN
    256           ! Get dayly values from monthly fluxes
    257           CALL read_map2D('fl_co2_ocean.nc','CO2_OCN',mth_cur,.FALSE.,fco2_ocn_day)
    258           CALL read_map2D('fl_co2_land.nc','CO2_LAND', mth_cur,.FALSE.,fco2_land_day)
    259           CALL read_map2D('fl_co2_land_use.nc','CO2_LAND_USE',mth_cur,.FALSE.,fco2_lu_day)
    260        END IF
    261     END IF
    262    
    263 
    264 
    265 ! -) Update tracers at beginning of a new day. Beginning of a new day correspond to a new coupling period in cpl_mod.
    266     IF (newday) THEN
     331
     332! 2)  For each carbon tracer find out if it is time to inject (update)
     333! --------------------------------------------------------------------
     334    DO it = 1, ntr_co2
     335       IF ( MOD(nstep,INT(co2trac(it)%updatefreq/pdtphys)) == 1 ) THEN
     336          co2trac(it)%updatenow = .TRUE.
     337       ELSE
     338          co2trac(it)%updatenow = .FALSE.
     339       END IF
     340    END DO
     341
     342! 3) Get tracer update
     343! --------------------------------------
     344    DO it = 1, ntr_co2
     345       IF ( co2trac(it)%updatenow ) THEN
     346          IF ( co2trac(it)%cpl ) THEN
     347             ! Get tracer from coupled model
     348             SELECT CASE(co2trac(it)%name)
     349             CASE('fCO2_land')     ! from ORCHIDEE
     350                dtr_add(:,it) = fco2_land_inst(:)*pctsrf(:,is_ter)*fact ! [ppm/m2/day]
     351             CASE('fCO2_land_use') ! from ORCHIDEE
     352                dtr_add(:,it) = fco2_lu_inst(:)  *pctsrf(:,is_ter)*fact ! [ppm/m2/day]
     353             CASE('fCO2_ocn')      ! from PISCES
     354                dtr_add(:,it) = fco2_ocn_day(:)  *pctsrf(:,is_oce)*fact ! [ppm/m2/day]
     355             CASE DEFAULT
     356                WRITE(lunout,*) 'Error with tracer ',co2trac(it)%name
     357                CALL abort_gcm('carbon_cycle', 'No coupling implemented for this tracer',1)
     358             END SELECT
     359          ELSE
     360             ! Read tracer from file
     361             co2trac(it)%readstep = co2trac(it)%readstep + 1 ! increment time step in file
     362! Patricia   CALL read_map2D(co2trac(it)%file,'fco2',co2trac(it)%readstep,.FALSE.,dtr_add(:,it))
     363             CALL read_map2D(co2trac(it)%file,'fco2',co2trac(it)%readstep,.TRUE.,dtr_add(:,it))
     364
     365             ! Converte from kgC/m2/h to kgC/m2/s
     366             dtr_add(:,it) = dtr_add(:,it)/3600
     367             ! Add individual treatment of values read from file
     368             SELECT CASE(co2trac(it)%name)
     369             CASE('fCO2_land')
     370                dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_ter)
     371             CASE('fCO2_land_use')
     372                dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_ter)
     373             CASE('fCO2_ocn')
     374                dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_oce)
     375! Patricia :
     376!             CASE('fCO2_fos_fuel')
     377!                dtr_add(:,it) = dtr_add(:,it)/mth_len
     378!                co2trac(it)%readstep = 0 ! Always read same value for fossil fuel(Cadule case)
     379             END SELECT
     380          END IF
     381       END IF
     382    END DO
     383
     384! 4) Update co2 tracers :
     385!    Loop over all carbon tracers and add source
     386! ------------------------------------------------------------------
     387    IF (carbon_cycle_tr) THEN
     388       DO it = 1, ntr_co2
     389          IF (.FALSE.) THEN
     390             tr_seri(1:klon,1,co2trac(it)%id) = tr_seri(1:klon,1,co2trac(it)%id) + dtr_add(1:klon,it)
     391             source(1:klon,co2trac(it)%id) = 0.
     392          ELSE
     393             source(1:klon,co2trac(it)%id) = dtr_add(1:klon,it)
     394          END IF
     395       END DO
     396    END IF
     397
     398
     399! 5) Calculations for new CO2 value for the radiation scheme(instead of reading value from .def)
     400! ----------------------------------------------------------------------------------------------
     401    IF (RCO2_inter) THEN
     402       ! Cumulate fluxes from ORCHIDEE at each timestep
     403       IF (.NOT. carbon_cycle_tr .AND. carbon_cycle_cpl) THEN
     404          IF (newday) THEN ! Reset cumulative variables once a day
     405             fco2_land_day(1:klon) = 0.
     406             fco2_lu_day(1:klon)   = 0.
     407          END IF
     408          fco2_land_day(1:klon) = fco2_land_day(1:klon) + fco2_land_inst(1:klon) ![gC/m2/day]
     409          fco2_lu_day(1:klon)   = fco2_lu_day(1:klon)   + fco2_lu_inst(1:klon)   ![gC/m2/day]
     410       END IF
     411
     412       ! At the end of a new day, calculate a mean scalare value of CO2
     413       ! JG : Ici on utilise uniquement le traceur du premier couche du modele. Est-ce que c'est correcte ?
     414       IF (endday) THEN
     415
     416          IF (carbon_cycle_tr) THEN
     417             ! Sum all co2 tracers to get the total delta CO2 flux
     418             fco2_tmp(:) = 0.
     419             DO it = 1, ntr_co2
     420                fco2_tmp(1:klon) = fco2_tmp(1:klon) + tr_seri(1:klon,1,co2trac(it)%id)
     421             END DO
     422             
     423          ELSE IF (carbon_cycle_cpl) THEN ! no carbon_cycle_tr
     424             ! Sum co2 fluxes comming from coupled models and parameter for fossil fuel
     425             fco2_tmp(1:klon) = fos_fuel_s + ((fco2_lu_day(1:klon) + fco2_land_day(1:klon))*pctsrf(1:klon,is_ter) &
     426                  + fco2_ocn_day(:)*pctsrf(:,is_oce)) * fact
     427          END IF
     428
     429          ! Calculate a global mean value of delta CO2 flux
     430          fco2_tmp(1:klon) = fco2_tmp(1:klon) * airephy(1:klon)
     431          CALL reduce_sum(SUM(fco2_tmp),sumtmp)
     432          CALL bcast(sumtmp)
     433          delta_co2_ppm = sumtmp/airetot
     434         
     435          ! Add initial value for co2_ppm and delta value
     436          co2_ppm = co2_ppm0 + delta_co2_ppm
     437         
     438          ! Transformation of atmospheric CO2 concentration for the radiation code
     439          RCO2 = co2_ppm * 1.0e-06  * 44.011/28.97
     440         
     441          WRITE(lunout,*) 'RCO2 is now updated! RCO2 = ', RCO2
     442       END IF ! endday
     443
     444    END IF ! RCO2_inter
     445
     446
     447! 6) Calculate CO2 flux to send to ocean and land models : PISCES and ORCHIDEE         
     448! ----------------------------------------------------------------------------
     449    IF (carbon_cycle_cpl) THEN
    267450
    268451       IF (carbon_cycle_tr) THEN
    269 
    270           ! Update tracers
    271           IF (ntr_co2 == 1) THEN
    272              ! Calculate the new flux CO2
    273              tr_seri(:,1,id_fco2_tot) = tr_seri(:,1,id_fco2_tot) + &
    274                   (fos_fuel(:) + &
    275                   fco2_lu_day(:)  * pctsrf(:,is_ter) + &
    276                   fco2_land_day(:)* pctsrf(:,is_ter) + &
    277                   fco2_ocn_day(:) * pctsrf(:,is_oce)) * fact
    278 
    279           ELSE ! ntr_co2 == 4
    280              tr_seri(:,1,id_fco2_fos_fuel)  = tr_seri(:,1,id_fco2_fos_fuel) + fos_fuel(:) * fact ! [ppm/m2/day]
    281 
    282              tr_seri(:,1,id_fco2_land_use)  = tr_seri(:,1,id_fco2_land_use) + &
    283                   fco2_lu_day(:)  *pctsrf(:,is_ter)*fact ! [ppm/m2/day]
    284 
    285              tr_seri(:,1,id_fco2_land)      = tr_seri(:,1,id_fco2_land)     + &
    286                   fco2_land_day(:)*pctsrf(:,is_ter)*fact ! [ppm/m2/day]
    287 
    288              tr_seri(:,1,id_fco2_ocn)       = tr_seri(:,1,id_fco2_ocn)      + &
    289                   fco2_ocn_day(:) *pctsrf(:,is_oce)*fact ! [ppm/m2/day]
    290 
    291           END IF
    292 
    293        ELSE ! no transport
    294           IF (carbon_cycle_cpl) THEN
    295              IF (carbon_cycle_emis_comp) THEN
    296                 ! Calcul emission compatible a partir des champs 2D et co2_ppm
    297                 !!! TO DO!!
    298                 CALL abort_gcm('carbon_cycle', ' Option carbon_cycle_emis_comp not yet implemented',1)
    299              END IF
    300           END IF
    301 
    302        END IF ! carbon_cycle_tr
    303    
    304        ! Reset cumluative variables
    305        IF (carbon_cycle_cpl) THEN
    306           fco2_land_day(:) = 0.
    307           fco2_lu_day(:)   = 0.
    308        END IF
    309    
    310     END IF ! newday
    311        
    312 
    313 
    314 ! -) Cumulate fluxes from ORCHIDEE at each timestep
    315     IF (carbon_cycle_cpl) THEN
    316        fco2_land_day(:) = fco2_land_day(:) + fco2_land_inst(:)
    317        fco2_lu_day(:)   = fco2_lu_day(:)   + fco2_lu_inst(:)
    318     END IF
    319 
    320 
    321 
    322 ! -)  At the end of a new day, calculate a mean scalare value of CO2 to be used by
    323 !     the radiation scheme (instead of reading value from .def)
    324 
    325 ! JG : Ici on utilise uniquement le traceur du premier couche du modele. Est-ce que c'est correcte ?
    326     IF (endday) THEN
    327        ! Calculte total area of the earth surface
    328        CALL reduce_sum(SUM(airephy),airetot)
    329        
    330 
    331        IF (carbon_cycle_tr) THEN
    332           IF (ntr_co2 == 1) THEN
    333          
    334              ! Calculate mean value of tracer CO2 to get an scalare value to be used in the
    335              ! radiation scheme (instead of reading value from .def)
    336              ! Mean value weighted with the grid cell area
    337              
    338              ! Calculate mean value
    339              fco2_tmp(:) = tr_seri(:,1,id_fco2_tot) * airephy(:)
    340              CALL reduce_sum(SUM(fco2_tmp),sumtmp)
    341              co2_ppm = sumtmp/airetot + co2_ppm0
    342              
    343           ELSE ! ntr_co2 == 4
    344              
    345              ! Calculate the delta CO2 flux
    346              tr_seri_sum(:) = tr_seri(:,1,id_fco2_fos_fuel) + tr_seri(:,1,id_fco2_land_use) + &
    347                   tr_seri(:,1,id_fco2_land) + tr_seri(:,1,id_fco2_ocn)
    348              
    349              ! Calculate mean value of delta CO2 flux
    350              fco2_tmp(:) = tr_seri_sum(:) * airephy(:)
    351              CALL reduce_sum(SUM(fco2_tmp),sumtmp)
    352              delta_co2_ppm = sumtmp/airetot
    353              
    354              ! Add initial value for co2_ppm to delta value
    355              co2_ppm = delta_co2_ppm + co2_ppm0
    356           END IF
    357 
    358        ELSE IF (carbon_cycle_cpl) THEN ! no carbon_cycle_tr
    359 
    360           ! Calculate the total CO2 flux and integrate to get scalare value for the radiation scheme
    361           fco2_tmp(:) = (fos_fuel(:) + (fco2_lu_day(:) + fco2_land_day(:))*pctsrf(:,is_ter) &
    362                + fco2_ocn_day(:)*pctsrf(:,is_oce)) * fact
    363          
    364           ! Calculate mean value
    365           fco2_tmp(:) = fco2_tmp(:) * airephy(:)
    366           CALL reduce_sum(SUM(fco2_tmp),sumtmp)
    367           delta_co2_ppm = sumtmp/airetot
    368 
    369           ! Update current value of the atmospheric co2_ppm
    370           co2_ppm = co2_ppm + delta_co2_ppm
    371          
    372        END IF ! carbon_cycle_tr
    373 
    374        ! transformation of the atmospheric CO2 concentration for the radiation code
    375        RCO2 = co2_ppm * 1.0e-06  * 44.011/28.97
    376 
    377     END IF
    378 
    379     ! Calculate CO2 flux to send to ocean and land models : PISCES and ORCHIDEE         
    380     IF (endday .AND. carbon_cycle_cpl) THEN
    381        
    382        IF (carbon_cycle_tr) THEN
    383           IF (ntr_co2==1) THEN
    384 
    385              co2_send(:) = tr_seri(:,1,id_fco2_tot) + co2_ppm0
    386 
    387           ELSE ! ntr_co2 == 4
    388 
    389              co2_send(:) = tr_seri_sum(:) + co2_ppm0
    390 
    391           END IF
     452          ! Sum all co2 tracers to get the total delta CO2 flux at first model layer
     453          fco2_tmp(:) = 0.
     454          DO it = 1, ntr_co2
     455             fco2_tmp(1:klon) = fco2_tmp(1:klon) + tr_seri(1:klon,1,co2trac(it)%id)
     456          END DO
     457          co2_send(1:klon) = fco2_tmp(1:klon) + co2_ppm0
    392458       ELSE
    393459          ! Send a scalare value in 2D variable to ocean and land model (PISCES and ORCHIDEE)
    394           co2_send(:) = co2_ppm
     460          co2_send(1:klon) = co2_ppm
    395461       END IF
    396462

  • trunk/libf/phylmd/change_srf_frac_mod.F90

    r1 r66  
    99!
    1010! Change Surface Fractions
    11 !
     11! Author J Ghattas 2008
     12
    1213  SUBROUTINE change_srf_frac(itime, dtime, jour, &
    1314       pctsrf, alb1, alb2, tsurf, u10m, v10m, pbl_tke)
    … …  
    7677    END SELECT
    7778
    78     IF (is_modified) THEN
     79
    7980!****************************************************************************************
    8081! 2)
    … …  
    8485!
    8586!****************************************************************************************
     87    IF (is_modified) THEN
    8688 
    8789! Test and exit if a fraction is negative
    … …  
    150152       CALL pbl_surface_newfrac(itime, pctsrf, pctsrf_old, tsurf, alb1, alb2, u10m, v10m, pbl_tke)
    151153
     154    ELSE
     155       ! No modifcation should be done
     156       pctsrf(:,:) = pctsrf_old(:,:)
     157
    152158    END IF ! is_modified
    153159

  • trunk/libf/phylmd/conf_phys.F90

    r1 r66  
    11
    22!
    3 ! $Id: conf_phys.F90 1423 2010-08-02 14:52:29Z jghattas $
     3! $Id: conf_phys.F90 1486 2011-02-11 12:07:39Z fairhead $
    44!
    55!
    … …  
    1313  subroutine conf_phys(ok_journe, ok_mensuel, ok_instan, ok_hf, &
    1414                       ok_LES,&
     15                       callstats,&
    1516                       solarlong0,seuil_inversion, &
    1617                       fact_cldcon, facttemps,ok_newmicro,iflag_radia,&
    … …  
    6667  logical              :: ok_journe, ok_mensuel, ok_instan, ok_hf
    6768  logical              :: ok_LES
     69  LOGICAL              :: callstats
    6870  LOGICAL              :: ok_ade, ok_aie, aerosol_couple
    6971  INTEGER              :: flag_aerosol
    … …  
    7981  logical,SAVE        :: ok_journe_omp, ok_mensuel_omp, ok_instan_omp, ok_hf_omp       
    8082  logical,SAVE        :: ok_LES_omp   
     83  LOGICAL,SAVE        :: callstats_omp
    8184  LOGICAL,SAVE        :: ok_ade_omp, ok_aie_omp, aerosol_couple_omp
    8285  INTEGER, SAVE       :: flag_aerosol_omp
    … …  
    14181421  ok_LES_omp = .false.                                             
    14191422  call getin('OK_LES', ok_LES_omp)                                 
     1423
     1424!Config Key  = callstats                                               
     1425!Config Desc = Pour des sorties callstats                                 
     1426!Config Def  = .false.                                             
     1427!Config Help = Pour creer le fichier stats contenant les sorties 
     1428!              stats                                                 
     1429!                                                                   
     1430  callstats_omp = .false.                                             
     1431  call getin('callstats', callstats_omp)                                 
    14201432!
    14211433!Config Key  = ecrit_LES
    … …  
    15811593    ok_hines = ok_hines_omp
    15821594    ok_LES = ok_LES_omp
     1595    callstats = callstats_omp
    15831596    ecrit_LES = ecrit_LES_omp
    15841597    carbon_cycle_tr = carbon_cycle_tr_omp

  • trunk/libf/phylmd/cpl_mod.F90

    r1 r66  
    100100  SUBROUTINE cpl_init(dtime, rlon, rlat)
    101101    USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, fco2_ocn_day
     102    USE surface_data
    102103
    103104    INCLUDE "dimensions.h"
    … …  
    270271    ENDIF    ! is_sequential
    271272   
     273
     274!*************************************************************************************
     275! compatibility test
     276!
     277!*************************************************************************************
     278    IF (carbon_cycle_cpl .AND. version_ocean=='opa8') THEN
     279       abort_message='carbon_cycle_cpl does not work with opa8'
     280       CALL abort_gcm(modname,abort_message,1)
     281    END IF
     282
    272283  END SUBROUTINE cpl_init
    273284 

  • trunk/libf/phylmd/cv3_routines.F

    r1 r66  
    11!
    2 ! $Id: cv3_routines.F 1403 2010-07-01 09:02:53Z fairhead $
     2! $Id: cv3_routines.F 1467 2010-12-17 11:17:22Z idelkadi $
    33!
    44c
    … …  
    8989
    9090      betad=10.0   ! original value (from convect 4.3)
     91
     92      OPEN(99,file='conv_param.data',status='old',
     93     $          form='formatted',err=9999)
     94      READ(99,*,end=9998) dpbase
     95      READ(99,*,end=9998) pbcrit
     96      READ(99,*,end=9998) ptcrit
     97      READ(99,*,end=9998) sigdz
     98      READ(99,*,end=9998) spfac
     999998  Continue
     100      CLOSE(99)
     1019999  Continue
     102      WRITE(*,*)'dpbase=',dpbase
     103      WRITE(*,*)'pbcrit=',pbcrit
     104      WRITE(*,*)'ptcrit=',ptcrit
     105      WRITE(*,*)'sigdz=',sigdz
     106      WRITE(*,*)'spfac=',spfac
    91107
    92108      return

  • trunk/libf/phylmd/phyetat0.F

    r1 r66  
    11!
    2 ! $Id: phyetat0.F 1403 2010-07-01 09:02:53Z fairhead $
     2! $Id: phyetat0.F 1457 2010-11-22 15:19:25Z musat $
    33!
    44c
    … …  
    2121      USE infotrac
    2222      USE traclmdz_mod,    ONLY : traclmdz_from_restart
    23       USE carbon_cycle_mod,ONLY : carbon_cycle_tr, carbon_cycle_cpl
     23      USE carbon_cycle_mod,ONLY :
     24     &     carbon_cycle_tr, carbon_cycle_cpl, co2_send
    2425
    2526      IMPLICIT none
    … …  
    133134
    134135       
    135 
    136          IF( clesphy0(1).NE.tab_cntrl( 5 ) )  THEN
    137              clesphy0(1)=tab_cntrl( 5 )
    138          ENDIF
    139 
    140          IF( clesphy0(2).NE.tab_cntrl( 6 ) )  THEN
    141              clesphy0(2)=tab_cntrl( 6 )
    142          ENDIF
    143 
    144          IF( clesphy0(3).NE.tab_cntrl( 7 ) )  THEN
    145              clesphy0(3)=tab_cntrl( 7 )
    146          ENDIF
    147 
    148          IF( clesphy0(4).NE.tab_cntrl( 8 ) )  THEN
    149              clesphy0(4)=tab_cntrl( 8 )
    150          ENDIF
    151 
    152          IF( clesphy0(5).NE.tab_cntrl( 9 ) )  THEN
    153              clesphy0(5)=tab_cntrl( 9 )
    154          ENDIF
    155 
    156          IF( clesphy0(6).NE.tab_cntrl( 10 ) )  THEN
    157              clesphy0(6)=tab_cntrl( 10 )
    158          ENDIF
    159 
    160          IF( clesphy0(7).NE.tab_cntrl( 11 ) )  THEN
    161              clesphy0(7)=tab_cntrl( 11 )
    162          ENDIF
    163 
    164          IF( clesphy0(8).NE.tab_cntrl( 12 ) )  THEN
    165              clesphy0(8)=tab_cntrl( 12 )
    166          ENDIF
    167 
     136      clesphy0(1)=tab_cntrl( 5 )
     137      clesphy0(2)=tab_cntrl( 6 )
     138      clesphy0(3)=tab_cntrl( 7 )
     139      clesphy0(4)=tab_cntrl( 8 )
     140      clesphy0(5)=tab_cntrl( 9 )
     141      clesphy0(6)=tab_cntrl( 10 )
     142      clesphy0(7)=tab_cntrl( 11 )
     143      clesphy0(8)=tab_cntrl( 12 )
    168144
    169145c
    … …  
    841817          ENDIF
    842818          WRITE(str2,'(i2.2)') nsrf
    843           CALL get_field("TKE"//str2,pbl_tke(:,1:klev,nsrf),found)
     819          CALL get_field("TKE"//str2,pbl_tke(:,1:klev+1,nsrf),found)
    844820          IF (.NOT. found) THEN
    845821            PRINT*, "phyetat0: <TKE"//str2//"> est absent"
    … …  
    848824          xmin = 1.0E+20
    849825          xmax = -1.0E+20
    850           DO k = 1, klev
     826          DO k = 1, klev+1
    851827            DO i = 1, klon
    852828              xmin = MIN(pbl_tke(i,k,nsrf),xmin)
    … …  
    10781054
    10791055         END DO
    1080          
    10811056         CALL traclmdz_from_restart(trs)
     1057
     1058         IF (carbon_cycle_cpl) THEN
     1059            ALLOCATE(co2_send(klon), stat=ierr)
     1060            IF (ierr /= 0) CALL abort_gcm
     1061     &           ('phyetat0','pb allocation co2_send',1)
     1062            CALL get_field("co2_send",co2_send,found)
     1063            IF (.NOT. found) THEN
     1064               PRINT*,"phyetat0: Le champ <co2_send> est absent"
     1065               PRINT*,"Initialisation uniforme a co2_ppm=",co2_ppm
     1066               co2_send(:) = co2_ppm
     1067            END IF
     1068         END IF
    10821069      END IF
    10831070

  • trunk/libf/phylmd/phyredem.F

    r1 r66  
    11!
    2 ! $Id: phyredem.F 1403 2010-07-01 09:02:53Z fairhead $
     2! $Id: phyredem.F 1457 2010-11-22 15:19:25Z musat $
    33!
    44c
    … …  
    1515      USE infotrac
    1616      USE control_mod
    17 
     17      USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, co2_send
    1818
    1919      IMPLICIT none
    … …  
    289289            WRITE(str2,'(i2.2)') nsrf
    290290            CALL put_field("TKE"//str2,"Energ. Cineti. Turb."//str2,
    291      .                     pbl_tke(:,1:klev,nsrf))
     291     .                     pbl_tke(:,1:klev+1,nsrf))
    292292          ELSE
    293293            PRINT*, "Trop de sous-mailles"
    … …  
    337337            CALL put_field("trs_"//tname(iiq),"",trs(:,it))
    338338         END DO
     339         IF (carbon_cycle_cpl) THEN
     340            IF (.NOT. ALLOCATED(co2_send)) THEN
     341               ! This is the case of create_etat0_limit, ce0l
     342               ALLOCATE(co2_send(klon))
     343               co2_send(:) = co2_ppm0
     344            END IF
     345            CALL put_field("co2_send","co2_ppm for coupling",co2_send)
     346         END IF
    339347      END IF
    340348

  • trunk/libf/phylmd/phys_output_mod.F90

    r1 r66  
    1 ! $Id: phys_output_mod.F90 1424 2010-09-01 09:20:28Z jghattas $
     1! $Id: phys_output_mod.F90 1473 2011-01-12 15:07:43Z fairhead $
    22!
    33! Abderrahmane 12 2007
    … …  
    427427  type(ctrl_out),save :: o_pres         = ctrl_out((/ 2, 3, 10, 10, 1 /),'pres')
    428428  type(ctrl_out),save :: o_paprs        = ctrl_out((/ 2, 3, 10, 10, 1 /),'paprs')
     429  type(ctrl_out),save :: o_mass        = ctrl_out((/ 2, 3, 10, 10, 1 /),'mass')
     430
    429431  type(ctrl_out),save :: o_rneb         = ctrl_out((/ 2, 5, 10, 10, 1 /),'rneb')
    430432  type(ctrl_out),save :: o_rnebcon      = ctrl_out((/ 2, 5, 10, 10, 1 /),'rnebcon')
    … …  
    11081110 CALL histdef3d(iff,o_pres%flag,o_pres%name, "Air pressure", "Pa" )
    11091111 CALL histdef3d(iff,o_paprs%flag,o_paprs%name, "Air pressure Inter-Couches", "Pa" )
     1112 CALL histdef3d(iff,o_mass%flag,o_mass%name, "Masse Couches", "kg/m2" )
    11101113 CALL histdef3d(iff,o_rneb%flag,o_rneb%name, "Cloud fraction", "-")
    11111114 CALL histdef3d(iff,o_rnebcon%flag,o_rnebcon%name, "Convective Cloud Fraction", "-")

  • trunk/libf/phylmd/phys_output_write.h

    r1 r66  
    104104     s                   o_psol%name,itau_w,zx_tmp_fi2d)
    105105       ENDIF
     106
     107       IF (o_mass%flag(iff)<=lev_files(iff)) THEN
     108      CALL histwrite_phy(nid_files(iff),o_mass%name,itau_w,zmasse)
     109        ENDIF
     110
    106111
    107112       IF (o_qsurf%flag(iff)<=lev_files(iff)) THEN

  • trunk/libf/phylmd/physiq.F

    r1 r66  
    1 ! $Id: physiq.F 1428 2010-09-13 08:43:37Z fairhead $
     1! $Id: physiq.F 1486 2011-02-11 12:07:39Z fairhead $
    22c#define IO_DEBUG
    33
    … …  
    158158      save ok_LES                           
    159159c$OMP THREADPRIVATE(ok_LES)                 
     160c
     161      LOGICAL callstats ! sortir le fichier stats
     162      save callstats                           
     163c$OMP THREADPRIVATE(callstats)                 
    160164c
    161165      LOGICAL ok_region ! sortir le fichier regional
    … …  
    606610      real, save :: ale_bl_prescr=0.
    607611
    608       real, save :: ale_max=100.
     612      real, save :: ale_max=1000.
    609613      real, save :: alp_max=2.
    610614
    … …  
    11501154!     and 360
    11511155
     1156      INTEGER ierr
    11521157#include "YOMCST.h"
    11531158#include "YOETHF.h"
    … …  
    12221227     .     ok_instan, ok_hf,
    12231228     .     ok_LES,
     1229     .     callstats,
    12241230     .     solarlong0,seuil_inversion,
    12251231     .     fact_cldcon, facttemps,ok_newmicro,iflag_radia,
    … …  
    12501256cym Attention pbase pas initialise dans concvl !!!!
    12511257          pbase=0
    1252           paire_ter(:)=0.   
    12531258cIM 180608
    12541259c         pmflxr=0.
    … …  
    18471852! doit donc etre placé avant radlwsw et pbl_surface
    18481853
     1854!!!   jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1855      call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
     1856      day_since_equinox = (jD_cur + jH_cur) - jD_eq
     1857!
     1858!   choix entre calcul de la longitude solaire vraie ou valeur fixee a
     1859!   solarlong0
     1860      if (solarlong0<-999.) then
     1861       if (new_orbit) then
    18491862! calcul selon la routine utilisee pour les planetes
    1850       if (new_orbit) then
    1851         call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
    1852         day_since_equinox = (jD_cur + jH_cur) - jD_eq
    1853 !        day_since_equinox = (jD_cur) - jD_eq
    18541863        call solarlong(day_since_equinox, zlongi, dist)
    1855       else     
     1864       else   
    18561865! calcul selon la routine utilisee pour l'AR4
    1857 !   choix entre calcul de la longitude solaire vraie ou valeur fixee a
    1858 !   solarlong0
    1859         if (solarlong0<-999.) then
    1860            CALL orbite(REAL(days_elapsed+1),zlongi,dist)
    1861         else
     1866        CALL orbite(REAL(days_elapsed+1),zlongi,dist)
     1867       endif   
     1868      else   
    18621869           zlongi=solarlong0  ! longitude solaire vraie
    1863            dist=1.            ! distance au soleil / moyenne
    1864         endif
    1865       endif
     1870           dist=1.            ! distance au soleil / moyenne
     1871      endif   
    18661872      if(prt_level.ge.1)                                                &
    18671873     &    write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
    … …  
    33763382     I                   cdragh,coefh,u1,v1,ftsol,pctsrf,
    33773383     I                   frac_impa, frac_nucl,
    3378      I                   pphis,airephy,dtime,itap)
     3384     I                   pphis,airephy,dtime,itap,
     3385     I                   rlon,rlat,qx(:,:,ivap),da,phi,mp,upwd,dnwd)
    33793386
    33803387
    … …  
    36943701c====================================================================
    36953702c
    3696      
     3703
     3704c        -----------------------------------------------------------------
     3705c        WSTATS: Saving statistics
     3706c        -----------------------------------------------------------------
     3707c        ("stats" stores and accumulates 8 key variables in file "stats.nc"
     3708c        which can later be used to make the statistic files of the run:
     3709c        "stats")          only possible in 3D runs !
     3710
     3711         
     3712         IF (callstats) THEN
     3713
     3714           call wstats(klon,o_psol%name,"Surface pressure","Pa"
     3715     &                 ,2,paprs(:,1))
     3716           call wstats(klon,o_tsol%name,"Surface temperature","K",
     3717     &                 2,zxtsol)
     3718           zx_tmp_fi2d(:) = rain_fall(:) + snow_fall(:)
     3719           call wstats(klon,o_precip%name,"Precip Totale liq+sol",
     3720     &                 "kg/(s*m2)",2,zx_tmp_fi2d)
     3721           zx_tmp_fi2d(:) = rain_lsc(:) + snow_lsc(:)
     3722           call wstats(klon,o_plul%name,"Large-scale Precip",
     3723     &                 "kg/(s*m2)",2,zx_tmp_fi2d)
     3724           zx_tmp_fi2d(:) = rain_con(:) + snow_con(:)
     3725           call wstats(klon,o_pluc%name,"Convective Precip",
     3726     &                 "kg/(s*m2)",2,zx_tmp_fi2d)
     3727           call wstats(klon,o_sols%name,"Solar rad. at surf.",
     3728     &                 "W/m2",2,solsw)
     3729           call wstats(klon,o_soll%name,"IR rad. at surf.",
     3730     &                 "W/m2",2,sollw)
     3731          zx_tmp_fi2d(:) = topsw(:)-toplw(:)
     3732          call wstats(klon,o_nettop%name,"Net dn radiatif flux at TOA",
     3733     &                 "W/m2",2,zx_tmp_fi2d)
     3734
     3735
     3736
     3737           call wstats(klon,o_temp%name,"Air temperature","K",
     3738     &                 3,t_seri)
     3739           call wstats(klon,o_vitu%name,"Zonal wind","m.s-1",
     3740     &                 3,u_seri)
     3741           call wstats(klon,o_vitv%name,"Meridional wind",
     3742     &                "m.s-1",3,v_seri)
     3743           call wstats(klon,o_vitw%name,"Vertical wind",
     3744     &                "m.s-1",3,omega)
     3745           call wstats(klon,o_ovap%name,"Specific humidity", "kg/kg",
     3746     &                 3,q_seri)
     3747 
     3748
     3749
     3750           IF(lafin) THEN
     3751             write (*,*) "Writing stats..."
     3752             call mkstats(ierr)
     3753           ENDIF
     3754
     3755         ENDIF !if callstats
     3756     
    36973757
    36983758      IF (lafin) THEN

  • trunk/libf/phylmd/phytrac.F90

    r1 r66  
    6666!--------
    6767  REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
    68   REAL,DIMENSION(klon,klev),INTENT(IN)   :: u       !
    69   REAL,DIMENSION(klon,klev),INTENT(IN)   :: v       !
     68  REAL,DIMENSION(klon,klev),INTENT(IN)   :: u       ! variable not used
     69  REAL,DIMENSION(klon,klev),INTENT(IN)   :: v       ! variable not used
    7070  REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique
    7171  REAL,DIMENSION(klon,klev),INTENT(IN)   :: rh      ! humidite relative
    … …  
    118118!--------------
    119119!
    120   REAL,DIMENSION(klon,klev),INTENT(IN) :: cdragh ! coeff drag pour T et Q
     120  REAL,DIMENSION(klon),INTENT(IN)      :: cdragh ! coeff drag pour T et Q
    121121  REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh  ! coeff melange CL (m**2/s)
    122122  REAL,DIMENSION(klon),INTENT(IN)      :: yu1    ! vents au premier niveau
    … …  
    213213     SELECT CASE(type_trac)
    214214     CASE('lmdz')
    215 !IM ajout t_seri, pplay, sh    CALL traclmdz_init(pctsrf, ftsol, tr_seri, aerosol, lessivage)
    216         CALL traclmdz_init(pctsrf, ftsol, tr_seri, t_seri, pplay, sh, aerosol, lessivage)
     215        CALL traclmdz_init(pctsrf, ftsol, tr_seri, t_seri, pplay, sh, pdtphys, aerosol, lessivage)
    217216     CASE('inca')
    218217        source(:,:)=0.

  • trunk/libf/phylmd/read_map2D.F90

    r1 r66  
    1111
    1212! Input arguments
    13   CHARACTER(len=20), INTENT(IN) :: filename     ! name of file to read
    14   CHARACTER(len=20), INTENT(IN) :: varname      ! name of variable in file
     13  CHARACTER(len=*), INTENT(IN) :: filename     ! name of file to read
     14  CHARACTER(len=*), INTENT(IN) :: varname      ! name of variable in file
    1515  INTEGER, INTENT(IN)           :: timestep     ! actual timestep
    1616  LOGICAL, INTENT(IN)           :: inverse      ! TRUE if latitude needs to be inversed
    … …  
    2727  REAL, DIMENSION(nbp_lon,nbp_lat) :: var_glo2D_tmp ! 2D global
    2828  REAL, DIMENSION(klon_glo)        :: var_glo1D     ! 1D global
    29 
     29  INCLUDE "iniprint.h"
    3030
    3131! Read variable from file. Done by master process MPI and master thread OpenMP
    3232  IF (is_mpi_root .AND. is_omp_root) THEN
    33      ierr = NF90_OPEN (filename, NF90_NOWRITE, nid)
    34      IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Problem in opening file '//filename,1)
     33     ierr = NF90_OPEN(trim(filename), NF90_NOWRITE, nid)
     34     IF (ierr /= NF90_NOERR) CALL write_err_mess('Problem in opening file')
    3535
    36      ierr = NF90_INQ_VARID(nid, varname, nvarid)
    37      IF (ierr /= NF90_NOERR) CALL abort_gcm(modname, 'The variable '//varname//' is absent in file',1)
     36     ierr = NF90_INQ_VARID(nid, trim(varname), nvarid)
     37     IF (ierr /= NF90_NOERR) CALL write_err_mess('The variable is absent in file')
    3838     
    3939     start=(/1,1,timestep/)
    4040     count=(/nbp_lon,nbp_lat,1/)
    4141     ierr = NF90_GET_VAR(nid, nvarid, var_glo2D,start,count)
    42      IF (ierr /= NF90_NOERR) CALL abort_gcm(modname, 'Problem in reading varaiable '//varname,1)
     42     IF (ierr /= NF90_NOERR) CALL write_err_mess('Problem in reading varaiable')
     43
     44     ierr = NF90_CLOSE(nid)
     45     IF (ierr /= NF90_NOERR) CALL write_err_mess('Problem in closing file')
    4346
    4447     ! Inverse latitude order
    … …  
    5356     CALL grid2Dto1D_glo(var_glo2D,var_glo1D)
    5457
     58     WRITE(lunout,*) 'in read_map2D, filename = ', trim(filename)
     59     WRITE(lunout,*) 'in read_map2D, varname  = ', trim(varname)
     60     WRITE(lunout,*) 'in read_map2D, timestep = ', timestep
    5561  ENDIF
    5662
    … …  
    5864  CALL scatter(var_glo1D, varout)
    5965
     66  CONTAINS
     67    SUBROUTINE write_err_mess(err_mess)
     68
     69      CHARACTER(len=*), INTENT(IN) :: err_mess
     70      INCLUDE "iniprint.h"
     71     
     72      WRITE(lunout,*) 'Error in read_map2D, filename = ', trim(filename)
     73      WRITE(lunout,*) 'Error in read_map2D, varname  = ', trim(varname)
     74      WRITE(lunout,*) 'Error in read_map2D, timestep = ', timestep
     75
     76      CALL abort_gcm(modname, err_mess, 1)
     77
     78    END SUBROUTINE write_err_mess
     79
    6080END SUBROUTINE read_map2D

  • trunk/libf/phylmd/readaerosol.F90

    r1 r66  
    1 ! $Id: readaerosol.F90 1403 2010-07-01 09:02:53Z fairhead $
     1! $Id: readaerosol.F90 1484 2011-02-09 15:44:57Z jghattas $
    22!
    33MODULE readaerosol_mod
    … …  
    77CONTAINS
    88
    9 SUBROUTINE readaerosol(name_aero, type, iyr_in, klev_src, pt_ap, pt_b, pt_out, psurf, load)
     9SUBROUTINE readaerosol(name_aero, type, filename, iyr_in, klev_src, pt_ap, pt_b, pt_out, psurf, load)
    1010
    1111!****************************************************************************************
    … …  
    2727  ! Input arguments
    2828  CHARACTER(len=7), INTENT(IN) :: name_aero
    29   CHARACTER(len=*), INTENT(IN) :: type  ! correspond to aer_type in clesphys.h
     29  CHARACTER(len=*), INTENT(IN) :: type  ! actuel, annuel, scenario or preind
     30  CHARACTER(len=8), INTENT(IN) :: filename
    3031  INTEGER, INTENT(IN)          :: iyr_in
    3132
    … …  
    5859     ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
    5960     ! pt_out has dimensions (klon, klev_src, 12)
    60      CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
     61     CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
    6162     
    6263
    … …  
    6768     ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
    6869     ! pt_out has dimensions (klon, klev_src, 12)
    69      CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
     70     CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
    7071     
    7172  ELSE IF (type == 'annuel') THEN
    … …  
    7677     ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month
    7778     ! pt_out has dimensions (klon, klev_src, 12)
    78      CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
     79     CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
    7980     
    8081  ELSE IF (type == 'scenario') THEN
    … …  
    8687        ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
    8788        ! pt_out has dimensions (klon, klev_src, 12)
    88         CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
     89        CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
    8990       
    9091     ELSE IF (iyr_in .GE. 2100) THEN
    … …  
    9394        ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
    9495        ! pt_out has dimensions (klon, klev_src, 12)
    95         CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
     96        CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
    9697       
    9798     ELSE
    … …  
    113114        ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
    114115        ! pt_out has dimensions (klon, klev_src, 12)
    115         CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
     116        CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
    116117       
    117118        ! If to read two decades:
    … …  
    125126           ! get_aero_fromfile returns pt_2 allocated and initialized with data for 12 month
    126127           ! pt_2 has dimensions (klon, klev_src, 12)
    127            CALL get_aero_fromfile(name_aero, cyear, klev_src2, pt_ap, pt_b, p0, pt_2, psurf2, load2)
     128           CALL get_aero_fromfile(name_aero, cyear, filename, klev_src2, pt_ap, pt_b, p0, pt_2, psurf2, load2)
    128129           ! Test for same number of vertical levels
    129130           IF (klev_src /= klev_src2) THEN
    … …  
    160161
    161162  ELSE
    162      WRITE(lunout,*)'This option is not implemented : aer_type = ', type
     163     WRITE(lunout,*)'This option is not implemented : aer_type = ', type,' name_aero=',name_aero
    163164     CALL abort_gcm('readaerosol','Error : aer_type parameter not accepted',1)
    164165  END IF ! type
    … …  
    168169
    169170
    170   SUBROUTINE get_aero_fromfile(varname, cyr, klev_src, pt_ap, pt_b, p0, pt_year, psurf_out, load_out)
     171  SUBROUTINE get_aero_fromfile(varname, cyr, filename, klev_src, pt_ap, pt_b, p0, pt_year, psurf_out, load_out)
    171172!****************************************************************************************
    172173! Read 12 month aerosol from file and distribute to local process on physical grid.
    … …  
    200201    CHARACTER(len=7), INTENT(IN)          :: varname
    201202    CHARACTER(len=4), INTENT(IN)          :: cyr
     203    CHARACTER(len=8), INTENT(IN)          :: filename
    202204
    203205! Output arguments
    … …  
    213215! Local variables
    214216    CHARACTER(len=30)     :: fname
    215     CHARACTER(len=8)      :: filename='aerosols'
    216217    CHARACTER(len=30)     :: cvar
    217218    INTEGER               :: ncid, dimid, varid
    … …  
    242243! 1) Open file
    243244!****************************************************************************************
    244        fname = filename//cyr//'.nc'
     245! Add suffix to filename
     246       fname = trim(filename)//cyr//'.nc'
    245247 
    246        WRITE(lunout,*) 'reading ', TRIM(fname)
     248       WRITE(lunout,*) 'reading variable ',TRIM(varname),' in file ', TRIM(fname)
    247249       CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, ncid) )
    248250
    … …  
    283285          CALL abort_gcm('get_aero_fromfile', 'latitudes do not correspond between file and model',1)
    284286       END IF
    285 
    286 ! 1.5) Check number of month in file opened
    287 !
    288 !**************************************************************************************************
    289        ierr = nf90_inq_dimid(ncid, 'TIME',dimid)
    290        CALL check_err( nf90_inquire_dimension(ncid, dimid, len = nbr_tsteps) )
    291 !       IF (nbr_tsteps /= 12 .AND. nbr_tsteps /= 14) THEN
    292        IF (nbr_tsteps /= 12 ) THEN
    293          CALL abort_gcm('get_aero_fromfile', 'not the right number of months in aerosol file read (should be 12 for the moment)',1)
    294        ENDIF
    295287
    296288
    … …  
    335327
    336328       IF (new_file) THEN
     329! ++) Check number of month in file opened
     330!**************************************************************************************************
     331       ierr = nf90_inq_dimid(ncid, 'TIME',dimid)
     332       CALL check_err( nf90_inquire_dimension(ncid, dimid, len = nbr_tsteps) )
     333!       IF (nbr_tsteps /= 12 .AND. nbr_tsteps /= 14) THEN
     334       IF (nbr_tsteps /= 12 ) THEN
     335         CALL abort_gcm('get_aero_fromfile', 'not the right number of months in aerosol file read (should be 12 for the moment)',1)
     336       ENDIF
    337337
    338338! ++) Read the aerosol concentration month by month and concatenate to total variable varyear

  • trunk/libf/phylmd/readaerosol_interp.F90

    r1 r66  
    9292  LOGICAL,SAVE       :: debug=.FALSE.! Debugging in this subroutine
    9393!$OMP THREADPRIVATE(vert_interp, debug)
     94  CHARACTER(len=8)      :: type
     95  CHARACTER(len=8)      :: filename
    9496
    9597
    … …  
    173175     ! Reading values corresponding to the closest year taking into count the choice of aer_type.
    174176     ! For aer_type=scenario interpolation between 2 data sets is done in readaerosol.
    175      CALL readaerosol(name_aero(id_aero), aer_type, iyr, klev_src, pt_ap, pt_b, pt_tmp, &
     177     ! If aer_type=mix1 or mix2, the run type and file name depends on the aerosol.
     178     IF (aer_type=='preind' .OR. aer_type=='actuel' .OR. aer_type=='annuel' .OR. aer_type=='scenario') THEN
     179        ! Standard case
     180        filename='aerosols'
     181        type=aer_type
     182     ELSE IF (aer_type == 'mix1') THEN
     183        ! Special case using a mix of decenal sulfate file and annual aerosols(all aerosols except sulfate)
     184        IF (name_aero(id_aero) == 'SO4') THEN
     185           filename='so4.run '
     186           type='scenario'
     187        ELSE
     188           filename='aerosols'
     189           type='annuel'
     190        END IF
     191     ELSE  IF (aer_type == 'mix2') THEN
     192        ! Special case using a mix of decenal sulfate file and natrual aerosols
     193        IF (name_aero(id_aero) == 'SO4') THEN
     194           filename='so4.run '
     195           type='scenario'
     196        ELSE
     197           filename='aerosols'
     198           type='preind'
     199        END IF
     200     ELSE
     201        CALL abort_gcm('readaerosol_interp', 'this aer_type not supported',1)
     202     END IF
     203
     204     CALL readaerosol(name_aero(id_aero), type, filename, iyr, klev_src, pt_ap, pt_b, pt_tmp, &
    176205          psurf_year(:,:,id_aero), load_year(:,:,id_aero))
    177206     IF (.NOT. ALLOCATED(var_year)) THEN
    … …  
    182211
    183212     ! Reading values corresponding to the preindustrial concentrations.
    184      CALL readaerosol(name_aero(id_aero), 'preind', iyr, pi_klev_src, pt_ap, pt_b, pt_tmp, &
     213     type='preind'
     214     CALL readaerosol(name_aero(id_aero), type, filename, iyr, pi_klev_src, pt_ap, pt_b, pt_tmp, &
    185215          pi_psurf_year(:,:,id_aero), pi_load_year(:,:,id_aero))
    186216

  • trunk/libf/phylmd/surf_land_orchidee_mod.F90

    r1 r66  
    4343    USE mod_surf_para
    4444    USE mod_synchro_omp
    45    
    46 USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, fco2_land_inst, fco2_lu_inst
     45    USE carbon_cycle_mod, ONLY : carbon_cycle_cpl
    4746
    4847!   
    … …  
    138137    INTEGER                                   :: error
    139138    REAL, DIMENSION(klon)                     :: swdown_vrai
    140     REAL, DIMENSION(klon)                     :: fco2_land_comp  ! sur grille compresse
    141     REAL, DIMENSION(klon)                     :: fco2_lu_comp    ! sur grille compresse
    142139    CHARACTER (len = 20)                      :: modname = 'surf_land_orchidee'
    143140    CHARACTER (len = 80)                      :: abort_message
    … …  
    348345       ENDIF
    349346!
    350 ! Allocate variables needed for carbon_cycle_mod
     347! carbon_cycle_cpl not possible with this interface and version of ORHCHIDEE
    351348!
    352349       IF (carbon_cycle_cpl) THEN
    353           IF (.NOT. ALLOCATED(fco2_land_inst)) THEN
    354              ALLOCATE(fco2_land_inst(klon),stat=error)
    355              IF (error /= 0)  CALL abort_gcm(modname,'Pb in allocation fco2_land_inst',1)
    356              
    357              ALLOCATE(fco2_lu_inst(klon),stat=error)
    358              IF(error /=0) CALL abort_gcm(modname,'Pb in allocation fco2_lu_inst',1)
    359           END IF
     350          abort_message='carbon_cycle_cpl not yet possible with this interface of ORCHIDEE'
     351          CALL abort_gcm(modname,abort_message,1)
    360352       END IF
    361353       
    … …  
    464456
    465457   
    466 ! JG : TEMPORAIRE!!!! Les variables fco2_land_comp et fco2_lu_comp seront plus tard en sortie d'ORCHIDEE
    467 !      ici mis a valeur quelquonque pour test. Ces variables sont sur la grille compresse avec uniquement des points de terres
    468 
    469     fco2_land_comp(:) = 1.
    470     fco2_lu_comp(:)   = 10.
    471 
    472 ! Decompress variables for the module carbon_cycle_mod
    473     IF (carbon_cycle_cpl) THEN
    474        fco2_land_inst(:)=0.
    475        fco2_lu_inst(:)=0.
    476        
    477        DO igrid = 1, knon
    478           ireal = knindex(igrid)
    479           fco2_land_inst(ireal) = fco2_land_comp(igrid)
    480           fco2_lu_inst(ireal)   = fco2_lu_comp(igrid)
    481        END DO
    482     END IF
    483 
    484458  END SUBROUTINE surf_land_orchidee
    485459!

  • trunk/libf/phylmd/surf_land_orchidee_noopenmp_mod.F90

    r1 r66  
    138138    INTEGER                                   :: ij, jj, igrid, ireal, index
    139139    INTEGER                                   :: error
     140    INTEGER, SAVE                             :: nb_fields_cpl ! number of fields for the climate-carbon coupling (between ATM and ORCHIDEE).
     141    REAL, SAVE, ALLOCATABLE, DIMENSION(:,:)   :: fields_cpl    ! Fluxes for the climate-carbon coupling
    140142    REAL, DIMENSION(klon)                     :: swdown_vrai
    141     REAL, DIMENSION(klon)                     :: fco2_land_comp  ! sur grille compresse
    142     REAL, DIMENSION(klon)                     :: fco2_lu_comp    ! sur grille compresse
    143143    CHARACTER (len = 20)                      :: modname = 'surf_land_orchidee'
    144144    CHARACTER (len = 80)                      :: abort_message
    … …  
    210210 
    211211    IF (debut) THEN
     212! Test de coherence
     213#ifndef ORCH_NEW
     214       ! Compilation avec orchidee nouvelle version necessaire avec carbon_cycle_cpl=y
     215       IF (carbon_cycle_cpl) THEN
     216          abort_message='You must define preprossing key ORCH_NEW when running carbon_cycle_cpl=y'
     217          CALL abort_gcm(modname,abort_message,1)
     218       END IF
     219#endif
    212220       ALLOCATE(ktindex(knon))
    213221       IF ( .NOT. ALLOCATED(albedo_keep)) THEN
    … …  
    342350!
    343351! Allocate variables needed for carbon_cycle_mod
    344 !
     352       IF ( carbon_cycle_cpl ) THEN
     353          nb_fields_cpl=2
     354       ELSE
     355          nb_fields_cpl=1
     356       END IF
     357
     358
    345359       IF (carbon_cycle_cpl) THEN
    346           IF (.NOT. ALLOCATED(fco2_land_inst)) THEN
    347              ALLOCATE(fco2_land_inst(klon),stat=error)
    348              IF (error /= 0)  CALL abort_gcm(modname,'Pb in allocation fco2_land_inst',1)
    349              
    350              ALLOCATE(fco2_lu_inst(klon),stat=error)
    351              IF(error /=0) CALL abort_gcm(modname,'Pb in allocation fco2_lu_inst',1)
    352           END IF
     360          ALLOCATE(fco2_land_inst(klon),stat=error)
     361          IF (error /= 0)  CALL abort_gcm(modname,'Pb in allocation fco2_land_inst',1)
     362         
     363          ALLOCATE(fco2_lu_inst(klon),stat=error)
     364          IF(error /=0) CALL abort_gcm(modname,'Pb in allocation fco2_lu_inst',1)
    353365       END IF
     366
     367       ALLOCATE(fields_cpl(klon,nb_fields_cpl), stat = error)
     368       IF (error /= 0) CALL abort_gcm(modname,'Pb in allocation fields_cpl',1)
    354369
    355370    ENDIF                          ! (fin debut)
    … …  
    406421               evap, fluxsens, fluxlat, coastalflow, riverflow, &
    407422               tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, &
    408                lon_scat, lat_scat, q2m, t2m)
     423               lon_scat, lat_scat, q2m, t2m &
     424#ifdef ORCH_NEW
     425               , nb_fields_cpl, fields_cpl)
     426#else
     427               )
     428#endif
    409429
    410430#else         
    411           ! Interface for ORCHIDEE version 1.9 or later(1.9.2, 1.9.3, 1.9.4) compiled in parallel mode(with preprocessing flag CPP_MPI)
     431          ! Interface for ORCHIDEE version 1.9 or later(1.9.2, 1.9.3, 1.9.4, 1.9.5) compiled in parallel mode(with preprocessing flag CPP_MPI)
    412432          CALL intersurf_main (itime+itau_phy-1, iim, jjm+1, offset, knon, ktindex, &
    413433               orch_comm, dtime, lrestart_read, lrestart_write, lalo, &
    … …  
    418438               evap(1:knon), fluxsens(1:knon), fluxlat(1:knon), coastalflow(1:knon), riverflow(1:knon), &
    419439               tsol_rad(1:knon), tsurf_new(1:knon), qsurf(1:knon), albedo_out(1:knon,:), emis_new(1:knon), z0_new(1:knon), &
    420                lon_scat, lat_scat, q2m, t2m)
     440               lon_scat, lat_scat, q2m, t2m &
     441#ifdef ORCH_NEW
     442               , nb_fields_cpl, fields_cpl(1:knon,:))
     443#else
     444               )
     445#endif
    421446#endif
    422447         
    … …  
    431456
    432457    IF (knon /=0) THEN
    433    
    434458#ifndef CPP_MPI
    435459       ! Interface for ORCHIDEE compiled in sequential mode(without preprocessing flag CPP_MPI)
    … …  
    442466            evap, fluxsens, fluxlat, coastalflow, riverflow, &
    443467            tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, &
    444             lon_scat, lat_scat, q2m, t2m)
    445        
     468            lon_scat, lat_scat, q2m, t2m &
     469#ifdef ORCH_NEW
     470            , nb_fields_cpl, fields_cpl)
     471#else
     472            )
     473#endif
    446474#else
    447475       ! Interface for ORCHIDEE version 1.9 or later compiled in parallel mode(with preprocessing flag CPP_MPI)
    … …  
    454482            evap(1:knon), fluxsens(1:knon), fluxlat(1:knon), coastalflow(1:knon), riverflow(1:knon), &
    455483            tsol_rad(1:knon), tsurf_new(1:knon), qsurf(1:knon), albedo_out(1:knon,:), emis_new(1:knon), z0_new(1:knon), &
    456             lon_scat, lat_scat, q2m, t2m)
    457 #endif
    458        
     484            lon_scat, lat_scat, q2m, t2m &
     485#ifdef ORCH_NEW
     486            , nb_fields_cpl, fields_cpl(1:knon,:))
     487#else
     488            )
     489#endif
     490#endif
    459491    ENDIF
    460492
    … …  
    478510
    479511    IF (debut) lrestart_read = .FALSE.
    480 
    481 
    482 ! JG : TEMPORAIRE!!!! Les variables fco2_land_comp et fco2_lu_comp seront plus tard en sortie d'ORCHIDEE
    483 !      ici mis a valeur quelquonque pour test. Ces variables sont sur la grille compresse avec uniquement des points de terres
    484 
    485     fco2_land_comp(:) = 1.
    486     fco2_lu_comp(:)   = 10.
    487512
    488513! Decompress variables for the module carbon_cycle_mod
    … …  
    493518       DO igrid = 1, knon
    494519          ireal = knindex(igrid)
    495           fco2_land_inst(ireal) = fco2_land_comp(igrid)
    496           fco2_lu_inst(ireal)   = fco2_lu_comp(igrid)
     520          fco2_land_inst(ireal) = fields_cpl(igrid,1)
     521          fco2_lu_inst(ireal)   = fields_cpl(igrid,2)
    497522       END DO
    498523    END IF

  • trunk/libf/phylmd/traclmdz_mod.F90

    r1 r66  
    8484
    8585
    86   SUBROUTINE traclmdz_init(pctsrf, ftsol, tr_seri, t_seri, pplay, sh, aerosol, lessivage)
     86  SUBROUTINE traclmdz_init(pctsrf, ftsol, tr_seri, t_seri, pplay, sh, pdtphys, aerosol, lessivage)
    8787    ! This subroutine allocates and initialize module variables and control variables.
    8888    ! Initialization of the tracers should be done here only for those not found in the restart file.
    … …  
    104104    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
    105105    REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique
     106    REAL,INTENT(IN)                        :: pdtphys ! Pas d'integration pour la physique (seconde) 
    106107
    107108! Output variables
    … …  
    226227! ----------------------------------------------
    227228    IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
    228        CALL carbon_cycle_init(tr_seri, aerosol, radio)
     229       CALL carbon_cycle_init(tr_seri, pdtphys, aerosol, radio)
    229230    END IF
    230231
    … …  
    312313!--------------
    313314!
    314     REAL,DIMENSION(klon,klev),INTENT(IN) :: cdragh     ! coeff drag pour T et Q
     315    REAL,DIMENSION(klon),INTENT(IN)      :: cdragh     ! coeff drag pour T et Q
    315316    REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh      ! coeff melange CL (m**2/s)
    316317    REAL,DIMENSION(klon),INTENT(IN)      :: yu1        ! vents au premier niveau
    … …  
    546547!======================================================================
    547548    IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
    548        CALL carbon_cycle(nstep, pdtphys, pctsrf, tr_seri)
     549       CALL carbon_cycle(nstep, pdtphys, pctsrf, tr_seri, source)
    549550    END IF
    550551
Note: See TracChangeset for help on using the changeset viewer.