Ignore:
Timestamp:
Apr 17, 2009, 5:34:01 PM (15 years ago)
Author:
jghattas
Message:

Ajoute possiblilite de forcer avec d'autre aersols avec le meme principe que pour les so4.

Par default la convergence est rempu avec flag_aersosol=1 + ok_ade=ok_aie=y. Exactement les memes resultats peuvent etre retrouver avec new_aod=n. Les resultats sont exactement les memes sans prise en compte aucun aerosol, avec ok_ade=ok_aie=n.

  • flag_aerosol indique quel aersosol ou combinasion a utiliser (=1 uniquement les SO4 comme avant)
  • les fichiers d'entree so4.runXXXX.cdf change du nom pour majescule SO4.runXXXX.cdf
  • readsulfate change du nom pour readaerosol qui trait tous les aerosols
  • radlwsw change du nom pour radlwsw_aero.
  • aeropt_2bands.F90 et aeropt_5wv.F90 correspond a un reecriture de aeropt.F (premier et deuxieme moitie respectivement)
  • aeropt.F est gardé pour retrouver la convergence si demande avec new_aod=false (=true default)
  • sw_aero est un version evolue de sw_lmdar4 (dans fichier radiation_ar4.F)

ACo + JG

File:
1 moved

Legend:

Unmodified
Added
Removed

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

    r1149 r1150  
    1 !
    2 ! $Header$
    3 !
    4       SUBROUTINE radlwsw(dist, rmu0, fract,
    5      .                  paprs, pplay,tsol,alb1, alb2, t,q,wo,
    6      .                  cldfra, cldemi, cldtaupd,
    7      .                  heat,heat0,cool,cool0,radsol,albpla,
    8      .                  topsw,toplw,solsw,sollw,
    9      .                  sollwdown,
    10      .                  topsw0,toplw0,solsw0,sollw0,
    11      .                  lwdn0, lwdn, lwup0, lwup,
    12      .                  swdn0, swdn, swup0, swup,
    13      .                  ok_ade, ok_aie,
    14      .                  tau_ae, piz_ae, cg_ae,
    15      .                  topswad, solswad,
    16      .                  cldtaupi, topswai, solswai,qsat,flwc,fiwc)
    17 c     
    18       USE dimphy
    19       IMPLICIT none
    20 c======================================================================
    21 c Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
    22 c Objet: interface entre le modele et les rayonnements
    23 c Arguments:
    24 c dist-----input-R- distance astronomique terre-soleil
    25 c rmu0-----input-R- cosinus de l'angle zenithal
    26 c fract----input-R- duree d'ensoleillement normalisee
    27 c co2_ppm--input-R- concentration du gaz carbonique (en ppm)
    28 c solaire--input-R- constante solaire (W/m**2)
    29 c paprs----input-R- pression a inter-couche (Pa)
    30 c pplay----input-R- pression au milieu de couche (Pa)
    31 c tsol-----input-R- temperature du sol (en K)
    32 c alb1-----input-R- albedo du sol(entre 0 et 1) dans l'interval visible
    33 c alb2-----input-R- albedo du sol(entre 0 et 1) dans l'interval proche infra-rouge
    34 c t--------input-R- temperature (K)
    35 c q--------input-R- vapeur d'eau (en kg/kg)
    36 c wo-------input-R- contenu en ozone (en kg/kg) correction MPL 100505
    37 c cldfra---input-R- fraction nuageuse (entre 0 et 1)
    38 c cldtaupd---input-R- epaisseur optique des nuages dans le visible (present-day value)
    39 c cldemi---input-R- emissivite des nuages dans l'IR (entre 0 et 1)
    40 c ok_ade---input-L- apply the Aerosol Direct Effect or not?
    41 c ok_aie---input-L- apply the Aerosol Indirect Effect or not?
    42 c tau_ae, piz_ae, cg_ae-input-R- aerosol optical properties (calculated in aeropt.F)
    43 c cldtaupi-input-R- epaisseur optique des nuages dans le visible
    44 c                   calculated for pre-industrial (pi) aerosol concentrations, i.e. with smaller
    45 c                   droplet concentration, thus larger droplets, thus generally cdltaupi cldtaupd
    46 c                   it is needed for the diagnostics of the aerosol indirect radiative forcing     
    47 c
    48 c heat-----output-R- echauffement atmospherique (visible) (K/jour)
    49 c cool-----output-R- refroidissement dans l'IR (K/jour)
    50 c radsol---output-R- bilan radiatif net au sol (W/m**2) (+ vers le bas)
    51 c albpla---output-R- albedo planetaire (entre 0 et 1)
    52 c topsw----output-R- flux solaire net au sommet de l'atm.
    53 c toplw----output-R- ray. IR montant au sommet de l'atmosphere
    54 c solsw----output-R- flux solaire net a la surface
    55 c sollw----output-R- ray. IR montant a la surface
    56 c solswad---output-R- ray. solaire net absorbe a la surface (aerosol dir)
    57 c topswad---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol dir)
    58 c solswai---output-R- ray. solaire net absorbe a la surface (aerosol ind)
    59 c topswai---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol ind)
    60 c
    61 c ATTENTION: swai and swad have to be interpreted in the following manner:
    62 c ---------
    63 c ok_ade=F & ok_aie=F -both are zero
    64 c ok_ade=T & ok_aie=F -aerosol direct forcing is F_{AD} = topsw-topswad
    65 c                        indirect is zero
    66 c ok_ade=F & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
    67 c                        direct is zero
    68 c ok_ade=T & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
    69 c                        aerosol direct forcing is F_{AD} = topswai-topswad
    70 c
    71      
    72 c======================================================================
    73 cym#include "dimensions.h"
    74 cym#include "dimphy.h"
    75 cym#include "raddim.h"
    76 #include "YOETHF.h"
    77 c
    78       real rmu0(klon), fract(klon), dist
    79 cIM   real co2_ppm
    80 cIM   real solaire
    81 #include "clesphys.h"
    82 c
    83       real paprs(klon,klev+1), pplay(klon,klev)
    84       real alb1(klon), alb2(klon), tsol(klon)
    85       real t(klon,klev), q(klon,klev), wo(klon,klev)
    86       real cldfra(klon,klev), cldemi(klon,klev), cldtaupd(klon,klev)
    87       real heat(klon,klev), cool(klon,klev)
    88       real heat0(klon,klev), cool0(klon,klev)
    89       real radsol(klon), topsw(klon), toplw(klon)
    90       real solsw(klon), sollw(klon), albpla(klon)
    91       real topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
    92       real sollwdown(klon)
    93 cIM output 3D
    94       REAL*8 ZFSUP(KDLON,KFLEV+1)
    95       REAL*8 ZFSDN(KDLON,KFLEV+1)
    96       REAL*8 ZFSUP0(KDLON,KFLEV+1)
    97       REAL*8 ZFSDN0(KDLON,KFLEV+1)
    98 c
    99       REAL*8 ZFLUP(KDLON,KFLEV+1)
    100       REAL*8 ZFLDN(KDLON,KFLEV+1)
    101       REAL*8 ZFLUP0(KDLON,KFLEV+1)
    102       REAL*8 ZFLDN0(KDLON,KFLEV+1)
    103 c
    104       REAL*8 zx_alpha1, zx_alpha2
    105 c
    106 #include "YOMCST.h"
    107 c
    108       INTEGER k, kk, i, j, iof, nb_gr
    109       EXTERNAL LW_LMDAR4,SW_LMDAR4
    110 c
    111 cIM ctes ds clesphys.h  REAL*8 RCO2, RCH4, RN2O, RCFC11, RCFC12
    112       REAL*8 PSCT
    113 c
    114       REAL*8 PALBD(kdlon,2), PALBP(kdlon,2)
    115       REAL*8 PEMIS(kdlon), PDT0(kdlon), PVIEW(kdlon)
    116       REAL*8 PPSOL(kdlon), PDP(kdlon,klev)
    117       REAL*8 PTL(kdlon,kflev+1), PPMB(kdlon,kflev+1)
    118       REAL*8 PTAVE(kdlon,kflev)
    119       REAL*8 PWV(kdlon,kflev), PQS(kdlon,kflev), POZON(kdlon,kflev)
    120       REAL*8 PAER(kdlon,kflev,5)
    121       REAL*8 PCLDLD(kdlon,kflev)
    122       REAL*8 PCLDLU(kdlon,kflev)
    123       REAL*8 PCLDSW(kdlon,kflev)
    124       REAL*8 PTAU(kdlon,2,kflev)
    125       REAL*8 POMEGA(kdlon,2,kflev)
    126       REAL*8 PCG(kdlon,2,kflev)
    127 c
    128       REAL*8 zfract(kdlon), zrmu0(kdlon), zdist
    129 c
    130       REAL*8 zheat(kdlon,kflev), zcool(kdlon,kflev)
    131       REAL*8 zheat0(kdlon,kflev), zcool0(kdlon,kflev)
    132       REAL*8 ztopsw(kdlon), ztoplw(kdlon)
    133       REAL*8 zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon)
    134 cIM
    135       REAL*8 zsollwdown(kdlon)
    136 c
    137       REAL*8 ztopsw0(kdlon), ztoplw0(kdlon)
    138       REAL*8 zsolsw0(kdlon), zsollw0(kdlon)
    139       REAL*8 zznormcp
    140 cIM output 3D : SWup, SWdn, LWup, LWdn
    141       REAL swdn(klon,kflev+1),swdn0(klon,kflev+1)
    142       REAL swup(klon,kflev+1),swup0(klon,kflev+1)
    143       REAL lwdn(klon,kflev+1),lwdn0(klon,kflev+1)
    144       REAL lwup(klon,kflev+1),lwup0(klon,kflev+1)
    145       REAL qsat(klon,klev),flwc(klon,klev),fiwc(klon,klev)
    146 c-OB
    147 cjq the following quantities are needed for the aerosol radiative forcings
    148 
    149       real topswad(klon), solswad(klon) ! output: aerosol direct forcing at TOA and surface
    150       real topswai(klon), solswai(klon) ! output: aerosol indirect forcing atTOA and surface
    151       real tau_ae(klon,klev,2), piz_ae(klon,klev,2), cg_ae(klon,klev,2) ! aerosol optical properties (see aeropt.F)
    152       real cldtaupi(klon,klev)  ! cloud optical thickness for pre-industrial aerosol concentrations
    153                                 ! (i.e., with a smaller droplet concentrationand thus larger droplet radii)
    154       logical ok_ade, ok_aie    ! switches whether to use aerosol direct (indirect) effects or not
    155       real*8 tauae(kdlon,kflev,2) ! aer opt properties
    156       real*8 pizae(kdlon,kflev,2)
    157       real*8 cgae(kdlon,kflev,2)
    158       REAL*8 PTAUA(kdlon,2,kflev) ! present-day value of cloud opt thickness (PTAU is pre-industrial value), local use
    159       REAL*8 POMEGAA(kdlon,2,kflev) ! dito for single scatt albedo
    160       REAL*8 ztopswad(kdlon), zsolswad(kdlon) ! Aerosol direct forcing at TOAand surface
    161       REAL*8 ztopswai(kdlon), zsolswai(kdlon) ! dito, indirect
    162 cjq-end
    163 !rv
    164       tauae(:,:,:)=0.
    165       pizae(:,:,:)=0.
    166       cgae(:,:,:)=0.
    167 !rv
    168      
    169 c
    170 c-------------------------------------------
    171       nb_gr = klon / kdlon
    172       IF (nb_gr*kdlon .NE. klon) THEN
    173          PRINT*, "kdlon mauvais:", klon, kdlon, nb_gr
    174          CALL abort
    175       ENDIF
    176       IF (kflev .NE. klev) THEN
    177           PRINT*, "kflev differe de klev, kflev, klev"
    178           CALL abort
    179       ENDIF
    180 c-------------------------------------------
    181       DO k = 1, klev
    182       DO i = 1, klon
    183          heat(i,k)=0.
    184          cool(i,k)=0.
    185          heat0(i,k)=0.
    186          cool0(i,k)=0.
    187       ENDDO
    188       ENDDO
    189 c
    190       zdist = dist
    191 c
    192 cIM anciennes valeurs
    193 c     RCO2 = co2_ppm * 1.0e-06  * 44.011/28.97
    194 c
    195 cIM : on met RCO2, RCH4, RN2O, RCFC11 et RCFC12 dans clesphys.h /lecture ds conf_phys.F90
    196 c     RCH4 = 1.65E-06* 16.043/28.97
    197 c     RN2O = 306.E-09* 44.013/28.97
    198 c     RCFC11 = 280.E-12* 137.3686/28.97
    199 c     RCFC12 = 484.E-12* 120.9140/28.97
    200 cIM anciennes valeurs
    201 c     RCH4 = 1.72E-06* 16.043/28.97
    202 c     RN2O = 310.E-09* 44.013/28.97
    203 c
    204 c     PRINT*,'IMradlwsw : solaire, co2= ', solaire, co2_ppm
    205       PSCT = solaire/zdist/zdist
    206 c
    207       DO 99999 j = 1, nb_gr
    208       iof = kdlon*(j-1)
    209 c
     1SUBROUTINE radlwsw_aero( &
     2   dist, rmu0, fract, solaire, &
     3   paprs, pplay,tsol,albedo, alblw, &
     4   t,q,wo,&
     5   cldfra, cldemi, cldtaupd,&
     6   ok_ade, ok_aie,&
     7   tau_aero, piz_aero, cg_aero,&
     8   cldtaupi, new_aod, &
     9   heat,heat0,cool,cool0,radsol,albpla,&
     10   topsw,toplw,solsw,sollw,&
     11   sollwdown,&
     12   topsw0,toplw0,solsw0,sollw0,&
     13   lwdn0, lwdn, lwup0, lwup,&
     14   swdn0, swdn, swup0, swup,&
     15   topswad_aero, solswad_aero,&
     16   topswai_aero, solswai_aero, &
     17   topswad0_aero, solswad0_aero,&
     18   topsw_aero, topsw0_aero,&
     19   solsw_aero, solsw0_aero)
     20
     21
     22
     23  USE DIMPHY
     24
     25  IMPLICIT NONE
     26  !======================================================================
     27  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
     28  ! Objet: interface entre le modele et les rayonnements
     29  ! Arguments:
     30  ! dist-----input-R- distance astronomique terre-soleil
     31  ! rmu0-----input-R- cosinus de l'angle zenithal
     32  ! fract----input-R- duree d'ensoleillement normalisee
     33  ! co2_ppm--input-R- concentration du gaz carbonique (en ppm)
     34  ! solaire--input-R- constante solaire (W/m**2)
     35  ! paprs----input-R- pression a inter-couche (Pa)
     36  ! pplay----input-R- pression au milieu de couche (Pa)
     37  ! tsol-----input-R- temperature du sol (en K)
     38  ! albedo---input-R- albedo du sol (entre 0 et 1)
     39  ! t--------input-R- temperature (K)
     40  ! q--------input-R- vapeur d'eau (en kg/kg)
     41  ! wo-------input-R- contenu en ozone (en kg/kg) correction MPL 100505
     42  ! cldfra---input-R- fraction nuageuse (entre 0 et 1)
     43  ! cldtaupd---input-R- epaisseur optique des nuages dans le visible (present-day value)
     44  ! cldemi---input-R- emissivite des nuages dans l'IR (entre 0 et 1)
     45  ! ok_ade---input-L- apply the Aerosol Direct Effect or not?
     46  ! ok_aie---input-L- apply the Aerosol Indirect Effect or not?
     47  ! tau_ae, piz_ae, cg_ae-input-R- aerosol optical properties (calculated in aeropt.F)
     48  ! cldtaupi-input-R- epaisseur optique des nuages dans le visible
     49  !                   calculated for pre-industrial (pi) aerosol concentrations, i.e. with smaller
     50  !                   droplet concentration, thus larger droplets, thus generally cdltaupi cldtaupd
     51  !                   it is needed for the diagnostics of the aerosol indirect radiative forcing     
     52  !
     53  ! heat-----output-R- echauffement atmospherique (visible) (K/jour)
     54  ! cool-----output-R- refroidissement dans l'IR (K/jour)
     55  ! radsol---output-R- bilan radiatif net au sol (W/m**2) (+ vers le bas)
     56  ! albpla---output-R- albedo planetaire (entre 0 et 1)
     57  ! topsw----output-R- flux solaire net au sommet de l'atm.
     58  ! toplw----output-R- ray. IR montant au sommet de l'atmosphere
     59  ! solsw----output-R- flux solaire net a la surface
     60  ! sollw----output-R- ray. IR montant a la surface
     61  ! solswad---output-R- ray. solaire net absorbe a la surface (aerosol dir)
     62  ! topswad---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol dir)
     63  ! solswai---output-R- ray. solaire net absorbe a la surface (aerosol ind)
     64  ! topswai---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol ind)
     65  !
     66  ! ATTENTION: swai and swad have to be interpreted in the following manner:
     67  ! ---------
     68  ! ok_ade=F & ok_aie=F -both are zero
     69  ! ok_ade=T & ok_aie=F -aerosol direct forcing is F_{AD} = topsw-topswad
     70  !                        indirect is zero
     71  ! ok_ade=F & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
     72  !                        direct is zero
     73  ! ok_ade=T & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
     74  !                        aerosol direct forcing is F_{AD} = topswai-topswad
     75  !
     76 
     77  !======================================================================
     78 
     79  ! ====================================================================
     80  ! Adapte au modele de chimie INCA par Celine Deandreis & Anne Cozic -- 2009
     81  ! 1 = ZERO   
     82  ! 2 = AER total   
     83  ! 3 = NAT   
     84  ! 4 = BC   
     85  ! 5 = SO4   
     86  ! 6 = POM   
     87  ! 7 = DUST   
     88  ! 8 = SS   
     89  ! 9 = NO3   
     90  !
     91  ! ====================================================================
     92  include "YOETHF.h"
     93  include "YOMCST.h"
     94
     95! Input arguments
     96  REAL,    INTENT(in)  :: solaire
     97  REAL,    INTENT(in)  :: dist
     98  REAL,    INTENT(in)  :: rmu0(KLON), fract(KLON)
     99  REAL,    INTENT(in)  :: paprs(KLON,KLEV+1), pplay(KLON,KLEV)
     100  REAL,    INTENT(in)  :: albedo(KLON), alblw(KLON), tsol(KLON)
     101  REAL,    INTENT(in)  :: t(KLON,KLEV), q(KLON,KLEV), wo(KLON,KLEV)
     102  LOGICAL, INTENT(in)  :: ok_ade, ok_aie                                 ! switches whether to use aerosol direct (indirect) effects or not
     103  REAL,    INTENT(in)  :: cldfra(KLON,KLEV), cldemi(KLON,KLEV), cldtaupd(KLON,KLEV)
     104  REAL,    INTENT(in)  :: tau_aero(KLON,KLEV,9,2)                        ! aerosol optical properties (see aeropt.F)
     105  REAL,    INTENT(in)  :: piz_aero(KLON,KLEV,9,2)                        ! aerosol optical properties (see aeropt.F)
     106  REAL,    INTENT(in)  :: cg_aero(KLON,KLEV,9,2)                         ! aerosol optical properties (see aeropt.F)
     107  REAL,    INTENT(in)  :: cldtaupi(KLON,KLEV)                            ! cloud optical thickness for pre-industrial aerosol concentrations
     108  LOGICAL, INTENT(in)  :: new_aod                                        ! flag pour retrouver les resultats exacts de l'AR4 dans le cas ou l'on ne travaille qu'avec les sulfates
     109
     110! Output arguments
     111  REAL,    INTENT(out) :: heat(KLON,KLEV), cool(KLON,KLEV)
     112  REAL,    INTENT(out) :: heat0(KLON,KLEV), cool0(KLON,KLEV)
     113  REAL,    INTENT(out) :: radsol(KLON), topsw(KLON), toplw(KLON)
     114  REAL,    INTENT(out) :: solsw(KLON), sollw(KLON), albpla(KLON)
     115  REAL,    INTENT(out) :: topsw0(KLON), toplw0(KLON), solsw0(KLON), sollw0(KLON)
     116  REAL,    INTENT(out) :: sollwdown(KLON)
     117  REAL,    INTENT(out) :: swdn(KLON,kflev+1),swdn0(KLON,kflev+1)
     118  REAL,    INTENT(out) :: swup(KLON,kflev+1),swup0(KLON,kflev+1)
     119  REAL,    INTENT(out) :: lwdn(KLON,kflev+1),lwdn0(KLON,kflev+1)
     120  REAL,    INTENT(out) :: lwup(KLON,kflev+1),lwup0(KLON,kflev+1)
     121  REAL,    INTENT(out) :: topswad_aero(KLON), solswad_aero(KLON)         ! output: aerosol direct forcing at TOA and surface
     122  REAL,    INTENT(out) :: topswai_aero(KLON), solswai_aero(KLON)         ! output: aerosol indirect forcing atTOA and surface
     123  REAL, DIMENSION(klon), INTENT(out)    :: topswad0_aero
     124  REAL, DIMENSION(klon), INTENT(out)    :: solswad0_aero
     125  REAL, DIMENSION(kdlon,9), INTENT(out) :: topsw_aero
     126  REAL, DIMENSION(kdlon,9), INTENT(out) :: topsw0_aero
     127  REAL, DIMENSION(kdlon,9), INTENT(out) :: solsw_aero
     128  REAL, DIMENSION(kdlon,9), INTENT(out) :: solsw0_aero
     129
     130! Local variables
     131  REAL*8 ZFSUP(KDLON,KFLEV+1)
     132  REAL*8 ZFSDN(KDLON,KFLEV+1)
     133  REAL*8 ZFSUP0(KDLON,KFLEV+1)
     134  REAL*8 ZFSDN0(KDLON,KFLEV+1)
     135  REAL*8 ZFLUP(KDLON,KFLEV+1)
     136  REAL*8 ZFLDN(KDLON,KFLEV+1)
     137  REAL*8 ZFLUP0(KDLON,KFLEV+1)
     138  REAL*8 ZFLDN0(KDLON,KFLEV+1)
     139  REAL*8 zx_alpha1, zx_alpha2
     140  INTEGER k, kk, i, j, iof, nb_gr
     141  REAL*8 PSCT
     142  REAL*8 PALBD(kdlon,2), PALBP(kdlon,2)
     143  REAL*8 PEMIS(kdlon), PDT0(kdlon), PVIEW(kdlon)
     144  REAL*8 PPSOL(kdlon), PDP(kdlon,KLEV)
     145  REAL*8 PTL(kdlon,kflev+1), PPMB(kdlon,kflev+1)
     146  REAL*8 PTAVE(kdlon,kflev)
     147  REAL*8 PWV(kdlon,kflev), PQS(kdlon,kflev), POZON(kdlon,kflev)
     148  REAL*8 PAER(kdlon,kflev,5)
     149  REAL*8 PCLDLD(kdlon,kflev)
     150  REAL*8 PCLDLU(kdlon,kflev)
     151  REAL*8 PCLDSW(kdlon,kflev)
     152  REAL*8 PTAU(kdlon,2,kflev)
     153  REAL*8 POMEGA(kdlon,2,kflev)
     154  REAL*8 PCG(kdlon,2,kflev)
     155  REAL*8 zfract(kdlon), zrmu0(kdlon), zdist
     156  REAL*8 zheat(kdlon,kflev), zcool(kdlon,kflev)
     157  REAL*8 zheat0(kdlon,kflev), zcool0(kdlon,kflev)
     158  REAL*8 ztopsw(kdlon), ztoplw(kdlon)
     159  REAL*8 zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon)
     160  REAL*8 zsollwdown(kdlon)
     161  REAL*8 ztopsw0(kdlon), ztoplw0(kdlon)
     162  REAL*8 zsolsw0(kdlon), zsollw0(kdlon)
     163  REAL*8 zznormcp
     164  REAL*8 tauaero(kdlon,kflev,9,2)                     ! aer opt properties
     165  REAL*8 pizaero(kdlon,kflev,9,2)
     166  REAL*8 cgaero(kdlon,kflev,9,2)
     167  REAL*8 PTAUA(kdlon,2,kflev)                         ! present-day value of cloud opt thickness (PTAU is pre-industrial value), local use
     168  REAL*8 POMEGAA(kdlon,2,kflev)                       ! dito for single scatt albedo
     169  REAL*8 ztopswadaero(kdlon), zsolswadaero(kdlon)     ! Aerosol direct forcing at TOAand surface
     170  REAL*8 ztopswad0aero(kdlon), zsolswad0aero(kdlon)   ! Aerosol direct forcing at TOAand surface
     171  REAL*8 ztopswaiaero(kdlon), zsolswaiaero(kdlon)     ! dito, indirect
     172  REAL*8 ztopsw_aero(kdlon,9), ztopsw0_aero(kdlon,9)
     173  REAL*8 zsolsw_aero(kdlon,9), zsolsw0_aero(kdlon,9)
     174
     175  ! initialisation
     176  tauaero(:,:,:,:)=0.
     177  pizaero(:,:,:,:)=0.
     178  cgaero(:,:,:,:)=0.
     179 
     180  !
     181  !-------------------------------------------
     182  nb_gr = KLON / kdlon
     183  IF (nb_gr*kdlon .NE. KLON) THEN
     184      PRINT*, "kdlon mauvais:", KLON, kdlon, nb_gr
     185      CALL abort
     186  ENDIF
     187  IF (kflev .NE. KLEV) THEN
     188      PRINT*, "kflev differe de KLEV, kflev, KLEV"
     189      CALL abort
     190  ENDIF
     191  !-------------------------------------------
     192  DO k = 1, KLEV
     193    DO i = 1, KLON
     194      heat(i,k)=0.
     195      cool(i,k)=0.
     196      heat0(i,k)=0.
     197      cool0(i,k)=0.
     198    ENDDO
     199  ENDDO
     200  !
     201  zdist = dist
     202  !
     203  PSCT = solaire/zdist/zdist
     204  DO j = 1, nb_gr
     205    iof = kdlon*(j-1)
     206    DO i = 1, kdlon
     207      zfract(i) = fract(iof+i)
     208      zrmu0(i) = rmu0(iof+i)
     209      PALBD(i,1) = albedo(iof+i)
     210      PALBD(i,2) = alblw(iof+i)
     211      PALBP(i,1) = albedo(iof+i)
     212      PALBP(i,2) = alblw(iof+i)
     213      PEMIS(i) = 1.0
     214      PVIEW(i) = 1.66
     215      PPSOL(i) = paprs(iof+i,1)
     216      zx_alpha1 = (paprs(iof+i,1)-pplay(iof+i,2))/(pplay(iof+i,1)-pplay(iof+i,2))
     217      zx_alpha2 = 1.0 - zx_alpha1
     218      PTL(i,1) = t(iof+i,1) * zx_alpha1 + t(iof+i,2) * zx_alpha2
     219      PTL(i,KLEV+1) = t(iof+i,KLEV)
     220      PDT0(i) = tsol(iof+i) - PTL(i,1)
     221    ENDDO
     222    DO k = 2, kflev
    210223      DO i = 1, kdlon
    211          zfract(i) = fract(iof+i)
    212          zrmu0(i) = rmu0(iof+i)
    213          PALBD(i,1) = alb1(iof+i)
    214 !         PALBD(i,2) = alb1(iof+i)
    215          PALBD(i,2) = alb2(iof+i)
    216          PALBP(i,1) = alb1(iof+i)
    217 !         PALBP(i,2) = alb1(iof+i)
    218          PALBP(i,2) = alb2(iof+i)
    219 cIM cf. JLD pour etre en accord avec ORCHIDEE il faut mettre PEMIS(i) = 0.96
    220          PEMIS(i) = 1.0
    221          PVIEW(i) = 1.66
    222          PPSOL(i) = paprs(iof+i,1)
    223          zx_alpha1 = (paprs(iof+i,1)-pplay(iof+i,2))
    224      .             / (pplay(iof+i,1)-pplay(iof+i,2))
    225          zx_alpha2 = 1.0 - zx_alpha1
    226          PTL(i,1) = t(iof+i,1) * zx_alpha1 + t(iof+i,2) * zx_alpha2
    227          PTL(i,klev+1) = t(iof+i,klev)
    228          PDT0(i) = tsol(iof+i) - PTL(i,1)
    229       ENDDO
    230       DO k = 2, kflev
     224        PTL(i,k) = (t(iof+i,k)+t(iof+i,k-1))*0.5
     225      ENDDO
     226    ENDDO
     227    DO k = 1, kflev
    231228      DO i = 1, kdlon
    232          PTL(i,k) = (t(iof+i,k)+t(iof+i,k-1))*0.5
    233       ENDDO
    234       ENDDO
     229        PDP(i,k) = paprs(iof+i,k)-paprs(iof+i,k+1)
     230        PTAVE(i,k) = t(iof+i,k)
     231        PWV(i,k) = MAX (q(iof+i,k), 1.0e-12)
     232        PQS(i,k) = PWV(i,k)
     233        ! wo:    cm.atm (epaisseur en cm dans la situation standard)
     234        ! POZON: kg/kg
     235        POZON(i,k) = MAX(wo(iof+i,k),1.0e-12)*RG/46.6968 &
     236           /(paprs(iof+i,k)-paprs(iof+i,k+1))&
     237           *(paprs(iof+i,1)/101325.0)
     238        PCLDLD(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
     239        PCLDLU(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
     240        PCLDSW(i,k) = cldfra(iof+i,k)
     241        PTAU(i,1,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! 1e-12 serait instable
     242        PTAU(i,2,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! pour 32-bit machines
     243        POMEGA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAU(i,1,k))
     244        POMEGA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i,2,k))
     245        PCG(i,1,k) = 0.865
     246        PCG(i,2,k) = 0.910
     247        !-
     248        ! Introduced for aerosol indirect forcings.
     249        ! The following values use the cloud optical thickness calculated from
     250        ! present-day aerosol concentrations whereas the quantities without the
     251        ! "A" at the end are for pre-industial (natural-only) aerosol concentrations
     252        !
     253        PTAUA(i,1,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! 1e-12 serait instable
     254        PTAUA(i,2,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! pour 32-bit machines
     255        POMEGAA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAUA(i,1,k))
     256        POMEGAA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i,2,k))
     257      ENDDO
     258    ENDDO
     259    !
     260    DO k = 1, kflev+1
     261      DO i = 1, kdlon
     262        PPMB(i,k) = paprs(iof+i,k)/100.0
     263      ENDDO
     264    ENDDO
     265    !
     266    DO kk = 1, 5
    235267      DO k = 1, kflev
     268        DO i = 1, kdlon
     269          PAER(i,k,kk) = 1.0E-15
     270        ENDDO
     271      ENDDO
     272    ENDDO
     273    DO k = 1, kflev
    236274      DO i = 1, kdlon
    237          PDP(i,k) = paprs(iof+i,k)-paprs(iof+i,k+1)
    238          PTAVE(i,k) = t(iof+i,k)
    239          PWV(i,k) = MAX (q(iof+i,k), 1.0e-12)
    240          PQS(i,k) = PWV(i,k)
    241 c wo:    cm.atm (epaisseur en cm dans la situation standard)
    242 c POZON: kg/kg
    243          POZON(i,k) = MAX(wo(iof+i,k),1.0e-12)*RG/46.6968
    244      .               /(paprs(iof+i,k)-paprs(iof+i,k+1))
    245      .               *(paprs(iof+i,1)/101325.0)
    246          PCLDLD(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
    247          PCLDLU(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
    248          PCLDSW(i,k) = cldfra(iof+i,k)
    249          PTAU(i,1,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! 1e-12 serait instable
    250          PTAU(i,2,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! pour 32-bit machines
    251          POMEGA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAU(i,1,k))
    252          POMEGA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i,2,k))
    253          PCG(i,1,k) = 0.865
    254          PCG(i,2,k) = 0.910
    255 c-OB
    256 cjq Introduced for aerosol indirect forcings.
    257 cjq The following values use the cloud optical thickness calculated from
    258 cjq present-day aerosol concentrations whereas the quantities without the
    259 cjq "A" at the end are for pre-industial (natural-only) aerosol concentrations
    260 cjq
    261          PTAUA(i,1,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! 1e-12 serait instable
    262          PTAUA(i,2,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! pour 32-bit machines
    263          POMEGAA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAUA(i,1,k))
    264          POMEGAA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i,2,k))
    265 cjq-end
    266       ENDDO
    267       ENDDO
    268 c
     275        tauaero(i,k,:,1)=tau_aero(iof+i,k,:,1)
     276        pizaero(i,k,:,1)=piz_aero(iof+i,k,:,1)
     277        cgaero(i,k,:,1) =cg_aero(iof+i,k,:,1)
     278        tauaero(i,k,:,2)=tau_aero(iof+i,k,:,2)
     279        pizaero(i,k,:,2)=piz_aero(iof+i,k,:,2)
     280        cgaero(i,k,:,2) =cg_aero(iof+i,k,:,2)
     281      ENDDO
     282    ENDDO
     283    !
     284    !======================================================================
     285    CALL LW_LMDAR4(&
     286       PPMB, PDP,&
     287       PPSOL,PDT0,PEMIS,&
     288       PTL, PTAVE, PWV, POZON, PAER,&
     289       PCLDLD,PCLDLU,&
     290       PVIEW,&
     291       zcool, zcool0,&
     292       ztoplw,zsollw,ztoplw0,zsollw0,&
     293       zsollwdown,&
     294       ZFLUP, ZFLDN, ZFLUP0,ZFLDN0)
     295
     296
     297    IF (.NOT. new_aod) THEN
     298       ! use old version
     299        CALL SW_LMDAR4(PSCT, zrmu0, zfract,&
     300           PPMB, PDP, &
     301           PPSOL, PALBD, PALBP,&
     302           PTAVE, PWV, PQS, POZON, PAER,&
     303           PCLDSW, PTAU, POMEGA, PCG,&
     304           zheat, zheat0,&
     305           zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,&
     306           ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,&
     307           tau_aero(:,:,5,:), piz_aero(:,:,5,:), cg_aero(:,:,5,:),&
     308           PTAUA, POMEGAA,&
     309           ztopswadaero,zsolswadaero,&
     310           ztopswaiaero,zsolswaiaero,&
     311           ok_ade, ok_aie)
     312    ELSE
     313
     314        CALL SW_AERO(PSCT, zrmu0, zfract,&
     315           PPMB, PDP,&
     316           PPSOL, PALBD, PALBP,&
     317           PTAVE, PWV, PQS, POZON, PAER,&
     318           PCLDSW, PTAU, POMEGA, PCG,&
     319           zheat, zheat0,&
     320           zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,&
     321           ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,&
     322           tauaero, pizaero, cgaero, &
     323           PTAUA, POMEGAA,&
     324           ztopswadaero,zsolswadaero,&
     325           ztopswad0aero,zsolswad0aero,&
     326           ztopswaiaero,zsolswaiaero, &
     327           ztopsw_aero,ztopsw0_aero,&
     328           zsolsw_aero,zsolsw0_aero,&
     329           ok_ade, ok_aie)
     330   
     331    ENDIF
     332    !======================================================================
     333    DO i = 1, kdlon
     334      radsol(iof+i) = zsolsw(i) + zsollw(i)
     335      topsw(iof+i) = ztopsw(i)
     336      toplw(iof+i) = ztoplw(i)
     337      solsw(iof+i) = zsolsw(i)
     338      sollw(iof+i) = zsollw(i)
     339      sollwdown(iof+i) = zsollwdown(i)
    269340      DO k = 1, kflev+1
     341        lwdn0 ( iof+i,k)   = ZFLDN0 ( i,k)
     342        lwdn  ( iof+i,k)   = ZFLDN  ( i,k)
     343        lwup0 ( iof+i,k)   = ZFLUP0 ( i,k)
     344        lwup  ( iof+i,k)   = ZFLUP  ( i,k)
     345      ENDDO
     346      topsw0(iof+i) = ztopsw0(i)
     347      toplw0(iof+i) = ztoplw0(i)
     348      solsw0(iof+i) = zsolsw0(i)
     349      sollw0(iof+i) = zsollw0(i)
     350      albpla(iof+i) = zalbpla(i)
     351
     352      DO k = 1, kflev+1
     353        swdn0 ( iof+i,k)   = ZFSDN0 ( i,k)
     354        swdn  ( iof+i,k)   = ZFSDN  ( i,k)
     355        swup0 ( iof+i,k)   = ZFSUP0 ( i,k)
     356        swup  ( iof+i,k)   = ZFSUP  ( i,k)
     357      ENDDO
     358    ENDDO
     359    !-transform the aerosol forcings, if they have
     360    ! to be calculated
     361    IF (ok_ade) THEN
     362        DO i = 1, kdlon
     363          topswad_aero(iof+i) = ztopswadaero(i)
     364          topswad0_aero(iof+i) = ztopswad0aero(i)
     365          solswad_aero(iof+i) = zsolswadaero(i)
     366          solswad0_aero(iof+i) = zsolswad0aero(i)
     367          topsw_aero(iof+i,:) = ztopsw_aero(iof+i,:)
     368          topsw0_aero(iof+i,:) = ztopsw0_aero(iof+i,:)
     369          solsw_aero(iof+i,:) = zsolsw_aero(iof+i,:)
     370          solsw0_aero(iof+i,:) = zsolsw0_aero(iof+i,:)
     371         
     372        ENDDO
     373    ELSE
     374        DO i = 1, kdlon
     375          topswad_aero(iof+i) = 0.0
     376          solswad_aero(iof+i) = 0.0
     377          topswad0_aero(iof+i) = 0.0
     378          solswad0_aero(iof+i) = 0.0
     379          topsw_aero(iof+i,:) = 0.
     380          topsw0_aero(iof+i,:) =0.
     381          solsw_aero(iof+i,:) = 0.
     382          solsw0_aero(iof+i,:) = 0.
     383        ENDDO
     384    ENDIF
     385    IF (ok_aie) THEN
     386        DO i = 1, kdlon
     387          topswai_aero(iof+i) = ztopswaiaero(i)
     388          solswai_aero(iof+i) = zsolswaiaero(i)
     389        ENDDO
     390    ELSE
     391        DO i = 1, kdlon
     392          topswai_aero(iof+i) = 0.0
     393          solswai_aero(iof+i) = 0.0
     394        ENDDO
     395    ENDIF
     396    DO k = 1, kflev
    270397      DO i = 1, kdlon
    271          PPMB(i,k) = paprs(iof+i,k)/100.0
    272       ENDDO
    273       ENDDO
    274 c
    275       DO kk = 1, 5
    276       DO k = 1, kflev
    277       DO i = 1, kdlon
    278          PAER(i,k,kk) = 1.0E-15
    279       ENDDO
    280       ENDDO
    281       ENDDO
    282 c-OB
    283       DO k = 1, kflev
    284       DO i = 1, kdlon
    285         tauae(i,k,1)=tau_ae(iof+i,k,1)
    286         pizae(i,k,1)=piz_ae(iof+i,k,1)
    287         cgae(i,k,1) =cg_ae(iof+i,k,1)
    288         tauae(i,k,2)=tau_ae(iof+i,k,2)
    289         pizae(i,k,2)=piz_ae(iof+i,k,2)
    290         cgae(i,k,2) =cg_ae(iof+i,k,2)
    291       ENDDO
    292       ENDDO
    293 c
    294 c===== si iflag_rrtm=0 ================================================
    295 cIM ctes ds clesphys.h   CALL LW(RCO2,RCH4,RN2O,RCFC11,RCFC12,
    296 cIM ctes ds clesphys.h   CALL SW(PSCT, RCO2, zrmu0, zfract,
    297 c     
    298       if (iflag_rrtm.eq.0) then
    299          CALL LW_LMDAR4(
    300      .        PPMB, PDP,
    301      .        PPSOL,PDT0,PEMIS,
    302      .        PTL, PTAVE, PWV, POZON, PAER,
    303      .        PCLDLD,PCLDLU,
    304      .        PVIEW,
    305      .        zcool, zcool0,
    306      .        ztoplw,zsollw,ztoplw0,zsollw0,
    307      .        zsollwdown,
    308      .        ZFLUP, ZFLDN, ZFLUP0,ZFLDN0)
    309          CALL SW_LMDAR4(PSCT, zrmu0, zfract,
    310      S        PPMB, PDP,
    311      S        PPSOL, PALBD, PALBP,
    312      S        PTAVE, PWV, PQS, POZON, PAER,
    313      S        PCLDSW, PTAU, POMEGA, PCG,
    314      S        zheat, zheat0,
    315      S        zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,
    316      S        ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,
    317      S        tauae, pizae, cgae, ! aerosol optical properties
    318      s        PTAUA, POMEGAA,
    319      s        ztopswad,zsolswad,ztopswai,zsolswai, ! diagnosed aerosol forcing
    320      J        ok_ade, ok_aie) ! apply aerosol effects or not?
    321       else
    322 c===== si iflag_rrtm=1, on passe dans SW via RECMWFL ===============
    323           PRINT*, "Cette option ne fonctionne pas encore !!!"
    324          CALL abort
    325          endif   ! if(iflag_rrtm=0)
    326 
    327 c======================================================================
    328       DO i = 1, kdlon
    329          radsol(iof+i) = zsolsw(i) + zsollw(i)
    330          topsw(iof+i) = ztopsw(i)
    331          toplw(iof+i) = ztoplw(i)
    332          solsw(iof+i) = zsolsw(i)
    333          sollw(iof+i) = zsollw(i)
    334          sollwdown(iof+i) = zsollwdown(i)
    335 cIM
    336          DO k = 1, kflev+1
    337          lwdn0 ( iof+i,k)   = ZFLDN0 ( i,k)
    338          lwdn  ( iof+i,k)   = ZFLDN  ( i,k)
    339          lwup0 ( iof+i,k)   = ZFLUP0 ( i,k)
    340          lwup  ( iof+i,k)   = ZFLUP  ( i,k)
    341          ENDDO
    342 c
    343          topsw0(iof+i) = ztopsw0(i)
    344          toplw0(iof+i) = ztoplw0(i)
    345          solsw0(iof+i) = zsolsw0(i)
    346          sollw0(iof+i) = zsollw0(i)
    347          albpla(iof+i) = zalbpla(i)
    348 cIM
    349          DO k = 1, kflev+1
    350          swdn0 ( iof+i,k)   = ZFSDN0 ( i,k)
    351          swdn  ( iof+i,k)   = ZFSDN  ( i,k)
    352          swup0 ( iof+i,k)   = ZFSUP0 ( i,k)
    353          swup  ( iof+i,k)   = ZFSUP  ( i,k)
    354          ENDDO !k=1, kflev+1
    355       ENDDO
    356 cjq-transform the aerosol forcings, if they have
    357 cjq to be calculated
    358       IF (ok_ade) THEN
    359       DO i = 1, kdlon
    360          topswad(iof+i) = ztopswad(i)
    361          solswad(iof+i) = zsolswad(i)
    362       ENDDO
    363       ELSE
    364       DO i = 1, kdlon
    365          topswad(iof+i) = 0.0
    366          solswad(iof+i) = 0.0
    367       ENDDO
    368       ENDIF
    369       IF (ok_aie) THEN
    370       DO i = 1, kdlon
    371          topswai(iof+i) = ztopswai(i)
    372          solswai(iof+i) = zsolswai(i)
    373       ENDDO
    374       ELSE
    375       DO i = 1, kdlon
    376          topswai(iof+i) = 0.0
    377          solswai(iof+i) = 0.0
    378       ENDDO
    379       ENDIF
    380 cjq-end
    381       DO k = 1, kflev
    382 c      DO i = 1, kdlon
    383 c         heat(iof+i,k) = zheat(i,k)
    384 c         cool(iof+i,k) = zcool(i,k)
    385 c         heat0(iof+i,k) = zheat0(i,k)
    386 c         cool0(iof+i,k) = zcool0(i,k)
    387 c      ENDDO
    388       DO i = 1, kdlon
    389 C        scale factor to take into account the difference between
    390 C        dry air and watter vapour scpecific heat capacity
    391          zznormcp=1.0+RVTMP2*PWV(i,k)
    392          heat(iof+i,k) = zheat(i,k)/zznormcp
    393          cool(iof+i,k) = zcool(i,k)/zznormcp
    394          heat0(iof+i,k) = zheat0(i,k)/zznormcp
    395          cool0(iof+i,k) = zcool0(i,k)/zznormcp
    396       ENDDO
    397       ENDDO
    398 c
    399 99999 CONTINUE
    400       RETURN
    401       END
     398        !        scale factor to take into account the difference between
     399        !        dry air and watter vapour scpecifi! heat capacity
     400        zznormcp=1.0+RVTMP2*PWV(i,k)
     401        heat(iof+i,k) = zheat(i,k)/zznormcp
     402        cool(iof+i,k) = zcool(i,k)/zznormcp
     403        heat0(iof+i,k) = zheat0(i,k)/zznormcp
     404        cool0(iof+i,k) = zcool0(i,k)/zznormcp
     405      ENDDO
     406    ENDDO
     407    !
     408  ENDDO
     409
     410
     411ENDSUBROUTINE radlwsw_aero
     412
     413
Note: See TracChangeset for help on using the changeset viewer.