Changeset 1179


Ignore:
Timestamp:
Jun 11, 2009, 4:18:47 PM (15 years ago)
Author:
jghattas
Message:
  • Ajout de l'interpolation vertical pour les nouveaux fichiers de forcage des aerosols. Utilisant les anciennes fichiers de SO4 pas d'interpolation possible. Convergence numerique avec la version precedente en utilisant les anciens fichiers des SO4. aerosol_optic.F90 change du nom pour readaerosol_optic.F90 (lecture d'aerosol + optic) Les fichiers de forcage aerosol doit maintenant avoir le suffix .nc.
  • Correction des bugs pour inca et certain diagnostiques optionelles de radlwsw.
  • Ajout de test pour le choix advection schema.
Location:
LMDZ4/branches/LMDZ4-dev/libf
Files:
1 added
10 edited
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ4/branches/LMDZ4-dev/libf/dyn3d/infotrac.F90

    r1114 r1179  
     1! $Id$
     2!
    13MODULE infotrac
    24
    … …  
    299301    END DO
    300302
     303!
     304! Test for advection schema.
     305! This version of LMDZ only garantees iadv=10 and iadv=14 (14 only for water vapour) .
     306!
     307    DO iq=1,nqtot
     308       IF (iadv(iq)/=10 .AND. iadv(iq)/=14) THEN
     309          WRITE(lunout,*)'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
     310          CALL abort_gcm('infotrac_init','In this version only iadv=10 and iadv=14 is tested!',1)
     311       ELSE IF (iadv(iq)==14 .AND. iq/=1) THEN
     312          WRITE(lunout,*)'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
     313          CALL abort_gcm('infotrac_init','In this version iadv=14 is only permitted for water vapour!',1)
     314       END IF
     315    END DO
     316
    301317!-----------------------------------------------------------------------
    302318! Finalize :

  • LMDZ4/branches/LMDZ4-dev/libf/dyn3d/pres2lev.F

    r1046 r1179  
    1 !
    2 ! $Header$
     1! $Id$
    32!
    43c******************************************************
    … …  
    2120c  ARGUMENTS
    2221c  """""""""
    23        LOGICAL ok_invertp
    24        INTEGER lmo ! dimensions ancienne couches (input)
    25        INTEGER lmn ! dimensions nouvelle couches (input)
    26        INTEGER lmomx ! dimensions ancienne couches (input)
    27        INTEGER lmnmx ! dimensions nouvelle couches (input)
     22       LOGICAL, INTENT(IN) :: ok_invertp
     23       INTEGER, INTENT(IN) :: lmo ! dimensions ancienne couches
     24       INTEGER, INTENT(IN) :: lmn ! dimensions nouvelle couches
     25       INTEGER lmomx ! dimensions ancienne couches
     26       INTEGER lmnmx ! dimensions nouvelle couches
    2827
    2928       parameter(lmomx=10000,lmnmx=10000)
    3029
    31         real po(ni,nj,lmo)! niveau de pression ancienne grille
    32         real pn(ni,nj,lmn) ! niveau de pression nouvelle grille
     30        real, INTENT(IN) :: po(ni,nj,lmo) ! niveau de pression ancienne grille
     31        real, INTENT(IN) :: pn(ni,nj,lmn) ! niveau de pression nouvelle grille
    3332
    34        INTEGER i,j,Nhoriz,ni,nj ! nombre de point horizontale (input)
     33       INTEGER, INTENT(IN) :: ni,nj ! nombre de point horizontale
    3534
    36        REAL varo(ni,nj,lmo) ! var dans l'ancienne grille (input)
    37        REAL varn(ni,nj,lmn) ! var dans la nouvelle grille (output)
     35       REAL, INTENT(IN)  :: varo(ni,nj,lmo) ! var dans l'ancienne grille
     36       REAL, INTENT(OUT) :: varn(ni,nj,lmn) ! var dans la nouvelle grille
    3837
    3938       real zvaro(lmomx),zpo(lmomx)
    … …  
    4140c Autres variables
    4241c """"""""""""""""
    43        INTEGER n, ln ,lo 
     42       INTEGER n, ln ,lo, i, j, Nhoriz
    4443       REAL coef
    4544

  • LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/infotrac.F90

    r1117 r1179  
     1! $Id$
     2!
    13MODULE infotrac
    24
    35! nqtot : total number of tracers and higher order of moment, water vapor and liquid included
    46  INTEGER, SAVE :: nqtot
    5 !!$OMP THREADPRIVATE(nqtot)   
    67
    78! nbtr : number of tracers not including higher order of moment or water vapor or liquid
    89!        number of tracers used in the physics
    910  INTEGER, SAVE :: nbtr
    10 !!$OMP THREADPRIVATE(nbtr)   
    1111
    1212! Name variables
    1313  CHARACTER(len=20), ALLOCATABLE, DIMENSION(:), SAVE :: tname ! tracer short name for restart and diagnostics
    1414  CHARACTER(len=23), ALLOCATABLE, DIMENSION(:), SAVE :: ttext ! tracer long name for diagnostics
    15 !!$OMP THREADPRIVATE(tname,ttext)   
    1615
    1716! iadv  : index of trasport schema for each tracer
    1817  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: iadv
    19 !!$OMP THREADPRIVATE(iadv)   
    2018
    2119! niadv : vector keeping the coorspondance between all tracers(nqtot) treated in the
    2220!         dynamic part of the code and the tracers (nbtr+2) used in the physics part of the code.
    2321  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: niadv ! equivalent dyn / physique
    24 !!$OMP THREADPRIVATE(niadv)   
    2522
    2623! Variables for INCA
    2724  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE  :: conv_flg
    2825  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE  :: pbl_flg
    29 !!$OMP THREADPRIVATE(conv_flg, pbl_flg)   
    3026
    3127CONTAINS
    … …  
    305301    END DO
    306302
     303!
     304! Test for advection schema.
     305! This version of LMDZ only garantees iadv=10 and iadv=14 (14 only for water vapour) .
     306!
     307    DO iq=1,nqtot
     308       IF (iadv(iq)/=10 .AND. iadv(iq)/=14) THEN
     309          WRITE(lunout,*)'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
     310          CALL abort_gcm('infotrac_init','In this version only iadv=10 and iadv=14 is tested!',1)
     311       ELSE IF (iadv(iq)==14 .AND. iq/=1) THEN
     312          WRITE(lunout,*)'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
     313          CALL abort_gcm('infotrac_init','In this version iadv=14 is only permitted for water vapour!',1)
     314       END IF
     315    END DO
     316
    307317!-----------------------------------------------------------------------
    308318! Finalize :

  • LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/pres2lev.F

    r1046 r1179  
    1 !
    2 ! $Header$
     1! $Id$
    32!
    43c******************************************************
    … …  
    2120c  ARGUMENTS
    2221c  """""""""
    23        LOGICAL ok_invertp
    24        INTEGER lmo ! dimensions ancienne couches (input)
    25        INTEGER lmn ! dimensions nouvelle couches (input)
    26        INTEGER lmomx ! dimensions ancienne couches (input)
    27        INTEGER lmnmx ! dimensions nouvelle couches (input)
     22       LOGICAL, INTENT(IN) :: ok_invertp
     23       INTEGER, INTENT(IN) :: lmo ! dimensions ancienne couches
     24       INTEGER, INTENT(IN) :: lmn ! dimensions nouvelle couches
     25       INTEGER lmomx ! dimensions ancienne couches
     26       INTEGER lmnmx ! dimensions nouvelle couches
    2827
    2928       parameter(lmomx=10000,lmnmx=10000)
    3029
    31         real po(ni,nj,lmo)! niveau de pression ancienne grille
    32         real pn(ni,nj,lmn) ! niveau de pression nouvelle grille
     30        real, INTENT(IN) :: po(ni,nj,lmo) ! niveau de pression ancienne grille
     31        real, INTENT(IN) :: pn(ni,nj,lmn) ! niveau de pression nouvelle grille
    3332
    34        INTEGER i,j,Nhoriz,ni,nj ! nombre de point horizontale (input)
     33       INTEGER, INTENT(IN) :: ni,nj ! nombre de point horizontale
    3534
    36        REAL varo(ni,nj,lmo) ! var dans l'ancienne grille (input)
    37        REAL varn(ni,nj,lmn) ! var dans la nouvelle grille (output)
     35       REAL, INTENT(IN)  :: varo(ni,nj,lmo) ! var dans l'ancienne grille
     36       REAL, INTENT(OUT) :: varn(ni,nj,lmn) ! var dans la nouvelle grille
    3837
    3938       real zvaro(lmomx),zpo(lmomx)
    … …  
    4140c Autres variables
    4241c """"""""""""""""
    43        INTEGER n, ln ,lo 
     42       INTEGER n, ln ,lo, i, j, Nhoriz
    4443       REAL coef
    4544

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/newmicro.F

    r1092 r1179  
    1 !
    2 ! $Header$
    3 !
     1! $Id$
     2!     
    43      SUBROUTINE newmicro (paprs, pplay,ok_newmicro,
    54     .                  t, pqlwp, pclc, pcltau, pclemi,
    … …  
    76     s                  xflwp, xfiwp, xflwc, xfiwc,
    87     e                  ok_aie,
    9      e                  sulfate, sulfate_pi,
     8     e                  mass_ins_aero, mass_ins_aero_pi,
    109     e                  bl95_b0, bl95_b1,
    1110     s                  cldtaupi, re, fl)
    … …  
    2221c
    2322c ok_aie--input-L-apply aerosol indirect effect or not
    24 c sulfate-input-R-sulfate aerosol mass concentration [um/m^3]
    25 c sulfate_pi-input-R-dito, pre-industrial value
     23c mass_ins_aero-----input-R-total mass concentration for all indissoluble aerosols[ug/m^3]
     24c mass_ins_aero_pi--input-R-dito, pre-industrial value
    2625c bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
    2726c bl95_b1-input-R-a parameter, may be varied for tests (    -"-      )
    … …  
    9493      LOGICAL ok_a1lwpdep       ! a1 LWP dependent?
    9594     
    96       REAL sulfate(klon, klev)  ! sulfate aerosol mass concentration [ug m-3]
     95      REAL mass_ins_aero(klon, klev)    ! total mass concentration for all indissoluble aerosols [ug m-3]
     96      REAL mass_ins_aero_pi(klon, klev) ! - " - (pre-industrial value)
    9797      REAL cdnc(klon, klev)     ! cloud droplet number concentration [m-3]
    9898      REAL re(klon, klev)       ! cloud droplet effective radius [um]
    99       REAL sulfate_pi(klon, klev)  ! sulfate aerosol mass concentration [ug m-3] (pre-industrial value)
    10099      REAL cdnc_pi(klon, klev)     ! cloud droplet number concentration [m-3] (pi value)
    101100      REAL re_pi(klon, klev)       ! cloud droplet effective radius [um] (pi value)
    … …  
    157156                                !             
    158157                  cdnc(i,k) = 10.**(bl95_b0+bl95_b1*
    159      &                 log(MAX(sulfate(i,k),1.e-4))/log(10.))*1.e6 !-m-3
     158     &                 log(MAX(mass_ins_aero(i,k),1.e-4))/log(10.))*1.e6 !-m-3
    160159                                ! Cloud droplet number concentration (CDNC) is restricted
    161160                                ! to be within [20, 1000 cm^3]
    … …  
    165164                                !
    166165                  cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1*
    167      &                 log(MAX(sulfate_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3
     166     &                 log(MAX(mass_ins_aero_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3
    168167                  cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k)))
    169168               ENDDO

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/nuage.F

    r766 r1179  
    1 !
    2 ! $Header$
     1! $Id$
    32!
    43      SUBROUTINE nuage (paprs, pplay,
    … …  
    65     .                  pch, pcl, pcm, pct, pctlwp,
    76     e                  ok_aie,
    8      e                  sulfate, sulfate_pi,
     7     e                  mass_ins_aero, mass_ins_aero_pi,
    98     e                  bl95_b0, bl95_b1,
    109     s                  cldtaupi, re, fl)
    … …  
    2019c pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
    2120c ok_aie--input-L-apply aerosol indirect effect or not
    22 c sulfate-input-R-sulfate aerosol mass concentration [um/m^3]
    23 c sulfate_pi-input-R-dito, pre-industrial value
     21c mass_ins_aero-----input-R-total mass concentration for all indissoluble aerosols[ug/m^3]
     22c mass_ins_aero_pi--input-R-dito, pre-industrial value
    2423c bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
    2524c bl95_b1-input-R-a parameter, may be varied for tests (    -"-      )
    … …  
    7473      LOGICAL ok_aie            ! Apply AIE or not?
    7574     
    76       REAL sulfate(klon, klev)  ! sulfate aerosol mass concentration [ug m-3]
     75      REAL mass_ins_aero(klon, klev)    ! total mass concentration for all indissoluble aerosols[ug m-3]
     76      REAL mass_ins_aero_pi(klon, klev) ! - " - pre-industrial value
    7777      REAL cdnc(klon, klev)     ! cloud droplet number concentration [m-3]
    7878      REAL re(klon, klev)       ! cloud droplet effective radius [um]
    79       REAL sulfate_pi(klon, klev)  ! sulfate aerosol mass concentration [ug m-3] (pre-industrial value)
    8079      REAL cdnc_pi(klon, klev)     ! cloud droplet number concentration [m-3] (pi value)
    8180      REAL re_pi(klon, klev)       ! cloud droplet effective radius [um] (pi value)
    … …  
    108107            !             
    109108            cdnc(i,k) = 10.**(bl95_b0+bl95_b1*
    110      .           log(MAX(sulfate(i,k),1.e-4))/log(10.))*1.e6 !-m-3
     109     .           log(MAX(mass_ins_aero(i,k),1.e-4))/log(10.))*1.e6 !-m-3
    111110            ! Cloud droplet number concentration (CDNC) is restricted
    112111            ! to be within [20, 1000 cm^3]
    … …  
    114113            cdnc(i,k)=MIN(1000.e6,MAX(20.e6,cdnc(i,k)))
    115114            cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1*
    116      .           log(MAX(sulfate_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3
     115     .           log(MAX(mass_ins_aero_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3
    117116            cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k)))
    118117            !           

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_local_var_mod.F90

    r1178 r1179  
    6969      REAL, SAVE, ALLOCATABLE :: d_ts(:,:), d_tr(:,:,:)
    7070      !$OMP THREADPRIVATE(d_ts, d_tr)
     71
     72! diagnostique pour le rayonnement
     73      REAL, SAVE, ALLOCATABLE :: topswad_aero(:),  solswad_aero(:)  ! diag
     74      !$OMP THREADPRIVATE(topswad_aero,solswad_aero)
     75      REAL, SAVE, ALLOCATABLE :: topswai_aero(:),  solswai_aero(:)  ! diag
     76      !$OMP THREADPRIVATE(topswai_aero,solswai_aero)
     77      REAL, SAVE, ALLOCATABLE :: topswad0_aero(:), solswad0_aero(:) ! pas utilise, eventuellment pour diag
     78      !$OMP THREADPRIVATE(topswad0_aero,solswad0_aero)
     79      REAL, SAVE, ALLOCATABLE :: topsw_aero(:,:),  solsw_aero(:,:)  ! pas utilise, eventuellment pour diag
     80      !$OMP THREADPRIVATE(topsw_aero,solsw_aero)
     81      REAL, SAVE, ALLOCATABLE :: topsw0_aero(:,:), solsw0_aero(:,:) ! pas utilise, eventuellment pour diag
     82      !$OMP THREADPRIVATE(topsw0_aero,solsw0_aero)
    7183CONTAINS
    7284
    … …  
    100112      allocate(d_u_lif(klon,klev),d_v_lif(klon,klev))
    101113      allocate(d_ts(klon,klev), d_tr(klon,klev,nbtr))
     114      allocate(topswad_aero(klon), solswad_aero(klon))
     115      allocate(topswai_aero(klon), solswai_aero(klon))
     116      allocate(topswad0_aero(klon), solswad0_aero(klon))
     117      allocate(topsw_aero(klon,9), solsw_aero(klon,9))
     118      allocate(topsw0_aero(klon,9), solsw0_aero(klon,9))
    102119      allocate(d_u_hin(klon,klev),d_v_hin(klon,klev),d_t_hin(klon,klev))
    103120
    … …  
    132149      deallocate(d_u_lif,d_v_lif)
    133150      deallocate(d_ts, d_tr)
     151      deallocate(topswad_aero,solswad_aero)
     152      deallocate(topswai_aero,solswai_aero)
     153      deallocate(topswad0_aero,solswad0_aero)
     154      deallocate(topsw_aero,solsw_aero)
     155      deallocate(topsw0_aero,solsw0_aero)
    134156      deallocate(d_u_hin,d_v_hin,d_t_hin)
    135157

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_state_var_mod.F90

    r1150 r1179  
    274274!$OMP THREADPRIVATE(topswai,solswai)
    275275
     276      REAL,SAVE,ALLOCATABLE :: tau_aero(:,:,:,:), piz_aero(:,:,:,:), cg_aero(:,:,:,:)
     277!$OMP THREADPRIVATE(tau_aero, piz_aero, cg_aero)
    276278      REAL,SAVE,ALLOCATABLE :: ccm(:,:,:)
    277279!$OMP THREADPRIVATE(ccm)
    … …  
    379381      ALLOCATE(topswad(klon), solswad(klon))
    380382      ALLOCATE(topswai(klon), solswai(klon))
    381 
     383      ALLOCATE(tau_aero(klon,klev,9,2),piz_aero(klon,klev,9,2),cg_aero(klon,klev,9,2))
    382384      ALLOCATE(ccm(klon,klev,2))
    383385
    … …  
    466468      deallocate(topswad, solswad)
    467469      deallocate(topswai, solswai)
    468 
     470      deallocate(tau_aero,piz_aero,cg_aero)
    469471      deallocate(ccm)
    470472       

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F

    r1176 r1179  
    1 !
    21! $Id$
    32!
    … …  
    537536c
    538537c Variables propres a la physique
    539 c
    540 c      INTEGER radpas
    541 c      SAVE radpas                 ! frequence d'appel rayonnement
    542 ccccccccc$OMP THREADPRIVATE(radpas)
    543 c
    544 cc      INTEGER iflag_con
    545 c
    546538      INTEGER itap
    547539      SAVE itap                   ! compteur pour la physique
    … …  
    10741066      CHARACTER*4, DIMENSION(9)      :: rfname
    10751067      REAL, DIMENSION(klon)          :: aerindex     ! POLDER aerosol index
    1076       REAL, DIMENSION(klon,klev)     :: maerosol     ! aerosol concentration [ug/m3]
    1077       REAL, DIMENSION(klon,klev)     :: maerosol_pi  ! aerosol concentration [ug/m3] (pre-industrial value)
    1078       REAL, DIMENSION(klon,klev,9,2) :: tau_aero, piz_aero, cg_aero
    1079       REAL, DIMENSION(klon)          :: topswad_aero, solswad_aero   ! diag
    1080       REAL, DIMENSION(klon)          :: topswai_aero, solswai_aero   ! diag
    1081       REAL, DIMENSION(klon)          :: topswad0_aero, solswad0_aero ! pas utilise, eventuellment pour diag
    1082       REAL, DIMENSION(klon,9)        :: topsw_aero, solsw_aero       ! pas utilise
    1083       REAL, DIMENSION(klon,9)        :: topsw0_aero, solsw0_aero     ! pas utilise
    1084 
     1068      REAL, DIMENSION(klon,klev)     :: mass_ins_aero! total mass concentration for all indissoluble aerosols[ug/m3]
     1069      REAL, DIMENSION(klon,klev)     :: mass_ins_aero_pi  ! - " - (pre-industrial value)
    10851070
    10861071      ! Parameters
    … …  
    12281213         tau_overturning_th(:)=0.
    12291214
    1230          IF (config_inca /= 'none') ccm(:,:,:) = 0.
     1215         IF (config_inca /= 'none') THEN
     1216            ! jg : initialisation jusqu'au ces variables sont dans restart
     1217            ccm(:,:,:) = 0.
     1218            tau_aero(:,:,:,:) = 0.
     1219            piz_aero(:,:,:,:) = 0.
     1220            cg_aero(:,:,:,:) = 0.
     1221         END IF
    12311222
    12321223         rnebcon0(:,:) = 0.0
    … …  
    26442635      IF (ok_ade.OR.ok_aie) THEN
    26452636         IF (.NOT. aerosol_couple)
    2646      &        CALL aerosol_optic(
     2637     &        CALL readaerosol_optic(
    26472638     &        debut, new_aod, flag_aerosol, rjourvrai, pdtphys,
    26482639     &        pplay, paprs, t_seri, rhcl,
    2649      &        maerosol, maerosol_pi,
     2640     &        mass_ins_aero, mass_ins_aero_pi,
    26502641     &        tau_aero, piz_aero, cg_aero )
    26512642      ELSE
    … …  
    28292820
    28302821      IF (aerosol_couple) THEN
    2831          maerosol(:,:)    = ccm(:,:,1)
    2832          maerosol_pi(:,:) = ccm(:,:,2)
     2822         mass_ins_aero(:,:)    = ccm(:,:,1)
     2823         mass_ins_aero_pi(:,:) = ccm(:,:,2)
    28332824      END IF
    28342825
    … …  
    28392830     .            flwp, fiwp, flwc, fiwc,
    28402831     e            ok_aie,
    2841      e            maerosol, maerosol_pi,
     2832     e            mass_ins_aero, mass_ins_aero_pi,
    28422833     e            bl95_b0, bl95_b1,
    28432834     s            cldtaupi, re, fl)
    … …  
    28472838     .            cldh, cldl, cldm, cldt, cldq,
    28482839     e            ok_aie,
    2849      e            maerosol, maerosol_pi,
     2840     e            mass_ins_aero, mass_ins_aero_pi,
    28502841     e            bl95_b0, bl95_b1,
    28512842     s            cldtaupi, re, fl)
    … …  
    29222913         
    29232914
    2924       ENDIF
     2915      ENDIF ! aerosol_couple
    29252916      itaprad = 0
    2926       ENDIF
     2917      ENDIF ! MOD(itaprad,radpas)
    29272918      itaprad = itaprad + 1
    29282919

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol.F90

    r1151 r1179  
     1! $Id$
    12!
    2 ! $Header$
     3MODULE readaerosol_mod
     4
     5  REAL, SAVE :: not_valid=-333.
     6
     7CONTAINS
     8
     9SUBROUTINE readaerosol(name_aero, type, iyr_in, klev_src, pt_ap, pt_b, pt_out, psurf, load)
     10
     11!****************************************************************************************
     12! This routine will read the aersosol from file.
    313!
    4 SUBROUTINE readaerosol (id_aero, r_day, first, massvar_p)
     14! Read a year data with get_aero_fromfile depending on aer_type :
     15! - actuel   : read year 1980
     16! - preind   : read natural data
     17! - scenario : read one or two years and do eventually linare time interpolation
     18!
     19! Return pointer, pt_out, to the year read or result from interpolation
     20!****************************************************************************************
     21  USE dimphy
     22
     23  IMPLICIT NONE
     24
     25 INCLUDE "iniprint.h"
     26
     27  ! Input arguments
     28  CHARACTER(len=7), INTENT(IN) :: name_aero
     29  CHARACTER(len=*), INTENT(IN) :: type  ! correspond to aer_type in clesphys.h
     30  INTEGER, INTENT(IN)          :: iyr_in
     31
     32  ! Output
     33  INTEGER, INTENT(OUT)            :: klev_src
     34  REAL, POINTER, DIMENSION(:)     :: pt_ap        ! Pointer for describing the vertical levels     
     35  REAL, POINTER, DIMENSION(:)     :: pt_b         ! Pointer for describing the vertical levels     
     36  REAL, POINTER, DIMENSION(:,:,:) :: pt_out       ! The massvar distributions, DIMENSION(klon, klev_src, 12)
     37  REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf  ! Surface pression for 12 months
     38  REAL, DIMENSION(klon,12), INTENT(OUT) :: load   ! Aerosol mass load in each column for 12 months
     39
     40  ! Local variables
     41  CHARACTER(len=4)                :: cyear
     42  REAL, POINTER, DIMENSION(:,:,:) :: pt_2
     43  REAL, DIMENSION(klon,12)        :: psurf2, load2
     44  REAL                            :: p0           ! Reference pressure
     45  INTEGER                         :: iyr1, iyr2, klev_src2
     46  INTEGER                         :: it, k, i
     47  LOGICAL, PARAMETER              :: lonlyone=.FALSE.
     48
     49!****************************************************************************************
     50! Read data depending on aer_type
     51!
     52!****************************************************************************************
     53
     54  IF (type == 'actuel') THEN
     55! Read and return data for year 1980
     56!****************************************************************************************
     57     cyear='1980'
     58     ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
     59     ! 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     
     62
     63  ELSE IF (type == 'preind') THEN
     64! Read and return data from file with suffix .nat
     65!****************************************************************************************     
     66     cyear='.nat'
     67     ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
     68     ! 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     
     71  ELSE IF (type == 'scenario') THEN
     72! Read data depending on actual year and interpolate if necessary
     73!****************************************************************************************
     74     IF (iyr_in .LT. 1850) THEN
     75        cyear='.nat'
     76        WRITE(lunout,*) 'get_aero 1 iyr_in=', iyr_in,'   ',cyear
     77        ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
     78        ! pt_out has dimensions (klon, klev_src, 12)
     79        CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
     80       
     81     ELSE IF (iyr_in .GE. 2100) THEN
     82        cyear='2100'
     83        WRITE(lunout,*) 'get_aero 2 iyr_in=', iyr_in,'   ',cyear
     84        ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
     85        ! pt_out has dimensions (klon, klev_src, 12)
     86        CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
     87       
     88     ELSE
     89        ! Read data from 2 decades and interpolate to actual year
     90        ! a) from actual 10-yr-period
     91        IF (iyr_in.LT.1900) THEN
     92           iyr1 = 1850
     93           iyr2 = 1900
     94        ELSE IF (iyr_in.GE.1900.AND.iyr_in.LT.1920) THEN
     95           iyr1 = 1900
     96           iyr2 = 1920
     97        ELSE
     98           iyr1 = INT(iyr_in/10)*10
     99           iyr2 = INT(1+iyr_in/10)*10
     100        ENDIF
     101       
     102        WRITE(cyear,'(I4)') iyr1
     103        WRITE(lunout,*) 'get_aero 3 iyr_in=', iyr_in,'   ',cyear
     104        ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
     105        ! pt_out has dimensions (klon, klev_src, 12)
     106        CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
     107       
     108        ! If to read two decades:
     109        IF (.NOT.lonlyone) THEN
     110           
     111           ! b) from the next following one
     112           WRITE(cyear,'(I4)') iyr2
     113           WRITE(lunout,*) 'get_aero 4 iyr_in=', iyr_in,'   ',cyear
     114           
     115           NULLIFY(pt_2)
     116           ! get_aero_fromfile returns pt_2 allocated and initialized with data for 12 month
     117           ! pt_2 has dimensions (klon, klev_src, 12)
     118           CALL get_aero_fromfile(name_aero, cyear, klev_src2, pt_ap, pt_b, p0, pt_2, psurf2, load2)
     119           ! Test for same number of vertical levels
     120           IF (klev_src /= klev_src2) THEN
     121              WRITE(lunout,*) 'Two aerosols files with different number of vertical levels is not allowded'
     122              CALL abort_gcm('readaersosol','Error in number of vertical levels',1)
     123           END IF
     124           
     125           ! Linare interpolate to the actual year:
     126           DO it=1,12
     127              DO k=1,klev_src
     128                 DO i = 1, klon
     129                    pt_out(i,k,it) = &
     130                         pt_out(i,k,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * &
     131                         (pt_out(i,k,it) - pt_2(i,k,it))
     132                 END DO
     133              END DO
     134
     135              DO i = 1, klon
     136                 psurf(i,it) = &
     137                      psurf(i,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * &
     138                      (psurf(i,it) - psurf2(i,it))
     139
     140                 load(i,it) = &
     141                      load(i,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * &
     142                      (load(i,it) - load2(i,it))
     143              END DO
     144           END DO
     145
     146           ! Deallocate pt_2 no more needed
     147           DEALLOCATE(pt_2)
     148           
     149        END IF ! lonlyone
     150     END IF ! iyr_in .LT. 1850
     151
     152  ELSE
     153     WRITE(lunout,*)'This option is not implemented : aer_type = ', type
     154     CALL abort_gcm('readaerosol','Error : aer_type parameter not accepted',1)
     155  END IF ! type
     156
     157
     158END SUBROUTINE readaerosol
     159
     160
     161  SUBROUTINE get_aero_fromfile(varname, cyr, klev_src, pt_ap, pt_b, p0, pt_year, psurf_out, load_out)
     162!****************************************************************************************
     163! Read 12 month aerosol from file and distribute to local process on physical grid.
     164! Vertical levels, klev_src, may differ from model levels if new file format.
     165!
     166! For mpi_root and master thread :
     167! 1) Open file
     168! 2) Find vertical dimension klev_src
     169! 3) Read field month by month
     170! 4) Close file 
     171! 5) Transform the global field from 2D(iim, jjp+1) to 1D(klon_glo)
     172!     - Also the levels and the latitudes have to be inversed
     173!
     174! For all processes and threads :
     175! 6) Scatter global field(klon_glo) to local process domain(klon)
     176! 7) Test for negative values
     177!****************************************************************************************
     178
     179    USE netcdf
     180    USE dimphy
     181    USE mod_grid_phy_lmdz
     182    USE mod_phys_lmdz_para
     183
     184    IMPLICIT NONE
     185     
     186    INCLUDE "dimensions.h"     
     187    INCLUDE "iniprint.h"
     188
     189! Input argumets
     190    CHARACTER(len=7), INTENT(IN)          :: varname
     191    CHARACTER(len=4), INTENT(IN)          :: cyr
     192
     193! Output arguments
     194    INTEGER, INTENT(OUT)                  :: klev_src     ! Number of vertical levels in file
     195    REAL, POINTER, DIMENSION(:)           :: pt_ap        ! Pointer for describing the vertical levels     
     196    REAL, POINTER, DIMENSION(:)           :: pt_b         ! Pointer for describing the vertical levels     
     197    REAL                                  :: p0           ! Reference pressure value
     198    REAL, POINTER, DIMENSION(:,:,:)       :: pt_year      ! Pointer-variabale from file, 12 month, grid : klon,klev_src
     199    REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf_out    ! Surface pression for 12 months
     200    REAL, DIMENSION(klon,12), INTENT(OUT) :: load_out     ! Aerosol mass load in each column
     201
     202! Local variables
     203    CHARACTER(len=30)     :: fname
     204    CHARACTER(len=30)     :: cvar
     205    INTEGER               :: ncid, dimid, varid
     206    INTEGER               :: imth, i, j, k, ierr
     207    REAL                  :: npole, spole
     208    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: varmth
     209    REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: varyear       ! Global variable read from file, 12 month
     210    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: varyear_glo1D !(klon_glo, klev_src, 12)
     211    REAL, ALLOCATABLE, DIMENSION(:)       :: varktmp
     212
     213    REAL, DIMENSION(iim,jjm+1,12)         :: psurf_glo2D   ! Surface pression for 12 months on dynamics global grid
     214    REAL, DIMENSION(klon_glo,12)          :: psurf_glo1D   ! -"- on physical global grid
     215    REAL, DIMENSION(iim,jjm+1,12)         :: load_glo2D    ! Load for 12 months on dynamics global grid
     216    REAL, DIMENSION(klon_glo,12)          :: load_glo1D    ! -"- on physical global grid
     217    REAL, DIMENSION(iim,jjm+1)            :: vartmp
     218    LOGICAL                               :: new_file
     219
     220
     221    ! Deallocate pointers
     222    IF (ASSOCIATED(pt_ap)) DEALLOCATE(pt_ap)
     223    IF (ASSOCIATED(pt_b))  DEALLOCATE(pt_b)
     224
     225!$OMP MASTER
     226    IF (is_mpi_root) THEN
     227
     228! 1) Open file
     229!****************************************************************************************
     230       fname = TRIM(varname)//'.run'//cyr//'.nc'
    5231 
    6   USE dimphy, ONLY : klev
    7   USE mod_grid_phy_lmdz, klon=>klon_glo
    8   USE mod_phys_lmdz_para
     232       WRITE(lunout,*) 'reading ', TRIM(fname)
     233       CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, ncid) )
     234
     235! 2) Check if old or new file is avalabale.
     236!    New type of file should contain the dimension 'lev'
     237!    Old type of file should contain the dimension 'PRESNIVS'
     238!****************************************************************************************
     239       ierr = nf90_inq_dimid(ncid, 'lev', dimid)
     240       IF (ierr /= NF90_NOERR) THEN
     241          ! Coordinate axe lev not found. Check for presnivs.
     242          ierr = nf90_inq_dimid(ncid, 'PRESNIVS', dimid)
     243          IF (ierr /= NF90_NOERR) THEN
     244             ! Dimension PRESNIVS not found either
     245             CALL abort_gcm('get_aero_fromfile', 'dimension lev or presnivs not in file',1)
     246          ELSE
     247             ! Old file found
     248             new_file=.FALSE.
     249             WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will not be done'
     250          END IF
     251       ELSE
     252          ! New file found
     253          new_file=.TRUE.
     254          WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will be done'
     255       END IF
     256       
     257! 2) Find vertical dimension klev_src
     258!****************************************************************************************
     259       CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src) )
     260       
     261     ! Allocate variables depending on the number of vertical levels
     262       ALLOCATE(varmth(iim, jjm+1, klev_src), varyear(iim, jjm+1, klev_src, 12), stat=ierr)
     263       IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 1',1)
     264
     265       ALLOCATE(pt_ap(klev_src), pt_b(klev_src), varktmp(klev_src), stat=ierr)
     266       IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 2',1)
     267
     268! 3) Read all variables from file
     269!    There is 2 options for the file structure :
     270!    new_file=TRUE  : read varyear, ps, pt_ap and pt_b
     271!    new_file=FALSE : read varyear month by month
     272!****************************************************************************************
     273
     274       IF (new_file) THEN
     275
     276! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
     277!****************************************************************************************
     278          ! Get variable id
     279          CALL check_err( nf90_inq_varid(ncid, TRIM(varname), varid) )
     280         
     281          ! Get the variable
     282          CALL check_err( nf90_get_var(ncid, varid, varyear(:,:,:,:)) )
     283         
     284! ++) Read surface pression, 12 month in one variable
     285!****************************************************************************************
     286          ! Get variable id
     287          CALL check_err( nf90_inq_varid(ncid, "ps", varid) )
     288          ! Get the variable
     289          CALL check_err( nf90_get_var(ncid, varid, psurf_glo2D) )
     290         
     291! ++) Read mass load, 12 month in one variable
     292!****************************************************************************************
     293          ! Get variable id
     294          CALL check_err( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) )
     295          ! Get the variable
     296          CALL check_err( nf90_get_var(ncid, varid, load_glo2D) )
     297         
     298! ++) Read ap
     299!****************************************************************************************
     300          ! Get variable id
     301          CALL check_err( nf90_inq_varid(ncid, "ap", varid) )
     302          ! Get the variable
     303          CALL check_err( nf90_get_var(ncid, varid, pt_ap) )
     304
     305! ++) Read b
     306!****************************************************************************************
     307          ! Get variable id
     308          CALL check_err( nf90_inq_varid(ncid, "b", varid) )
     309          ! Get the variable
     310          CALL check_err( nf90_get_var(ncid, varid, pt_b) )
     311
     312! ++) Read p0 : reference pressure
     313!****************************************************************************************
     314          ! Get variable id
     315          CALL check_err( nf90_inq_varid(ncid, "p0", varid) )
     316          ! Get the variable
     317          CALL check_err( nf90_get_var(ncid, varid, p0) )
     318         
     319
     320       ELSE  ! old file
     321
     322! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
     323!****************************************************************************************
     324          DO imth=1, 12
     325             IF (imth.EQ.1) THEN
     326                cvar=TRIM(varname)//'JAN'
     327             ELSE IF (imth.EQ.2) THEN
     328                cvar=TRIM(varname)//'FEB'
     329             ELSE IF (imth.EQ.3) THEN
     330                cvar=TRIM(varname)//'MAR'
     331             ELSE IF (imth.EQ.4) THEN
     332                cvar=TRIM(varname)//'APR'
     333             ELSE IF (imth.EQ.5) THEN
     334                cvar=TRIM(varname)//'MAY'
     335             ELSE IF (imth.EQ.6) THEN
     336                cvar=TRIM(varname)//'JUN'
     337             ELSE IF (imth.EQ.7) THEN
     338                cvar=TRIM(varname)//'JUL'
     339             ELSE IF (imth.EQ.8) THEN
     340                cvar=TRIM(varname)//'AUG'
     341             ELSE IF (imth.EQ.9) THEN
     342                cvar=TRIM(varname)//'SEP'
     343             ELSE IF (imth.EQ.10) THEN
     344                cvar=TRIM(varname)//'OCT'
     345             ELSE IF (imth.EQ.11) THEN
     346                cvar=TRIM(varname)//'NOV'
     347             ELSE IF (imth.EQ.12) THEN
     348                cvar=TRIM(varname)//'DEC'
     349             END IF
     350             
     351             ! Get variable id
     352             CALL check_err( nf90_inq_varid(ncid, TRIM(cvar), varid) )
     353             
     354             ! Get the variable
     355             CALL check_err( nf90_get_var(ncid, varid, varmth) )
     356             
     357             ! Store in variable for the whole year
     358             varyear(:,:,:,imth)=varmth(:,:,:)
     359             
     360          END DO
     361         
     362          ! Putting dummy
     363          psurf_glo2D(:,:,:) = not_valid
     364          load_glo2D(:,:,:)  = not_valid
     365          pt_ap(:) = not_valid
     366          pt_b(:)  = not_valid
     367
     368       END IF
     369
     370! 4) Close file 
     371!****************************************************************************************
     372       CALL check_err( nf90_close(ncid) )
     373     
     374
     375! 5) Transform the global field from 2D(iim, jjp+1) to 1D(klon_glo)
     376!****************************************************************************************
     377! Test if vertical levels have to be inversed
     378
     379       IF ((pt_b(1) < pt_b(klev_src)) .OR. .NOT. new_file) THEN
     380          WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' needs to be inversed'
     381          WRITE(lunout,*) 'before pt_ap = ', pt_ap
     382          WRITE(lunout,*) 'before pt_b = ', pt_b
     383         
     384          ! Inverse vertical levels for varyear
     385          DO imth=1, 12
     386             varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
     387             DO k=1, klev_src
     388                DO j=1, jjm+1
     389                   DO i=1,iim
     390                      varyear(i,j,k,imth) = varmth(i,j,klev_src+1-k)
     391                   END DO
     392                END DO
     393             END DO
     394          END DO
     395           
     396          ! Inverte vertical axes for pt_ap and pt_b
     397          varktmp(:) = pt_ap(:)
     398          DO k=1, klev_src
     399             pt_ap(k) = varktmp(klev_src+1-k)
     400          END DO
     401
     402          varktmp(:) = pt_b(:)
     403          DO k=1, klev_src
     404             pt_b(k) = varktmp(klev_src+1-k)
     405          END DO
     406          WRITE(lunout,*) 'after pt_ap = ', pt_ap
     407          WRITE(lunout,*) 'after pt_b = ', pt_b
     408
     409       ELSE
     410          WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' is ok, no vertical inversion is done'
     411          WRITE(lunout,*) 'pt_ap = ', pt_ap
     412          WRITE(lunout,*) 'pt_b = ', pt_b
     413       END IF
     414
     415!     - Also the latitudes have to be inversed
     416       DO imth=1, 12
     417          ! Inverse latitudes
     418          varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
     419          DO k=1,klev_src
     420             DO j=1,jjm+1
     421                DO i=1,iim
     422                   varyear(i,j,k,imth) = varmth(i,jjm+1+1-j,k)
     423                END DO
     424             END DO
     425          END DO
     426
     427          ! Inverse latitudes for surface pressure
     428          vartmp(:,:) = psurf_glo2D(:,:,imth)
     429          DO j=1, jjm+1
     430             DO i=1,iim
     431                psurf_glo2D(i,j,imth)= vartmp(i,jjm+1+1-j)
     432             END DO
     433          END DO
     434         
     435          ! Inverse latitudes for the load
     436          vartmp(:,:) = load_glo2D(:,:,imth)
     437          DO j=1, jjm+1
     438             DO i=1,iim
     439                load_glo2D(i,j,imth)= vartmp(i,jjm+1+1-j)
     440             END DO
     441          END DO
     442 
     443          ! Do zonal mead at poles and distribut at whole first and last latitude
     444          DO k=1, klev_src
     445             npole=0.  ! North pole, j=1
     446             spole=0.  ! South pole, j=jjm+1         
     447             DO i=1,iim
     448                npole = npole + varyear(i,1,k,imth)
     449                spole = spole + varyear(i,jjm+1,k,imth)
     450             END DO
     451             npole = npole/FLOAT(iim)
     452             spole = spole/FLOAT(iim)
     453             varyear(:,1,    k,imth) = npole
     454             varyear(:,jjm+1,k,imth) = spole
     455          END DO
     456       END DO ! imth
     457
     458       ALLOCATE(varyear_glo1D(klon_glo, klev_src, 12), stat=ierr)
     459       IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 3',1)
     460       
     461       ! Transform from 2D to 1D field
     462       CALL grid2Dto1D_glo(varyear,varyear_glo1D)
     463       CALL grid2Dto1D_glo(psurf_glo2D,psurf_glo1D)
     464       CALL grid2Dto1D_glo(load_glo2D,load_glo1D)
     465       
     466    END IF ! is_mpi_root
     467!$OMP END MASTER BARRIER
    9468 
    10   IMPLICIT NONE
    11      
    12 ! Content:
    13 ! --------
    14 ! This routine reads in monthly mean values of massvar aerosols and
    15 ! returns a linearly interpolated dayly-mean field.     
    16 !
    17 !
    18 ! Author:
    19 ! -------
    20 ! Johannes Quaas (quaas@lmd.jussieu.fr)
    21 ! Celine Deandreis & Anne Cozic
    22 ! 19/02/09
    23 !
    24  
    25 !     
    26 ! Include-files:
    27 ! --------------     
    28   INCLUDE "YOMCST.h"
    29   INCLUDE "chem.h"     
    30   INCLUDE "dimensions.h"     
    31   INCLUDE "temps.h"     
    32   INCLUDE "clesphys.h"
    33   INCLUDE "iniprint.h"
    34 !
    35 ! Input:
    36 ! ------
    37   INTEGER, INTENT(in) :: id_aero
    38   REAL, INTENT(in)    :: r_day        ! Day of integration
    39   LOGICAL, INTENT(in) :: first        ! First timestep
    40                                       ! (and therefore initialization necessary)
    41 !     
    42 ! Output:     
    43 ! -------     
    44   REAL, INTENT(out) ::   massvar_p(klon_omp,klev) ! Mass of massvar (monthly mean data,from file) [ug AIBCM/m3]
    45 
    46 !     
    47 ! Local Variables:
    48 ! ----------------     
    49   INTEGER                         ::  i, ig, k, it
    50   INTEGER                         :: j, iday, iyr, iyr1, iyr2
    51   INTEGER                         :: im, day1, day2, im2
    52   INTEGER, PARAMETER              :: ny=jjm+1
    53   CHARACTER(len=4)                :: cyear
    54   CHARACTER(len=7),DIMENSION(8)   :: name_aero
    55   REAL                            :: var_1(iim, jjm+1, klev, 12)
    56   REAL, DIMENSION(klon,klev)      ::  massvar
    57   REAL, DIMENSION(iim, jjm+1, klev, 12)        :: var_2  ! The massvar distributions
    58   REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE  :: var ! VAR in right dimension
    59   REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE     :: var_out
    60 !$OMP THREADPRIVATE(var,var_out)
    61 
    62   LOGICAL            :: lnewday
    63   LOGICAL, PARAMETER :: lonlyone=.FALSE.
    64   LOGICAL,SAVE       :: first2=.TRUE.
    65 !$OMP THREADPRIVATE(first2)
    66 
    67 ! Variable pour definir a partir d'un numero d'aerosol, son nom
    68   name_aero(1) = "SSSSM  "
    69   name_aero(2) = "ASSSM  "
    70   name_aero(3) = "ASBCM  "
    71   name_aero(4) = "ASPOMM "
    72   name_aero(5) = "SO4    "
    73   name_aero(6) = "CIDUSTM"
    74   name_aero(7) = "AIBCM  "
    75   name_aero(8) = "AIPOMM "
    76 
    77 !$OMP MASTER
    78   IF (first2) THEN
    79       ALLOCATE( var(klon, klev, 12,8) )
    80       ALLOCATE( var_out(klon, klev,8))
    81       first2=.FALSE.
    82 
    83       IF (aer_type /= 'actuel  ' .AND. aer_type /= 'preind  ' .AND.   &
    84            aer_type /= 'scenario') THEN
    85          WRITE(lunout,*)' *** Warning ***'
    86          WRITE(lunout,*)'Option aer_type non prevu pour les aerosols = ',       &
    87               aer_type
    88          CALL abort_gcm('readaerosol','Error : aer_type parameter not accepted',1)
    89       ENDIF
    90    ENDIF
    91 
    92   IF (is_mpi_root) THEN
    93      
    94       iday = INT(r_day)
    95       ! Get the year of the run
    96       iyr  = iday/360
    97       ! Get the day of the actual year:
    98       iday = iday-iyr*360
    99       ! 0.02 is about 0.5/24, namly less than half an hour
    100       lnewday = (r_day-FLOAT(iday).LT.0.02)
    101      
    102       !     ---------------------------------------------
    103       !     All has to be done only, if a new day begins!       
    104       !     ---------------------------------------------
    105          
    106       IF (lnewday.OR.first) THEN
    107          
    108           im = iday/30 +1     ! the actual month
    109           ! annee_ref is the initial year (defined in temps.h)
    110           iyr = iyr + annee_ref
    111          
    112           ! Do I have to read new data? (Is this the first day of a year?)
    113           IF (first.OR.iday.EQ.1.) THEN
    114               ! Initialize values
    115               DO it=1,12
    116                 DO k=1,klev
    117                   DO i=1,klon
    118                     var(i,k,it,id_aero)=0.
    119                   ENDDO
    120                 ENDDO
    121               ENDDO
    122              
    123 
    124               IF (aer_type == 'actuel  ') then
    125 
    126                  cyear='1980'
    127                  CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero))
    128 
    129               ELSE IF (aer_type == 'preind  ') THEN
    130 
    131                  cyear='.nat'
    132                  CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero))
    133 
    134               ELSE ! aer_type=scenario
    135 
    136                  IF (iyr .LT. 1850) THEN
    137                     cyear='.nat'
    138                     WRITE(lunout,*) 'get_aero  iyr=', iyr,'   ',cyear
    139                     CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero))
    140 
    141                  ELSE IF (iyr .GE. 2100) THEN
    142                     cyear='2100'
    143                     WRITE(lunout,*) 'get_aero  iyr=', iyr,'   ',cyear
    144                     CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero))
    145 
    146                  ELSE ! Read data from 2 decades
    147                     ! a) from actual 10-yr-period
    148                     IF (iyr.LT.1900) THEN
    149                        iyr1 = 1850
    150                        iyr2 = 1900
    151                     ELSE IF (iyr.GE.1900.AND.iyr.LT.1920) THEN
    152                        iyr1 = 1900
    153                        iyr2 = 1920
    154                     ELSE
    155                        iyr1 = INT(iyr/10)*10
    156                        iyr2 = INT(1+iyr/10)*10
    157                     ENDIF
    158                    
    159                     WRITE(cyear,'(I4)') iyr1
    160                     WRITE(lunout,*) 'get_aero  iyr=', iyr,'   ',cyear
    161                     CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero))
    162                    
    163                     ! If to read two decades:
    164                     IF (.NOT.lonlyone) THEN
    165                        
    166                        ! b) from the next following one
    167                        WRITE(cyear,'(I4)') iyr2
    168                        WRITE(lunout,*) 'get_aero  iyr=', iyr,'   ',cyear
    169                        
    170                        CALL get_aero_fromfile(cyear, var_2, name_aero(id_aero))
    171                        
    172                        ! Interpolate linarily to the actual year:
    173                        DO it=1,12
    174                           DO k=1,klev
    175                              DO j=1,jjm
    176                                 DO i=1,iim
    177                                    var_1(i,j,k,it) = &
    178                                         var_1(i,j,k,it) - FLOAT(iyr-iyr1)/FLOAT(iyr2-iyr1) * &
    179                                         (var_1(i,j,k,it) - var_2(i,j,k,it))
    180                                 ENDDO
    181                              ENDDO
    182                           ENDDO
    183                        ENDDO
    184                    
    185                     ENDIF ! lonlyone
    186                  ENDIF ! iyr .LT. 1850       
    187               ENDIF ! aer_type
    188                
    189               ! Transform the horizontal 2D-field into the physics-field
    190               ! (Also the levels and the latitudes have to be inversed)
    191              
    192               DO it=1,12
    193                 DO k=1, klev         
    194                   ! a) at the poles, use the zonal mean:
    195                   DO i=1,iim
    196                     ! North pole
    197                     var(1,k,it,id_aero) = &
    198                        var(1,k,it,id_aero)+ var_1(i,jjm+1,klev+1-k,it)
    199                     ! South pole
    200                     var(klon,k,it,id_aero)= &
    201                        var(klon,k,it,id_aero)+ var_1(i,1,klev+1-k,it)
    202                   ENDDO
    203                   var(1,k,it,id_aero)   = var(1,k,it,id_aero)/FLOAT(iim)
    204                   var(klon,k,it,id_aero)= var(klon,k,it,id_aero)/FLOAT(iim)
    205                    
    206                   ! b) the values between the poles:
    207                   ig=1
    208                   DO j=2,jjm
    209                     DO i=1,iim
    210                       ig=ig+1
    211 
    212                       IF (ig.GT.klon)  THEN
    213                           WRITE(lunout,*) 'Attention dans readaerosol', &
    214                              name_aero(id_aero), 'ig > klon', ig, klon
    215                       ENDIF
    216 
    217                       var(ig,k,it,id_aero) =  var_1(i,jjm+1+1-j,klev+1-k,it)
    218 
    219                     ENDDO
    220                   ENDDO
    221                   IF (ig.NE.klon-1) THEN
    222                       PRINT *,'Error in readaerosol', name_aero(id_aero), 'conversion'
    223                       CALL abort_gcm('readaerosol','Error in readaerosol 1',1)
    224                   ENDIF
    225                 ENDDO         ! Loop over k (vertical)
    226               ENDDO           ! Loop over it (months)
    227              
    228           ENDIF               ! Had to read new data?
    229            
    230      
    231           ! Interpolate to actual day:
    232           IF (iday.LT.im*30-15) THEN         
    233               ! in the first half of the month use month before and actual month
    234               im2=im-1
    235               day2 = im2*30-15
    236               day1 = im2*30+15
    237               IF (im2.LE.0) THEN
    238                   ! the month is january, thus the month before december
    239                   im2=12
    240               ENDIF
    241               DO k=1,klev
    242                 DO i=1,klon
    243                   massvar(i,k) = &
    244                      var(i,k,im2,id_aero)- FLOAT(iday-day2)/FLOAT(day1-day2) * &
    245                      (var(i,k,im2,id_aero) - var(i,k,im,id_aero))
    246 
    247                   IF (massvar(i,k).LT.0.) THEN
    248                       IF (iday-day2.LT.0.) WRITE(lunout,*) 'iday-day2',iday-day2
    249                       IF (var(i,k,im2,id_aero) - var(i,k,im,id_aero).LT.0.) THEN
    250                           PRINT*, name_aero(id_aero),'(i,k,im2)- ', &
    251                              name_aero(id_aero),'(i,k,im)',         &
    252                              var(i,k,im2,id_aero) - var(i,k,im,id_aero)
    253                       ENDIF
    254 
    255                       IF (day1-day2.LT.0.) WRITE(lunout,*) 'day1-day2',day1-day2
    256                       WRITE(lunout,*) 'stop',name_aero(id_aero)
    257                       CALL abort_gcm('readaerosol','Error in interpolation 2',1)
    258                   ENDIF
    259                 ENDDO
    260               ENDDO
    261           ELSE
    262               ! the second half of the month
    263               im2=im+1
    264               IF (im2.GT.12) THEN
    265                   ! the month is december, the following thus january
    266                   im2=1
    267               ENDIF
    268               day2 = im*30-15
    269               day1 = im*30+15
    270               DO k=1,klev
    271                 DO i=1,klon
    272                   massvar(i,k) = &
    273                      var(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * &
    274                      (var(i,k,im2,id_aero) - var(i,k,im,id_aero))
    275 
    276                   IF (massvar(i,k).LT.0.) THEN
    277                       IF (iday-day2.LT.0.) &
    278                          WRITE(lunout,*) 'iday-day2',iday-day2
    279                       IF (var(i,k,im2,id_aero) -  var(i,k,im,id_aero).LT.0.) THEN
    280                           WRITE(lunout,*) name_aero(id_aero),'(i,k,im2)-', &
    281                              name_aero(id_aero),'(i,k,im)',           &
    282                              var(i,k,im2,id_aero) - var(i,k,im,id_aero)
    283                       ENDIF
    284                       IF (day1-day2.LT.0.) &
    285                          WRITE(lunout,*) 'day1-day2',day1-day2
    286                       WRITE(lunout,*) 'stop',name_aero(id_aero)
    287                       CALL abort_gcm('readaerosol','Error in interpolation 3',1)
    288                   ENDIF
    289                 ENDDO
    290               ENDDO
    291           ENDIF
    292          
    293      
    294           ! The massvar concentration [molec cm-3] is read in.
    295           ! Convert it into mass [ug VAR/m3]
    296           ! masse_var in [g/mol], n_avogadro in [molec/mol]
    297           ! The sulfate mass [ug VAR/m3] is read in.
    298           DO k=1,klev
    299             DO i=1,klon
    300               var_out(i,k,id_aero) = massvar(i,k)
    301               IF (var_out(i,k,id_aero).LT.0) THEN
    302                   PRINT*, 'Attention il y a un probleme dans readaerosol', &
    303                      name_aero(id_aero), 'la masse est negative'
    304                   CALL abort_gcm('readaerosol','Error in readaerosol 4',1)
    305               ENDIF
    306             ENDDO
    307           ENDDO
    308       ELSE                    ! if no new day, use old data:
    309           DO k=1,klev
    310             DO i=1,klon
    311               massvar(i,k) = var_out(i,k,id_aero)
    312               IF (var_out(i,k,id_aero).LT.0)   THEN
    313                   PRINT*, 'Attention il y a un probleme dans readaerosol', &
    314                      name_aero(id_aero), 'la masse est negative'
    315                   CALL abort_gcm('readaerosol','Error in readaerosol 5',1)
    316               ENDIF
    317             ENDDO
    318           ENDDO
    319          
    320            
    321       ENDIF                   ! Did I have to do anything (was it a new day?)
    322      
    323   ENDIF                      ! phy_rank==0
    324  
    325 !$OMP END MASTER
    326   CALL Scatter(massvar,massvar_p)             
    327 
    328 END SUBROUTINE readaerosol
    329      
    330      
    331      
    332      
    333      
    334 !-----------------------------------------------------------------------------
    335 ! Read in /calculate pre-industrial values of sulfate     
    336 !-----------------------------------------------------------------------------
    337      
    338 SUBROUTINE readaerosol_preind (id_aero, r_day, first, pi_massvar_p)
    339  
    340   USE dimphy, ONLY : klev
    341   USE mod_grid_phy_lmdz, klon=>klon_glo
    342   USE mod_phys_lmdz_para
    343   IMPLICIT NONE
    344      
    345   ! Content:
    346   ! --------
    347   ! This routine reads in monthly mean values of massvar aerosols and
    348   ! returns a linearly interpolated dayly-mean field.     
    349   !
    350   ! It does so for the preindustriel values of the massvar, to a large part
    351   ! analogous to the routine readmassvar above.     
    352   !
    353   ! Only Pb: Variables must be saved and don t have to be overwritten!
    354   !     
    355   ! Author:
    356   ! -------
    357   ! Celine Deandreis & Anne Cozic (LSCE) 2009
    358   ! Johannes Quaas (quaas@lmd.jussieu.fr)
    359   ! 26/06/01
    360   !
    361   ! Modifications:
    362   ! --------------
    363   ! see above
    364 
    365      
    366   ! Include-files:
    367   ! --------------     
    368 
    369   INCLUDE "YOMCST.h"
    370   INCLUDE "chem.h"     
    371   INCLUDE "dimensions.h"     
    372   INCLUDE "temps.h"     
    373   INCLUDE "iniprint.h"
    374  
    375   ! Input:
    376   ! ------
    377   INTEGER, INTENT(in) ::  id_aero
    378   REAL, INTENT(in)    ::  r_day          ! Day of integration
    379   LOGICAL, INTENT(in) ::  first          ! First timestep (and therefore initialization necessary)
    380 
    381 
    382   !     
    383   ! Output:     
    384   ! -------     
    385   REAL, DIMENSION(klon_omp,klev), INTENT(out) ::  pi_massvar_p ! Number conc. massvar (monthly mean data, from file)
    386  
    387 
    388   !     
    389   ! Local Variables:
    390   ! ----------------     
    391   INTEGER                                      :: i, ig, k, it
    392   INTEGER                                      :: j, iday, iyr
    393   INTEGER, PARAMETER                           :: ny=jjm+1
    394   INTEGER                                      :: im, day1, day2, im2
    395 
    396   REAL, DIMENSION(iim, jjm+1, klev, 12)        :: pi_var_1
    397   REAL, DIMENSION(klon,klev)                   :: pi_massvar
    398   REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE  :: pi_var ! VAR in right dimension
    399   REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE    :: pi_var_out
    400 !$OMP THREADPRIVATE(pi_var,pi_var_out)           
    401 
    402   CHARACTER(len=4)                :: cyear
    403   CHARACTER(len=7),DIMENSION(8)   :: name_aero
    404   LOGICAL                         :: lnewday
    405   LOGICAL,SAVE                    :: first2=.TRUE.
    406 !$OMP THREADPRIVATE(first2)
    407  
    408 
    409 ! Variable pour definir a partir d'un numero d'aerosol, son nom
    410   name_aero(1) = "SSSSM  "
    411   name_aero(2) = "ASSSM  "
    412   name_aero(3) = "ASBCM  "
    413   name_aero(4) = "ASPOMM "
    414   name_aero(5) = "SO4    "
    415   name_aero(6) = "CIDUSTM"
    416   name_aero(7) = "AIBCM  "
    417   name_aero(8) = "AIPOMM "
    418 
    419 
    420 !$OMP MASTER
    421   IF (first2) THEN
    422       ALLOCATE( pi_var(klon, klev, 12,8) )
    423       ALLOCATE( pi_var_out(klon, klev,8))
    424       first2=.FALSE.
    425   ENDIF
    426 
    427   IF (is_mpi_root) THEN
    428       iday = INT(r_day)
    429       ! Get the year of the run
    430       iyr  = iday/360
    431       ! Get the day of the actual year:
    432       iday = iday-iyr*360
    433       ! 0.02 is about 0.5/24, namly less than half an hour
    434       lnewday = (r_day-FLOAT(iday).LT.0.02)
    435      
    436       !     ---------------------------------------------
    437       !     All has to be done only, if a new day begins!       
    438       !     ---------------------------------------------
    439      
    440       IF (lnewday.OR.first) THEN
    441          
    442           im = iday/30 +1     ! the actual month
    443           ! annee_ref is the initial year (defined in temps.h)
    444           iyr = iyr + annee_ref     
    445 
    446           IF (first) THEN
    447               cyear='.nat'
    448               CALL get_aero_fromfile(cyear,pi_var_1, name_aero(id_aero))
    449              
    450               ! Transform the horizontal 2D-field into the physics-field
    451               ! (Also the levels and the latitudes have to be inversed)
    452              
    453               ! Initialize field
    454               DO it=1,12
    455                 DO k=1,klev
    456                   DO i=1,klon
    457                     pi_var(i,k,it,id_aero)=0.
    458                   ENDDO
    459                 ENDDO
    460               ENDDO
    461                
    462               WRITE(lunout,*) 'preind: finished reading', FLOAT(iim)
    463               DO it=1,12
    464                 DO k=1, klev         
    465                   ! a) at the poles, use the zonal mean:
    466                   DO i=1,iim
    467                     ! North pole
    468                     pi_var(1,k,it,id_aero)    = &
    469                        pi_var(1,k,it,id_aero) + pi_var_1(i,jjm+1,klev+1-k,it)
    470                     ! South pole
    471                     pi_var(klon,k,it,id_aero) = &
    472                        pi_var(klon,k,it,id_aero) + pi_var_1(i,1,klev+1-k,it)
    473                   ENDDO
    474                   pi_var(1,k,it,id_aero)    = pi_var(1,k,it,id_aero)/FLOAT(iim)
    475                   pi_var(klon,k,it,id_aero) = pi_var(klon,k,it,id_aero)/FLOAT(iim)
    476                  
    477                   ! b) the values between the poles:
    478                   ig=1
    479                   DO j=2,jjm
    480                     DO i=1,iim
    481                       ig=ig+1
    482                       IF (ig.GT.klon)  THEN
    483                           WRITE(lunout,*) 'Attention dans readaerosol_preind', &
    484                              name_aero(id_aero), 'ig > klon', ig, klon
    485                       ENDIF
    486                       pi_var(ig,k,it,id_aero) = &
    487                          pi_var_1(i,jjm+1+1-j,klev+1-k,it)
    488                     ENDDO
    489                   ENDDO
    490                   IF (ig.NE.klon-1) THEN
    491                       WRITE(lunout,*) 'Error in readaerosol_preind (', name_aero(id_aero), ' conversion)'
    492                       CALL abort_gcm('readaerosol_preind','Error in readaerosol_preind 1',1)
    493                   ENDIF
    494                 ENDDO         ! Loop over k (vertical)
    495               ENDDO             ! Loop over it (months)
    496                
    497           ENDIF                ! Had to read new data?
    498            
    499            
    500                                 ! Interpolate to actual day:
    501           IF (iday.LT.im*30-15) THEN         
    502               ! in the first half of the month use month before and actual month
    503               im2=im-1
    504               day1 = im2*30+15
    505               day2 = im2*30-15
    506               IF (im2.LE.0) THEN 
    507                   ! the month is january, thus the month before december
    508                   im2=12
    509               ENDIF
    510               DO k=1,klev
    511                 DO i=1,klon
    512                   pi_massvar(i,k) = &
    513                      pi_var(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * &
    514                      (pi_var(i,k,im2,id_aero) -  pi_var(i,k,im,id_aero))
    515 
    516                   IF (pi_massvar(i,k).LT.0.) THEN
    517                       IF (iday-day2.LT.0.)  WRITE(lunout,*) 'iday-day2',iday-day2
    518                       IF (pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero).LT.0.) THEN
    519                           WRITE(lunout,*) 'pi_',name_aero(id_aero),'(i,k,im2) - pi_', &
    520                              name_aero(id_aero),'(i,k,im)', &
    521                              pi_var(i,k,im2,id_aero) -pi_var(i,k,im,id_aero)
    522                       ENDIF
    523                       IF (day1-day2.LT.0.)WRITE(lunout,*) 'day1-day2',day1-day2
    524                       PRINT *, 'pi_',name_aero(id_aero)
    525                       CALL abort_gcm('readaerosol_preind','Error in readaerosol_preind 2',1)
    526                   ENDIF
    527                 ENDDO
    528               ENDDO
    529           ELSE
    530 
    531               ! the second half of the month
    532               im2=im+1
    533               day1 = im*30+15
    534               IF (im2.GT.12) THEN
    535                   ! the month is december, the following thus january
    536                   im2=1
    537               ENDIF
    538               day2 = im*30-15
    539              
    540               DO k=1,klev
    541                 DO i=1,klon
    542                   pi_massvar(i,k) = &
    543                      pi_var(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * &
    544                      (pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero))
    545                   IF (pi_massvar(i,k).LT.0.) THEN
    546                       IF (iday-day2.LT.0.)  WRITE(lunout,*) 'iday-day2',iday-day2
    547                       IF (pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero).LT.0.) THEN
    548                           WRITE(lunout,*) 'pi_', name_aero(id_aero),'(i,k,im2) - pi_',&
    549                              name_aero(id_aero), '(i,k,im)',&
    550                              pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero)
    551                       ENDIF
    552                       IF (day1-day2.LT.0.) &
    553                          WRITE(lunout,*) 'day1-day2',day1-day2
    554                       PRINT *, 'pi_',name_aero(id_aero)
    555                       CALL abort_gcm('readaerosol_preind','Error in readaerosol_preind 3',1)
    556                   ENDIF
    557                 ENDDO
    558               ENDDO
    559           ENDIF
    560          
    561          
    562           ! The massvar concentration [molec cm-3] is read in.
    563           ! Convert it into mass [ug VAR/m3]
    564           ! masse_var in [g/mol], n_avogadro in [molec/mol]
    565           DO k=1,klev
    566             DO i=1,klon
    567               pi_var_out(i,k,id_aero) = pi_massvar(i,k)
    568             ENDDO
    569           ENDDO
    570          
    571       ELSE                   ! If no new day, use old data:
    572           DO k=1,klev
    573             DO i=1,klon
    574               pi_massvar(i,k) = pi_var_out(i,k,id_aero)           
    575             ENDDO
    576           ENDDO
    577       ENDIF                  ! Was this the beginning of a new day?
    578    ENDIF                     ! is_mpi_root
    579      
    580 !$OMP END MASTER
    581    CALL Scatter(pi_massvar,pi_massvar_p)           
    582      
    583  END SUBROUTINE readaerosol_preind
    584      
    585      
    586      
    587 !-----------------------------------------------------------------------------
    588 ! Routine for reading VAR data from files
    589 !-----------------------------------------------------------------------------
    590            
    591 
    592 SUBROUTINE get_aero_fromfile (cyr, var, varname)
    593   USE dimphy
    594   IMPLICIT NONE
    595      
    596   INCLUDE "netcdf.inc"
    597   INCLUDE "dimensions.h"     
    598   INCLUDE "iniprint.h"
    599 
    600   CHARACTER(len=4), INTENT(in)                       ::  cyr
    601   REAL, DIMENSION(iim, jjm+1, klev, 12), INTENT(out) ::  var
    602   CHARACTER(len=*), INTENT(in)                       ::  varname
    603  
    604 
    605   CHARACTER(len=30)     :: fname
    606   CHARACTER(len=30)     :: cvar
    607   INTEGER, DIMENSION(3) :: START, COUNT
    608   INTEGER               :: STATUS, NCID, VARID
    609   INTEGER               :: imth, i, j, k
    610   INTEGER, PARAMETER    :: ny=jjm+1
    611   REAL, DIMENSION(iim, ny, klev) :: varmth
    612 !
    613 !
    614 ! 
    615   fname = TRIM(varname)//'.run'//cyr//'.cdf'
    616  
    617   WRITE(lunout,*) 'reading ', TRIM(fname)
    618   STATUS = NF_OPEN (TRIM(fname), NF_NOWRITE, NCID)
    619   IF (STATUS .NE. NF_NOERR) THEN
    620       WRITE(lunout,*) 'err in open ',status
    621       CALL abort_gcm('get_aero_fromfile','Error in opening file',1)
    622   ENDIF
    623   DO imth=1, 12
    624     IF (imth.EQ.1) THEN
    625         cvar=TRIM(varname)//'JAN'
    626     ELSEIF (imth.EQ.2) THEN
    627         cvar=TRIM(varname)//'FEB'
    628     ELSEIF (imth.EQ.3) THEN
    629         cvar=TRIM(varname)//'MAR'
    630     ELSEIF (imth.EQ.4) THEN
    631         cvar=TRIM(varname)//'APR'
    632     ELSEIF (imth.EQ.5) THEN
    633         cvar=TRIM(varname)//'MAY'
    634     ELSEIF (imth.EQ.6) THEN
    635         cvar=TRIM(varname)//'JUN'
    636     ELSEIF (imth.EQ.7) THEN
    637         cvar=TRIM(varname)//'JUL'
    638     ELSEIF (imth.EQ.8) THEN
    639         cvar=TRIM(varname)//'AUG'
    640     ELSEIF (imth.EQ.9) THEN
    641         cvar=TRIM(varname)//'SEP'
    642     ELSEIF (imth.EQ.10) THEN
    643         cvar=TRIM(varname)//'OCT'
    644     ELSEIF (imth.EQ.11) THEN
    645         cvar=TRIM(varname)//'NOV'
    646     ELSEIF (imth.EQ.12) THEN
    647         cvar=TRIM(varname)//'DEC'
    648     ENDIF
    649     start(1)=1
    650     start(2)=1
    651     start(3)=1
    652     COUNT(1)=iim
    653     COUNT(2)=ny
    654     COUNT(3)=klev
    655     STATUS = NF_INQ_VARID (NCID, TRIM(cvar), VARID)
    656     WRITE(lunout,*) ncid,imth,TRIM(cvar), varid
    657    
    658     IF (STATUS .NE. NF_NOERR) WRITE(lunout,*) 'err in read ',status     
    659 
    660 #ifdef NC_DOUBLE
    661     status = NF_GET_VARA_DOUBLE(NCID, VARID, START, COUNT,varmth)
    662 #else
    663     status = NF_GET_VARA_REAL(NCID, VARID, START, COUNT, varmth)
    664 #endif
    665     IF (STATUS .NE. NF_NOERR) WRITE(lunout,*) 'err in read data',status
    666          
    667     DO k=1,klev
    668       DO j=1,jjm+1
    669         DO i=1,iim
    670           IF (varmth(i,j,k).LT.0.) THEN
    671               WRITE(lunout,*) 'this is shit'
    672               WRITE(lunout,*) varname,'(',i,j,k,') =',varmth(i,j,k)
    673           ENDIF
    674           var(i,j,k,imth)=varmth(i,j,k)
    675         ENDDO
    676       ENDDO
    677     ENDDO
    678   ENDDO
    679  
    680   STATUS = NF_CLOSE(NCID)
    681   IF (STATUS .NE. NF_NOERR) WRITE(lunout,*) 'err in closing file',status
    682  
    683 
    684 END SUBROUTINE get_aero_fromfile
    685      
    686 
    687 
    688 
    689 
    690 
    691 
    692 
    693 
    694 
    695 
    696 
    697 
    698 
    699 
    700 
     469! 6) Distribute to all processes
     470!    Scatter global field(klon_glo) to local process domain(klon)
     471!    and distribute klev_src to all processes
     472!****************************************************************************************
     473
     474    ! Distribute klev_src
     475    CALL bcast(klev_src)
     476
     477    ! Allocate and distribute pt_ap and pt_b
     478    IF (.NOT. ASSOCIATED(pt_ap)) THEN  ! if pt_ap is allocated also pt_b is allocated
     479       ALLOCATE(pt_ap(klev_src), pt_b(klev_src), stat=ierr)
     480       IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 4',1)
     481    END IF
     482    CALL bcast(pt_ap)
     483    CALL bcast(pt_b)
     484
     485    ! Allocate space for output pointer variable at local process
     486    IF (ASSOCIATED(pt_year)) DEALLOCATE(pt_year)
     487    ALLOCATE(pt_year(klon, klev_src, 12), stat=ierr)
     488    IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 5',1)
     489
     490    ! Scatter global field to local domain at local process
     491    CALL scatter(varyear_glo1D, pt_year)
     492    CALL scatter(psurf_glo1D, psurf_out)
     493    CALL scatter(load_glo1D,  load_out)
     494
     495! 7) Test for negative values
     496!****************************************************************************************
     497    IF (MINVAL(pt_year) < 0.) THEN
     498       WRITE(lunout,*) 'Warning! Negative values read from file :', fname
     499    END IF
     500
     501  END SUBROUTINE get_aero_fromfile
     502
     503
     504  SUBROUTINE check_err(status)
     505    USE netcdf
     506    IMPLICIT NONE
     507
     508    INCLUDE "iniprint.h"
     509    INTEGER, INTENT (IN) :: status
     510
     511    IF (status /= NF90_NOERR) THEN
     512       WRITE(lunout,*) 'Error in get_aero_fromfile ',status
     513       CALL abort_gcm('get_aero_fromfile',trim(nf90_strerror(status)),1)
     514    END IF
     515
     516  END SUBROUTINE check_err
     517
     518
     519END MODULE readaerosol_mod

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol_optic.F90

    r1160 r1179  
    1 SUBROUTINE aerosol_optic(debut, new_aod, flag_aerosol, rjourvrai, pdtphys, &
     1! $Id$
     2!
     3SUBROUTINE readaerosol_optic(debut, new_aod, flag_aerosol, rjourvrai, pdtphys, &
    24     pplay, paprs, t_seri, rhcl, &
    3      maerosol, maerosol_pi, &
     5     mass_ins_aero, mass_ins_aero_pi, &
    46     tau_aero, piz_aero, cg_aero )
     7
     8! This routine will :
     9! 1) recevie the aerosols(already read and interpolated) corresponding to flag_aerosol
     10! 2) calculate the optical properties for the aerosols
     11!
    512 
    613  USE dimphy
    … …  
    815
    916! Input arguments
     17!****************************************************************************************
    1018  LOGICAL, INTENT(IN)                      :: debut
    1119  LOGICAL, INTENT(IN)                      :: new_aod
    … …  
    1927
    2028! Output arguments
    21   REAL, DIMENSION(klon,klev), INTENT(OUT)     :: maerosol
    22   REAL, DIMENSION(klon,klev), INTENT(OUT)     :: maerosol_pi
    23   REAL, DIMENSION(klon,klev,9,2), INTENT(OUT) :: tau_aero
    24   REAL, DIMENSION(klon,klev,9,2), INTENT(OUT) :: piz_aero
    25   REAL, DIMENSION(klon,klev,9,2), INTENT(OUT) :: cg_aero
     29!****************************************************************************************
     30  REAL, DIMENSION(klon,klev), INTENT(OUT)     :: mass_ins_aero    ! Total mass for all indissoluble aerosols
     31  REAL, DIMENSION(klon,klev), INTENT(OUT)     :: mass_ins_aero_pi !     -"-     preindustrial values
     32  REAL, DIMENSION(klon,klev,9,2), INTENT(OUT) :: tau_aero    ! Aerosol optical thickness
     33  REAL, DIMENSION(klon,klev,9,2), INTENT(OUT) :: piz_aero    ! Single scattering albedo aerosol
     34  REAL, DIMENSION(klon,klev,9,2), INTENT(OUT) :: cg_aero     ! asymmetry parameter aerosol
    2635
    2736! Local variables
     37!****************************************************************************************
    2838  REAL, DIMENSION(klon)        :: aerindex ! POLDER aerosol index
    2939  REAL, DIMENSION(klon,10)     :: fractnat_allaer
    … …  
    4151  REAL, DIMENSION(klon,klev,8) :: m_allaer
    4252
    43   INTEGER :: k
     53  INTEGER :: k, i
    4454 
    45 !
    46 ! Read aersosols from file
    47 !
    48   maerosol(:,:) = 0.
    49   maerosol_pi(:,:) = 0.
    50 
    51   bcsol(:,:) = 0.
    52   bcins(:,:) = 0.
    53   pomsol(:,:) = 0.
    54   pomins(:,:) = 0.
    55   sulfate(:,:) = 0.
    56 
    57 ! Read sulfate
     55!****************************************************************************************
     56! 1) Get aerosol mass
     57!   
     58!****************************************************************************************
     59! Read and interpolate sulfate
    5860  IF ( flag_aerosol .EQ. 1 .OR. &
    5961       flag_aerosol .EQ. 4 .OR. &
    6062       flag_aerosol .EQ. 6 ) THEN
    6163
    62      CALL readaerosol(5,rjourvrai, debut, sulfate)
    63      CALL readaerosol_preind(5,rjourvrai, &
    64           debut, sulfate_pi)
     64     CALL readaerosol_interp(5, rjourvrai, debut, pplay, paprs, sulfate, sulfate_pi)
     65  ELSE
     66     sulfate(:,:) = 0. ; sulfate_pi(:,:) = 0.
     67  END IF
    6568
    66      maerosol(:,:) = maerosol(:,:) + sulfate(:,:)
    67      maerosol_pi(:,:) = maerosol_pi(:,:) + sulfate_pi(:,:)
    68   ENDIF
    69 
    70 
    71 ! Read bcscol and bcins
     69! Read and interpolate bcsol and bcins
    7270  IF ( flag_aerosol .EQ. 2 .OR. &
    7371       flag_aerosol .EQ. 4 .OR. &
    … …  
    7573
    7674     ! Get bc aerosol distribution
    77      CALL readaerosol(3,rjourvrai, debut,bcsol )
    78      CALL readaerosol_preind(3,rjourvrai, debut, bcsol_pi)
    79 
    80      CALL readaerosol(7,rjourvrai, debut,bcins )
    81      CALL readaerosol_preind(7,rjourvrai, debut, bcins_pi)
    82 
    83      maerosol(:,:) = maerosol(:,:) + bcsol(:,:)
    84      maerosol_pi(:,:) = maerosol_pi(:,:) + bcsol_pi(:,:)
    85   ENDIF
     75     CALL readaerosol_interp(3, rjourvrai, debut, pplay, paprs, bcsol, bcsol_pi )
     76     CALL readaerosol_interp(7, rjourvrai, debut, pplay, paprs, bcins, bcins_pi )
     77  ELSE
     78     bcsol(:,:) = 0. ; bcsol_pi(:,:) = 0.
     79     bcins(:,:) = 0. ; bcins_pi(:,:) = 0.
     80  END IF
    8681
    8782
    88 ! Read pomsol and pomins
     83! Read and interpolate pomsol and pomins
    8984  IF ( flag_aerosol .EQ. 3 .OR. &
    9085       flag_aerosol .EQ. 4 .OR. &
    … …  
    9287       flag_aerosol .EQ. 6 ) THEN
    9388
    94      CALL readaerosol(4,rjourvrai, debut,pomsol)
    95      CALL readaerosol_preind(4,rjourvrai,debut,pomsol_pi )
    96 
    97      CALL readaerosol(8,rjourvrai, debut,pomins)
    98      CALL readaerosol_preind(8,rjourvrai,debut,pomins_pi )
    99 
    100      maerosol(:,:) = maerosol(:,:) + pomsol
    101      maerosol_pi(:,:) = maerosol_pi(:,:) + pomsol_pi
    102   ENDIF
     89     CALL readaerosol_interp(4, rjourvrai, debut, pplay, paprs, pomsol, pomsol_pi)
     90     CALL readaerosol_interp(8, rjourvrai, debut, pplay, paprs, pomins, pomins_pi)
     91  ELSE
     92     pomsol(:,:) = 0. ; pomsol_pi(:,:) = 0.
     93     pomins(:,:) = 0. ; pomins_pi(:,:) = 0.
     94  END IF
    10395
    10496
    105 ! Save all aerosols in one variable
     97! Store all aerosols in one variable
    10698!
    10799! ACo pour couplage aerosol offline 07/04/2009
    … …  
    110102! tableaux lorsque les routines de lectures seront
    111103! ajoutees.
    112   m_allaer(:,:,1) = 0.                ! SSSSM || CSSSM
     104  m_allaer(:,:,1) = 0.                ! SSSSM || CSSSM ! Coarse Soluble Sea Salt Mass
    113105  m_allaer(:,:,2) = 0.                ! ASSSM
    114106  m_allaer(:,:,3) = bcsol(:,:)        ! ASBCM
    115107  m_allaer(:,:,4) = pomsol(:,:)       ! ASPOMM
    116108  m_allaer(:,:,5) = sulfate(:,:)      ! ASSO4M || CSSO4M   
    117   m_allaer(:,:,6) = 0.                ! CIDUSTM
     109  m_allaer(:,:,6) = 0.                ! CIDUSTM        ! Coarse Insoluble DUST Mass
    118110  m_allaer(:,:,7) = bcins(:,:)        ! AIBCM
    119111  m_allaer(:,:,8) = pomins(:,:)       ! AIPOMM
    120112
    121113!
    122 ! Calculate the optical properties for the aerosols
     114! Calculate the total mass of all indissoluble aersosols
    123115!
     116  mass_ins_aero(:,:)    = sulfate(:,:)    + bcsol(:,:)    + pomsol(:,:)
     117  mass_ins_aero_pi(:,:) = sulfate_pi(:,:) + bcsol_pi(:,:) + pomsol_pi(:,:)
     118
     119!****************************************************************************************
     120! 2) Calculate optical properties for the aerosols
     121!
     122!****************************************************************************************
    124123  DO k = 1, klev
    125      pdel(:,k) = paprs(:,k) - paprs (:,k+1)
     124     DO i = 1, klon
     125        pdel(i,k) = paprs(i,k) - paprs (i,k+1)
     126     END DO
    126127  END DO
    127128
    128   IF (.NOT. new_aod) THEN
    129      CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl, &
    130           tau_aero(:,:,5,:), piz_aero(:,:,5,:), cg_aero(:,:,5,:), aerindex)
    131   ELSE
     129  IF (new_aod) THEN
     130
    132131     fractnat_allaer(:,:) = 0.
    133132     CALL aeropt_2bands( &
    … …  
    136135          fractnat_allaer, flag_aerosol, &
    137136          pplay, t_seri)
    138 
     137     
     138     ! aeropt_5wv only for validation and diagnostics.
     139     ! In this version no diagnostics are set.
     140     ! jg : may be desactivated if no diagnostics added.
    139141     CALL aeropt_5wv( &
    140142          pdel, m_allaer, &
    141143          pdtphys, rhcl, aerindex, &
    142144          flag_aerosol, pplay, t_seri)
    143   ENDIF
     145  ELSE
    144146
    145 END SUBROUTINE aerosol_optic
     147     CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl, &
     148          tau_aero(:,:,5,:), piz_aero(:,:,5,:), cg_aero(:,:,5,:), aerindex)
     149     
     150  END IF
     151
     152END SUBROUTINE readaerosol_optic
Note: See TracChangeset for help on using the changeset viewer.