Ignore:
Timestamp:
Apr 28, 2021, 4:55:57 PM (4 years ago)
Author:
idelkadi
Message:

Online implementation of the radiative transfer code ECRAD in LMDZ.

  • Inclusion of the ecrad directory containing the sources of the ECRAD code
  • Adaptation of compilation scripts (CPP_ECRAD keys)
  • Call of ecrad in radlwsw_m.F90 under the logical key iflag_rrtm = 2
Location:
LMDZ6/branches/LMDZ-ECRAD/libf/phylmd
Files:
369 added
3 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/LMDZ-ECRAD/libf/phylmd/conf_phys_m.F90

    r3792 r3880  
    25162516      IF (NSW.NE.2.AND.NSW.NE.4.AND.NSW.NE.6) THEN
    25172517        WRITE(lunout,*) ' ERROR iflag_rrtm=1 and NSW<>2,4,6 not possible'
     2518        CALL abort_physic('conf_phys','choice NSW not valid',1)
     2519      ENDIF
     2520    ELSE IF (iflag_rrtm .EQ. 2) THEN
     2521      IF (NSW.NE.2.AND.NSW.NE.4.AND.NSW.NE.6) THEN
     2522        WRITE(lunout,*) ' ERROR iflag_rrtm=2 and NSW<>2,4,6 not possible'
    25182523        CALL abort_physic('conf_phys','choice NSW not valid',1)
    25192524      ENDIF

  • LMDZ6/branches/LMDZ-ECRAD/libf/phylmd/physiq_mod.F90

    r3792 r3880  
    39943994
    39953995       IF (ok_newmicro) then
    3996           IF (iflag_rrtm.NE.0) THEN
     3996!          IF (iflag_rrtm.NE.0) THEN
     3997          IF (iflag_rrtm.EQ.1) THEN
    39973998#ifdef CPP_RRTM
    39983999             IF (ok_cdnc.AND.NRADLP.NE.3) THEN
    … …  
    41144115       !input to radiation (DICE)
    41154116       !
    4116        IF (iflag_radia .ge. 2) THEN
     4117       IF (iflag_radia .gt. 2) THEN
    41174118          zsav_tsol (:) = zxtsol(:)
    41184119          CALL perturb_radlwsw(zxtsol,iflag_radia)

  • LMDZ6/branches/LMDZ-ECRAD/libf/phylmd/radlwsw_m.F90

    r3756 r3880  
    1919   flag_aerosol_strat, flag_aer_feedback, &
    2020   tau_aero, piz_aero, cg_aero,&
    21    tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,& ! rajoute par OB pour RRTM
    22    tau_aero_lw_rrtm, &                                   ! rajoute par C. Kleinschmitt pour RRTM
     21   tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,& ! rajoute par OB RRTM
     22   tau_aero_lw_rrtm, &              ! rajoute par C.Kleinschmitt pour RRTM
    2323   cldtaupi, &
    2424   qsat, flwc, fiwc, &
    … …  
    4545   ZSWFT0_i, ZFSDN0, ZFSUP0)
    4646
    47 
    48 
     47! Modules necessaires
    4948  USE DIMPHY
    5049  USE assert_m, ONLY : assert
    5150  USE infotrac_phy, ONLY : type_trac
    5251  USE write_field_phy
     52
    5353#ifdef REPROBUS
    5454  USE CHEM_REP, ONLY : solaireTIME, ok_SUNTIME, ndimozon
    5555#endif
     56
    5657#ifdef CPP_RRTM
    5758!    modules necessaires au rayonnement
    5859!    -----------------------------------------
    59 !     USE YOMCST   , ONLY : RG       ,RD       ,RTT      ,RPI
    60 !     USE YOERAD   , ONLY : NSW      ,LRRTM    ,LINHOM   , LCCNL,LCCNO,
    61 !     USE YOERAD   , ONLY : NSW      ,LRRTM    ,LCCNL    ,LCCNO ,&
    62 ! NSW mis dans .def MPL 20140211
    63 ! NLW ajoute par OB
    6460      USE YOERAD   , ONLY : NLW, LRRTM    ,LCCNL    ,LCCNO ,&
    6561          NRADIP   , NRADLP , NICEOPT, NLIQOPT ,RCCNLND  , RCCNSEA
    … …  
    7369          RFLDD1   ,RFLDD2   ,RFLDD3   ,RFUETA   ,RASWCA,&
    7470          RASWCB   ,RASWCC   ,RASWCD   ,RASWCE   ,RASWCF
    75 !    &    RASWCB   ,RASWCC   ,RASWCD   ,RASWCE   ,RASWCF, RLINLI
    7671      USE YOERDU   , ONLY : NUAER  ,NTRAER ,REPLOG ,REPSC  ,REPSCW ,DIFF
    77 !      USE YOETHF   , ONLY : RTICE
    7872      USE YOERRTWN , ONLY : DELWAVE   ,TOTPLNK     
    7973      USE YOMPHY3  , ONLY : RII0
    … …  
    8175      USE aero_mod
    8276
     77! AI 02.2021
     78! Besoin pour ECRAD de pctsrf, zmasq, longitude, altitude
     79#ifdef CPP_ECRAD
     80      USE geometry_mod, ONLY: latitude, longitude
     81      USE phys_state_var_mod, ONLY: pctsrf
     82      USE indice_sol_mod
     83      USE time_phylmdz_mod, only: current_time
     84      USE phys_cal_mod, only: day_cur
     85#endif
     86
    8387  !======================================================================
    8488  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
    8589  ! Objet: interface entre le modele et les rayonnements
    8690  ! Arguments:
    87   ! dist-----input-R- distance astronomique terre-soleil
    88   ! rmu0-----input-R- cosinus de l'angle zenithal
    89   ! fract----input-R- duree d'ensoleillement normalisee
    90   ! co2_ppm--input-R- concentration du gaz carbonique (en ppm)
    91   ! paprs----input-R- pression a inter-couche (Pa)
    92   ! pplay----input-R- pression au milieu de couche (Pa)
    93   ! tsol-----input-R- temperature du sol (en K)
    94   ! alb1-----input-R- albedo du sol(entre 0 et 1) dans l'interval visible
    95   ! alb2-----input-R- albedo du sol(entre 0 et 1) dans l'interval proche infra-rouge   
    96   ! t--------input-R- temperature (K)
    97   ! q--------input-R- vapeur d'eau (en kg/kg)
    98   ! cldfra---input-R- fraction nuageuse (entre 0 et 1)
    99   ! cldtaupd---input-R- epaisseur optique des nuages dans le visible (present-day value)
    100   ! cldemi---input-R- emissivite des nuages dans l'IR (entre 0 et 1)
    101   ! ok_ade---input-L- apply the Aerosol Direct Effect or not?
    102   ! ok_aie---input-L- apply the Aerosol Indirect Effect or not?
    103   ! ok_volcan-input-L- activate volcanic diags (SW heat & LW cool rate, SW & LW flux)
    104   ! flag_aerosol-input-I- aerosol flag from 0 to 6
    105   ! flag_aerosol_strat-input-I- use stratospheric aerosols flag (0, 1, 2)
    106   ! flag_aer_feedback-input-I- activate aerosol radiative feedback (T, F)
    107   ! tau_ae, piz_ae, cg_ae-input-R- aerosol optical properties (calculated in aeropt.F)
    108   ! cldtaupi-input-R- epaisseur optique des nuages dans le visible
     91  !                  INPUTS
     92  ! dist----- input-R- distance astronomique terre-soleil
     93  ! rmu0----- input-R- cosinus de l'angle zenithal
     94  ! fract---- input-R- duree d'ensoleillement normalisee
     95  ! co2_ppm-- input-R- concentration du gaz carbonique (en ppm)
     96  ! paprs---- input-R- pression a inter-couche (Pa)
     97  ! pplay---- input-R- pression au milieu de couche (Pa)
     98  ! tsol----- input-R- temperature du sol (en K)
     99  ! alb1----- input-R- albedo du sol(entre 0 et 1) dans l'interval visible
     100  ! alb2----- input-R- albedo du sol(entre 0 et 1) dans l'interval proche infra-rouge   
     101  ! t-------- input-R- temperature (K)
     102  ! q-------- input-R- vapeur d'eau (en kg/kg)
     103  ! cldfra--- input-R- fraction nuageuse (entre 0 et 1)
     104  ! cldtaupd- input-R- epaisseur optique des nuages dans le visible (present-day value)
     105  ! cldemi--- input-R- emissivite des nuages dans l'IR (entre 0 et 1)
     106  ! ok_ade--- input-L- apply the Aerosol Direct Effect or not?
     107  ! ok_aie--- input-L- apply the Aerosol Indirect Effect or not?
     108  ! ok_volcan input-L- activate volcanic diags (SW heat & LW cool rate, SW & LW flux)
     109  ! flag_aerosol input-I- aerosol flag from 0 to 6
     110  ! flag_aerosol_strat input-I- use stratospheric aerosols flag (0, 1, 2)
     111  ! flag_aer_feedback  input-I- activate aerosol radiative feedback (T, F)
     112  ! tau_ae, piz_ae, cg_ae input-R- aerosol optical properties (calculated in aeropt.F)
     113  ! cldtaupi  input-R- epaisseur optique des nuages dans le visible
    109114  !                   calculated for pre-industrial (pi) aerosol concentrations, i.e. with smaller
    110115  !                   droplet concentration, thus larger droplets, thus generally cdltaupi cldtaupd
    111116  !                   it is needed for the diagnostics of the aerosol indirect radiative forcing     
    112117  !
     118  !                  OUTPUTS
    113119  ! heat-----output-R- echauffement atmospherique (visible) (K/jour)
    114120  ! cool-----output-R- refroidissement dans l'IR (K/jour)
    … …  
    177183  !
    178184  ! ====================================================================
     185
     186! ==============
     187! DECLARATIONS
     188! ==============
    179189  include "YOETHF.h"
    180190  include "YOMCST.h"
    … …  
    286296  REAL(KIND=8) PTAVE(kdlon,kflev)
    287297  REAL(KIND=8) PWV(kdlon,kflev), PQS(kdlon,kflev)
     298
     299!!!!!!! Declarations specifiques pour ECRAD !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     300! AI 02.2021
     301#ifdef CPP_ECRAD
     302! ATTENTION les dimensions klon, kdlon ???
     303! INPUTS
     304  REAL(KIND=8) ZEMISW(klon), &              ! LW emissivity inside the window region
     305               ZEMIS(klon)                  ! LW emissivity outside the window region
     306  REAL(KIND=8) ZGELAM(klon), &              ! longitudes en rad
     307               ZGEMU(klon)                  ! sin(latitude)
     308  REAL(KIND=8) ZCO2(klon,klev), &           ! CO2 mass mixing ratios on full levels
     309               ZCH4(klon,klev), &           ! CH4 mass mixing ratios on full levels
     310               ZN2O(klon,klev), &           ! N2O mass mixing ratios on full levels
     311               ZNO2(klon,klev), &           ! NO2 mass mixing ratios on full levels
     312               ZCFC11(klon,klev), &         ! CFC11
     313               ZCFC12(klon,klev), &         ! CFC12
     314               ZHCFC22(klon,klev), &        ! HCFC22
     315               ZCCL4(klon,klev), &          ! CCL4
     316               ZO3_DP(klon,klev), ZO3_DP_i(klon,klev)            ! Ozone
     317  REAL(KIND=8) ZQ_RAIN(klon,klev), &        ! Rain cloud mass mixing ratio (kg/kg) ?
     318               ZQ_SNOW(klon,klev)           ! Snow cloud mass mixing ratio (kg/kg) ?
     319  REAL(KIND=8) ZAEROSOL_OLD(KLON,6,KLEV), &  !
     320               ZAEROSOL(KLON,KLEV,naero_tot) !
     321! OUTPUTS
     322  REAL(KIND=8) ZFLUX_DIR(klon), &           ! Direct compt of surf flux into horizontal plane
     323               ZFLUX_DIR_CLEAR(klon), &     ! CS Direct
     324               ZFLUX_DIR_INTO_SUN(klon), &  !
     325               ZFLUX_UV(klon), &            ! UV flux
     326               ZFLUX_PAR(klon), &           ! photosynthetically active radiation similarly
     327               ZFLUX_PAR_CLEAR(klon), &     ! CS photosynthetically
     328               ZFLUX_SW_DN_TOA(klon), &     ! DN SW flux at TOA
     329               ZEMIS_OUT(klon)              ! effective broadband emissivity
     330  REAL(KIND=8) ZLWDERIVATIVE(klon,klev+1)   ! LW derivatives
     331  REAL(KIND=8) ZSWDIFFUSEBAND(klon,NSW), &  ! SW DN flux in diffuse albedo band
     332               ZSWDIRECTBAND(klon,NSW)      ! SW DN flux in direct albedo band
     333#endif
     334!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    288335
    289336  REAL(kind=8) POZON(kdlon, kflev, size(wo, 3)) ! mass fraction of ozone
    … …  
    401448  REAL zdir, zdif
    402449
     450! =========  INITIALISATIONS ==============================================
     451  print*,'Entree dans radlwsw '
     452  print*,'************* INITIALISATIONS *****************************'
     453  print*,'klon, kdlon, klev, kflev =',klon, kdlon, klev, kflev
     454
    403455  CALL assert(size(wo, 1) == klon, size(wo, 2) == klev, "radlwsw wo")
    404   ! initialisation
     456 
    405457  ist=1
    406458  iend=klon
    407459  ktdia=1
    408460  kmode=ist
     461! Aeros
    409462  tauaero(:,:,:,:)=0.
    410463  pizaero(:,:,:,:)=0.
    … …  
    462515  ENDIF
    463516
     517  print*,'************** Debut boucle de 1 a ', nb_gr
    464518  DO j = 1, nb_gr
    465519    iof = kdlon*(j-1)
    466520    DO i = 1, kdlon
    467521      zfract(i) = fract(iof+i)
    468 !     zfract(i) = 1.     !!!!!!  essai MPL 19052010
    469522      zrmu0(i) = rmu0(iof+i)
    470523
    471524
    472 !albedo SB >>>
    473 !
    474525      IF (iflag_rrtm==0) THEN
    475 !
     526!     Albedo
    476527        PALBD(i,1)=alb_dif(iof+i,1)
    477528        PALBD(i,2)=alb_dif(iof+i,2)
    478529        PALBP(i,1)=alb_dir(iof+i,1)
    479530        PALBP(i,2)=alb_dir(iof+i,2)
    480 !
    481       ELSEIF (iflag_rrtm==1) THEn
    482 !
     531! AI 02.2021 cas iflag_rrtm=1 et 2
     532       ELSEIF (iflag_rrtm==1.OR.iflag_rrtm==2) THEN
    483533        DO kk=1,NSW
    484534          PALBD_NEW(i,kk)=alb_dif(iof+i,kk)
    … …  
    488538      ENDIF
    489539!albedo SB <<<
    490 
    491540
    492541      PEMIS(i) = 1.0    !!!!! A REVOIR (MPL)
    … …  
    513562        POZON(i,k, :) = wo(iof+i, k, :) * RG * dobson_u * 1e3 &
    514563             / (paprs(iof+i, k) - paprs(iof+i, k+1))
     564        ZO3_DP(i,k) = wo(iof+i, k, 1) * RG * dobson_u * 1e3
    515565!       A activer pour CCMVAL on prend l'ozone impose (MPL 07042010)
    516566!       POZON(i,k,:) = wo(i,k,:) 
    … …  
    569619      ENDDO
    570620    ENDDO
     621!
     622! AI 02.2021
     623#ifdef CPP_ECRAD
     624  ZEMIS = 1.0
     625  ZEMISW = 1.0
     626  ZGELAM = longitude
     627  ZGEMU = sin(latitude)
     628  ZCO2 = RCO2
     629  ZCH4 = RCH4
     630  ZN2O = RN2O
     631  ZNO2 = 0.0
     632  ZCFC11 = RCFC11
     633  ZCFC12 = RCFC12
     634  ZHCFC22 = 0.0
     635  ZCCL4 = 0.0
     636  ZQ_RAIN = 0.0
     637  ZQ_SNOW = 0.0
     638  ZAEROSOL_OLD = 0.0
     639  ZAEROSOL = 0.0
     640#endif
    571641!
    572642!===== iflag_rrtm ================================================
    … …  
    693763       ENDDO 
    694764!
    695     ELSE
     765    ELSE IF (iflag_rrtm == 1) then
    696766#ifdef CPP_RRTM
    697767!      if (prt_level.gt.10)write(lunout,*)'CPP_RRTM=.T.'
    … …  
    804874            ENDDO
    805875         ENDDO
     876
    806877!       print *,'RADLWSW: avant RECMWFL, RI0,rmu0=',solaire,rmu0
    807878
    … …  
    850921! Nouvel appel a RECMWF (celui du cy32t0)
    851922         CALL RECMWF_AERO (ist , iend, klon , ktdia  , klev   , kmode ,&
     923!                          KST,  KEND, KPROMA, KTDIA , KLEV,  KMODE
    852924         PALBD_NEW,PALBP_NEW, paprs_i , pplay_i , RCO2   , cldfra_i,&
     925!        PALBD ,   PALBP ,    PAPRS ,   PAPRSF ,  PCCO2 , PCLFR
    853926         POZON_i  , PAER_i  , PDP_i   , PEMIS   , rmu0   ,&
     927!        PQO3  ,    PAER  ,   PDP   ,   PEMIS  , PMU0
    854928          q_i     , qsat_i  , fiwc_i  , flwc_i  , zmasq  , t_i  ,tsol,&
     929!         PQ       , PQS   , PQIWP ,     PQLWP , PSLM   , PT    , PTS,
    855930         ref_liq_i, ref_ice_i, &
     931!        PREF_LIQ, PREF_ICE
    856932         ref_liq_pi_i, ref_ice_pi_i, &   ! rajoute par OB pour diagnostiquer effet indirect
     933!        PREF_LIQ_PI, PREF_ICE_PI
    857934         ZEMTD_i  , ZEMTU_i , ZTRSO_i ,&
     935!        PEMTD , PEMTU , PTRSO
    858936         ZTH_i    , ZCTRSO  , ZCEMTR  , ZTRSOD  ,&
     937!        PTH   ,     PCTRSO,   PCEMTR, PTRSOD
    859938         ZLWFC    , ZLWFT_i , ZSWFC   , ZSWFT_i ,&
     939!        PLWFC,     PLWFT,    PSWFC,    PSWFT,
    860940         PSFSWDIR , PSFSWDIF, PFSDNN  , PFSDNV  ,&
     941!        PSFSWDIR,  PSFSWDIF, PFSDNN,    PFSDNV
    861942         PPIZA_TOT, PCGA_TOT,PTAU_TOT,&
     943!        PPIZA_TOT, PCGA_TOT,PTAU_TOT
    862944         PPIZA_NAT, PCGA_NAT,PTAU_NAT,           &  ! rajoute par OB pour diagnostiquer effet direct
     945!        PPIZA_NAT,PCGA_NAT,PTAU_NAT
    863946         PTAU_LW_TOT, PTAU_LW_NAT,               &  ! rajoute par C. Kleinschmitt
     947!        PTAU_LW_TOT, PTAU_LW_NAT,
    864948         ZFLUX_i  , ZFLUC_i ,&
     949!        PFLUX,      PFLUC
    865950         ZFSDWN_i , ZFSUP_i , ZFCDWN_i, ZFCUP_i, ZFCCDWN_i, ZFCCUP_i, ZFLCCDWN_i, ZFLCCUP_i, &
     951!        PFSDN ,    PFSUP ,   PFSCDN , PFSCUP,   PFSCCDN,   PFSCCUP,  PFLCCDN,    PFLCCUP
    866952         ZTOPSWADAERO,ZSOLSWADAERO,&  ! rajoute par OB pour diagnostics
     953!        PTOPSWADAERO,PSOLSWADAERO
    867954         ZTOPSWAD0AERO,ZSOLSWAD0AERO,&
    868955         ZTOPSWAIAERO,ZSOLSWAIAERO, &
    … …  
    874961         ZLWADAERO, & !--NL
    875962         ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat, flag_aer_feedback) ! flags aerosols
     963
     964!--OB diagnostics
     965! & PTOPSWAIAERO,PSOLSWAIAERO,&
     966! & PTOPSWCFAERO,PSOLSWCFAERO,&
     967! & PSWADAERO,& !--NL
     968!!--LW diagnostics CK
     969! & PTOPLWADAERO,PSOLLWADAERO,&
     970! & PTOPLWAD0AERO,PSOLLWAD0AERO,&
     971! & PTOPLWAIAERO,PSOLLWAIAERO,&
     972! & PLWADAERO,& !--NL
     973!!..end
     974! & ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat,&
     975! & flag_aer_feedback)
     976
    876977           
    877978!        print *,'RADLWSW: apres RECMWF'
    … …  
    10611162    call abort_physic(modname, abort_message, 1)
    10621163#endif
    1063     ENDIF ! iflag_rrtm
     1164!======================================================================
     1165! AI fev 2021
     1166    ELSE IF(iflag_rrtm == 2) THEN
     1167    print*,'Traitement cas iflag_rrtm = ',iflag_rrtm
     1168    print*,'Mise a zero des flux '
     1169#ifdef CPP_ECRAD
     1170      DO k = 1, kflev+1
     1171      DO i = 1, kdlon
     1172        ZEMTD_i(i,k)=0.
     1173        ZEMTU_i(i,k)=0.
     1174        ZTRSO_i(i,k)=0.
     1175        ZTH_i(i,k)=0.
     1176        ZLWFT_i(i,k)=0.
     1177        ZSWFT_i(i,k)=0.
     1178        ZFLUX_i(i,1,k)=0.
     1179        ZFLUX_i(i,2,k)=0.
     1180        ZFLUC_i(i,1,k)=0.
     1181        ZFLUC_i(i,2,k)=0.
     1182        ZFSDWN_i(i,k)=0.
     1183        ZFCDWN_i(i,k)=0.
     1184        ZFCCDWN_i(i,k)=0.
     1185        ZFSUP_i(i,k)=0.
     1186        ZFCUP_i(i,k)=0.
     1187        ZFCCUP_i(i,k)=0.
     1188        ZFLCCDWN_i(i,k)=0.
     1189        ZFLCCUP_i(i,k)=0.
     1190      ENDDO
     1191      ENDDO
     1192!
     1193! Aerosols A REVOIR
     1194!      DO i = 1, kdlon
     1195!      DO k = 1, kflev
     1196!      DO kk=1, NSW
     1197!
     1198!      PTAU_TOT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,2,kk)
     1199!      PPIZA_TOT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,2,kk)
     1200!      PCGA_TOT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,2,kk)
     1201!
     1202!      PTAU_NAT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,1,kk)
     1203!      PPIZA_NAT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,1,kk)
     1204!      PCGA_NAT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,1,kk)
     1205!
     1206!      ENDDO
     1207!      ENDDO
     1208!      ENDDO
     1209!-end OB
     1210!
     1211!      DO i = 1, kdlon
     1212!      DO k = 1, kflev
     1213!      DO kk=1, NLW
     1214!
     1215!      PTAU_LW_TOT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,2,kk)
     1216!      PTAU_LW_NAT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,1,kk)
     1217!
     1218!      ENDDO
     1219!      ENDDO
     1220!      ENDDO
     1221!-end C. Kleinschmitt
     1222!     
     1223      DO i = 1, kdlon
     1224      ZCTRSO(i,1)=0.
     1225      ZCTRSO(i,2)=0.
     1226      ZCEMTR(i,1)=0.
     1227      ZCEMTR(i,2)=0.
     1228      ZTRSOD(i)=0.
     1229      ZLWFC(i,1)=0.
     1230      ZLWFC(i,2)=0.
     1231      ZSWFC(i,1)=0.
     1232      ZSWFC(i,2)=0.
     1233      PFSDNN(i)=0.
     1234      PFSDNV(i)=0.
     1235      DO kk = 1, NSW
     1236        PSFSWDIR(i,kk)=0.
     1237        PSFSWDIF(i,kk)=0.
     1238      ENDDO
     1239      ENDDO
     1240!----- Fin des mises a zero des tableaux output -------------------             
     1241
     1242! On met les donnees dans l'ordre des niveaux ecrad
     1243         print*,'On inverse sur la verticale '
     1244         paprs_i(:,1)=paprs(:,klev+1)
     1245         DO k=1,klev
     1246            paprs_i(1:klon,k+1) =paprs(1:klon,klev+1-k)
     1247            pplay_i(1:klon,k)   =pplay(1:klon,klev+1-k)
     1248            cldfra_i(1:klon,k)  =cldfra(1:klon,klev+1-k)
     1249            PDP_i(1:klon,k)     =PDP(1:klon,klev+1-k)
     1250            t_i(1:klon,k)       =t(1:klon,klev+1-k)
     1251            q_i(1:klon,k)       =q(1:klon,klev+1-k)
     1252            qsat_i(1:klon,k)    =qsat(1:klon,klev+1-k)
     1253            flwc_i(1:klon,k)    =flwc(1:klon,klev+1-k)
     1254            fiwc_i(1:klon,k)    =fiwc(1:klon,klev+1-k)
     1255            ref_liq_i(1:klon,k) =ref_liq(1:klon,klev+1-k)
     1256            ref_ice_i(1:klon,k) =ref_ice(1:klon,klev+1-k)
     1257!-OB
     1258            ref_liq_pi_i(1:klon,k) =ref_liq_pi(1:klon,klev+1-k)
     1259            ref_ice_pi_i(1:klon,k) =ref_ice_pi(1:klon,klev+1-k)
     1260         ENDDO
     1261         DO k=1,kflev
     1262            ZO3_DP_i(1:klon,k)=ZO3_DP(1:klon,kflev+1-k)
     1263            DO i=1,6
     1264            PAER_i(1:klon,k,i)=PAER(1:klon,kflev+1-k,i)
     1265            ENDDO
     1266         ENDDO
     1267! AI 02.2021
     1268! Calcul of ZTH_i
     1269         DO K=2,KLEV
     1270            ZTH_i(:,K)=&
     1271              & (t_i(:,K-1)*pplay_i(:,K-1)*(pplay_i(:,K)-paprs_i(:,K))&
     1272              & +t_i(:,K)*pplay_i(:,K)*(paprs_i(:,K)-pplay_i(:,K-1)))&
     1273              & *(1.0/(paprs_i(:,K)*(pplay_i(:,K)-pplay_i(:,K-1))))
     1274         ENDDO
     1275            ZTH_i(:,KLEV+1)=tsol(:)
     1276            ZTH_i(:,1)=t_i(:,1)-pplay_i(:,1)*(t_i(:,1)-ZTH_i(:,2))&
     1277                      & /(pplay_i(:,1)-paprs_i(:,2))
     1278
     1279      print *,'RADLWSW: avant RADIATION_SCHEME '
     1280      IF (lldebug) THEN
     1281        CALL writefield_phy('rmu0',rmu0,1)
     1282        CALL writefield_phy('tsol',tsol,1)
     1283        CALL writefield_phy('emissiv_out',ZEMIS,1)
     1284        CALL writefield_phy('emissiv_in',ZEMISW,1)
     1285        CALL writefield_phy('pctsrf_ter',pctsrf(:,is_ter),1)
     1286        CALL writefield_phy('pctsrf_oce',pctsrf(:,is_oce),1)
     1287        CALL writefield_phy('ZGELAM',ZGELAM,1)
     1288        CALL writefield_phy('ZGEMU',ZGEMU,1)
     1289        CALL writefield_phy('zmasq',zmasq,1)
     1290        CALL writefield_phy('paprs_i',paprs_i,klev+1)
     1291        CALL writefield_phy('pplay_i',pplay_i,klev)
     1292        CALL writefield_phy('t_i',t_i,klev)
     1293        CALL writefield_phy('ZTH_i',ZTH_i,klev+1)
     1294        CALL writefield_phy('cldfra_i',cldfra_i,klev)
     1295        CALL writefield_phy('paer_i',PAER_i,klev)
     1296        CALL writefield_phy('q_i',q_i,klev)
     1297        CALL writefield_phy('fiwc_i',fiwc_i,klev)
     1298        CALL writefield_phy('flwc_i',flwc_i,klev)
     1299        CALL writefield_phy('palbd_new',PALBD_NEW,NSW)
     1300        CALL writefield_phy('palbp_new',PALBP_NEW,NSW)
     1301        CALL writefield_phy('ZO3_DP',ZO3_DP,klev)
     1302      ENDIF
     1303
     1304      CALL RADIATION_SCHEME &
     1305      & (ist, iend, klon, klev, naero_tot, NSW, &
     1306! ??? naero_tot
     1307      & day_cur, current_time, &
     1308      & solaire, &
     1309      & rmu0, tsol, PALBD_NEW,PALBP_NEW, &   
     1310!       PEMIS_WINDOW (???), &
     1311      &  ZEMIS, ZEMISW, &
     1312!       PCCN_LAND, PCCN_SEA, & ???
     1313      & pctsrf(:,is_ter), pctsrf(:,is_oce), &
     1314!       longitude(rad), sin(latitude), PMASQ_ ???
     1315      & ZGELAM, ZGEMU, zmasq, &
     1316      & pplay_i, t_i, &
     1317!       PTEMPERATURE_H ?, &
     1318      & paprs_i, ZTH_i, q_i, qsat_i, &
     1319      & ZCO2, ZCH4, ZN2O, ZNO2, ZCFC11, ZCFC12, ZHCFC22, ZCCL4, ZO3_DP_i, &
     1320      & cldfra_i, flwc_i, fiwc_i, ZQ_RAIN, ZQ_SNOW, & 
     1321      & ref_liq_i, ref_ice_i, &
     1322!       aerosols
     1323      & ZAEROSOL_OLD, ZAEROSOL, &
     1324! Outputs
     1325!       Net flux
     1326      & ZSWFT_i, ZLWFT_i, ZSWFT0_i, ZLWFT0_i, &
     1327!       DN flux
     1328      & ZFSDWN_i, ZFLUX_i(:,2,:), ZFCDWN_i, ZFLUC_i(:,2,:), &
     1329!       UP flux
     1330      & ZFSUP_i, ZFLUX_i(:,1,:), ZFCUP_i, ZFLUC_i(:,1,:), &
     1331!       Surf Direct flux
     1332      & ZFLUX_DIR, ZFLUX_DIR_CLEAR, ZFLUX_DIR_INTO_SUN, &
     1333!       UV and para flux
     1334      & ZFLUX_UV, ZFLUX_PAR, ZFLUX_PAR_CLEAR, &
     1335!      & ZFLUX_SW_DN_TOA,
     1336      & ZEMIS_OUT, ZLWDERIVATIVE, &
     1337      & PSFSWDIF, PSFSWDIR)
     1338
     1339      print *,'========= RADLWSW: apres RADIATION_SCHEME ==================== '
     1340    if (0.eq.0) then
     1341      print *,' Net Flux '
     1342      print *,'ZSWFT_i =', ZSWFT_i
     1343      print *,'ZLWFT_i =', ZLWFT_i
     1344      print *,'ZSWFT0_i =', ZSWFT0_i
     1345      print*,'ZLWFT0_i =', ZLWFT0_i
     1346
     1347      print*,'DN Flux '
     1348      print*,'ZFSDWN_i =', ZFSDWN_i
     1349      print*,'ZFLUX_i(:,2,:)', ZFLUX_i(:,2,:)
     1350      print*,'ZFCDWN_i =', ZFCDWN_i
     1351      print*,'ZFLUC_i(:,2,:) =', ZFLUC_i(:,2,:)
     1352
     1353      print*,'UP Flux '
     1354      print*,'ZFSUP_i =', ZFSUP_i
     1355      print*,'ZFLUX_i(:,1,:) =', ZFLUX_i(:,1,:)
     1356      print*,'ZFCUP_i =', ZFCUP_i
     1357      print*,'ZFLUC_i(:,1,:) =', ZFLUC_i(:,1,:)
     1358
     1359      print*,'UV and para flux '
     1360      print*,'ZFLUX_DIR =', ZFLUX_DIR
     1361      print*,'ZFLUX_DIR_CLEAR', ZFLUX_DIR_CLEAR
     1362      print*,'ZFLUX_DIR_INTO_SUN', ZFLUX_DIR_INTO_SUN
     1363    endif
     1364
     1365      IF (lldebug) THEN
     1366        CALL writefield_phy('zlwft_i',ZLWFT_i,klev+1)
     1367        CALL writefield_phy('zlwft0_i',ZLWFT0_i,klev+1)
     1368        CALL writefield_phy('zswft_i',ZSWFT_i,klev+1)
     1369        CALL writefield_phy('zswft0_i',ZSWFT0_i,klev+1)
     1370        CALL writefield_phy('psfswdir',PSFSWDIR,6)
     1371        CALL writefield_phy('psfswdif',PSFSWDIF,6)
     1372        CALL writefield_phy('zflux_i',ZFLUX_i,klev+1)
     1373        CALL writefield_phy('zfluc_i',ZFLUC_i,klev+1)
     1374        CALL writefield_phy('zfsdwn_i',ZFSDWN_i,klev+1)
     1375        CALL writefield_phy('zfsup_i',ZFSUP_i,klev+1)
     1376        CALL writefield_phy('zfcdwn_i',ZFCDWN_i,klev+1)
     1377        CALL writefield_phy('zfcup_i',ZFCUP_i,klev+1)
     1378      ENDIF
     1379! ---------
     1380! On retablit l'ordre des niveaux lmd pour les tableaux de sortie
     1381! D autre part, on multiplie les resultats SW par fract pour etre coherent
     1382! avec l ancien rayonnement AR4. Si nuit, fract=0 donc pas de
     1383! rayonnement SW. (MPL 260609)
     1384      print*,'On retablit ordre des niveaux verticaux'
     1385      print*,'On multiplie les flux SW par fract ?'
     1386      DO k=0,klev
     1387         DO i=1,klon
     1388         ZEMTD(i,k+1)  = ZEMTD_i(i,k+1)
     1389         ZEMTU(i,k+1)  = ZEMTU_i(i,k+1)
     1390         ZTRSO(i,k+1)  = ZTRSO_i(i,k+1)
     1391         ZTH(i,k+1)    = ZTH_i(i,k+1)
     1392         ZFLUP(i,k+1)  = ZFLUX_i(i,1,k+1)
     1393         ZFLDN(i,k+1)  = ZFLUX_i(i,2,k+1)
     1394         ZFLUP0(i,k+1) = ZFLUC_i(i,1,k+1)
     1395         ZFLDN0(i,k+1) = ZFLUC_i(i,2,k+1)
     1396         ZFSDN(i,k+1)  = ZFSDWN_i(i,k+1)*fract(i)
     1397         ZFSDN0(i,k+1) = ZFCDWN_i(i,k+1)*fract(i)
     1398         ZFSDNC0(i,k+1)= ZFCCDWN_i(i,k+1)*fract(i)
     1399         ZFSUP (i,k+1) = ZFSUP_i(i,k+1)*fract(i)
     1400         ZFSUP0(i,k+1) = ZFCUP_i(i,k+1)*fract(i)
     1401         ZFSUPC0(i,k+1)= ZFCCUP_i(i,k+1)*fract(i)
     1402         ZFLDNC0(i,k+1)= ZFLCCDWN_i(i,k+1)
     1403         ZFLUPC0(i,k+1)= ZFLCCUP_i(i,k+1)
     1404         IF (ok_volcan) THEN
     1405            ZSWADAERO(i,k+1)=ZSWADAERO(i,k+1)*fract(i) !--NL
     1406         ENDIF
     1407         
     1408!   Nouveau calcul car visiblement ZSWFT et ZSWFC sont nuls dans RRTM cy32
     1409!   en sortie de radlsw.F90 - MPL 7.01.09
     1410         ZSWFT(i,k+1)  = (ZFSDWN_i(i,k+1)-ZFSUP_i(i,k+1))*fract(i)
     1411         ZSWFT0_i(i,k+1) = (ZFCDWN_i(i,k+1)-ZFCUP_i(i,k+1))*fract(i)
     1412         ZLWFT(i,k+1) =-ZFLUX_i(i,2,k+1)-ZFLUX_i(i,1,k+1)
     1413         ZLWFT0_i(i,k+1)=-ZFLUC_i(i,2,k+1)-ZFLUC_i(i,1,k+1)
     1414         ENDDO
     1415      ENDDO
     1416
     1417!--ajout OB
     1418      ZTOPSWADAERO(:) =ZTOPSWADAERO(:) *fract(:)
     1419      ZSOLSWADAERO(:) =ZSOLSWADAERO(:) *fract(:)
     1420      ZTOPSWAD0AERO(:)=ZTOPSWAD0AERO(:)*fract(:)
     1421      ZSOLSWAD0AERO(:)=ZSOLSWAD0AERO(:)*fract(:)
     1422      ZTOPSWAIAERO(:) =ZTOPSWAIAERO(:) *fract(:)
     1423      ZSOLSWAIAERO(:) =ZSOLSWAIAERO(:) *fract(:)
     1424      ZTOPSWCF_AERO(:,1)=ZTOPSWCF_AERO(:,1)*fract(:)
     1425      ZTOPSWCF_AERO(:,2)=ZTOPSWCF_AERO(:,2)*fract(:)
     1426      ZTOPSWCF_AERO(:,3)=ZTOPSWCF_AERO(:,3)*fract(:)
     1427      ZSOLSWCF_AERO(:,1)=ZSOLSWCF_AERO(:,1)*fract(:)
     1428      ZSOLSWCF_AERO(:,2)=ZSOLSWCF_AERO(:,2)*fract(:)
     1429      ZSOLSWCF_AERO(:,3)=ZSOLSWCF_AERO(:,3)*fract(:)
     1430
     1431! ---------
     1432! On renseigne les champs LMDz, pour avoir la meme chose qu'en sortie de
     1433! LW_LMDAR4 et SW_LMDAR4
     1434
     1435      !--fraction of diffuse radiation in surface SW downward radiation
     1436      DO i = 1, kdlon
     1437       IF (fract(i).GT.0.0) THEN
     1438         zdir=SUM(PSFSWDIR(i,:))
     1439         zdif=SUM(PSFSWDIF(i,:))
     1440         zsolswfdiff(i) = zdif/(zdir+zdif)
     1441       ELSE  !--night
     1442         zsolswfdiff(i) = 1.0
     1443       ENDIF
     1444      ENDDO
     1445!
     1446      DO i = 1, kdlon
     1447         zsolsw(i)    = ZSWFT(i,1)
     1448         zsolsw0(i)   = ZSWFT0_i(i,1)
     1449         ztopsw(i)    = ZSWFT(i,klev+1)
     1450         ztopsw0(i)   = ZSWFT0_i(i,klev+1)
     1451         zsollw(i)    = ZLWFT(i,1)
     1452         zsollw0(i)   = ZLWFT0_i(i,1)
     1453         ztoplw(i)    = ZLWFT(i,klev+1)*(-1)
     1454         ztoplw0(i)   = ZLWFT0_i(i,klev+1)*(-1)
     1455!         
     1456         zsollwdown(i)= -1.*ZFLDN(i,1)
     1457      ENDDO
     1458
     1459      DO k=1,kflev
     1460         DO i=1,kdlon
     1461           zheat(i,k)=(ZSWFT(i,k+1)-ZSWFT(i,k))*RDAY*RG/RCPD/PDP(i,k)
     1462           zheat0(i,k)=(ZSWFT0_i(i,k+1)-ZSWFT0_i(i,k))*RDAY*RG/RCPD/PDP(i,k)
     1463           zcool(i,k)=(ZLWFT(i,k)-ZLWFT(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
     1464           zcool0(i,k)=(ZLWFT0_i(i,k)-ZLWFT0_i(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
     1465           IF (ok_volcan) THEN
     1466              zheat_volc(i,k)=(ZSWADAERO(i,k+1)-ZSWADAERO(i,k))*RG/RCPD/PDP(i,k) !NL
     1467              zcool_volc(i,k)=(ZLWADAERO(i,k)-ZLWADAERO(i,k+1))*RG/RCPD/PDP(i,k) !NL
     1468           ENDIF
     1469         ENDDO
     1470      ENDDO
     1471#endif 
     1472  print*,'Fin traitement ECRAD'
     1473! Fin ECRAD
     1474  ENDIF        ! iflag_rrtm
     1475! ecrad
    10641476!======================================================================
    10651477
    … …  
    11021514          solswad_aero(iof+i) = zsolswadaero(i)
    11031515          solswad0_aero(iof+i) = zsolswad0aero(i)
    1104 ! MS the following lines seem to be wrong, why is iof on right hand side???
    1105 !          topsw_aero(iof+i,:) = ztopsw_aero(iof+i,:)
    1106 !          topsw0_aero(iof+i,:) = ztopsw0_aero(iof+i,:)
    1107 !          solsw_aero(iof+i,:) = zsolsw_aero(iof+i,:)
    1108 !          solsw0_aero(iof+i,:) = zsolsw0_aero(iof+i,:)
    11091516          topsw_aero(iof+i,:) = ztopsw_aero(i,:)
    11101517          topsw0_aero(iof+i,:) = ztopsw0_aero(i,:)
Note: See TracChangeset for help on using the changeset viewer.