Changeset 3908 for LMDZ6/trunk/libf


Ignore:
Timestamp:
May 20, 2021, 9:11:13 AM (3 years ago)
Author:
idelkadi
Message:

Online implementation of the radiative transfer code ECRAD in the LMDZ model.

  • Inclusion of the ecrad directory containing the sources of the ECRAD code
    • interface routine : radiation_scheme.F90
  • Adaptation of compilation scripts :
    • compilation under CPP key CPP_ECRAD
    • compilation with option "-rad ecard" or "-ecard true"
    • The "-rad old/rtm/ecran" build option will need to replace the "-rrtm true" and "-ecrad true" options in the future.
  • Runing LMDZ simulations with ecrad, you need :
    • logical key iflag_rrtm = 2 in physiq.def
    • namelist_ecrad (DefLists?)
    • the directory "data" containing the configuration files is temporarily placed in ../libfphylmd/ecrad/
  • Compilation and execution are tested in the 1D case. The repository under svn would allow to continue the implementation work: tests, verification of the results, ...
Location:
LMDZ6/trunk/libf/phylmd
Files:
392 added
3 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/conf_phys_m.F90

    r3901 r3908  
    26142614      ENDIF
    26152615    ELSE IF (iflag_rrtm .EQ. 1) THEN
     2616      IF (NSW.NE.2.AND.NSW.NE.4.AND.NSW.NE.6) THEN
     2617        WRITE(lunout,*) ' ERROR iflag_rrtm=1 and NSW<>2,4,6 not possible'
     2618        CALL abort_physic('conf_phys','choice NSW not valid',1)
     2619      ENDIF
     2620   ELSE IF (iflag_rrtm .EQ. 2) THEN
    26162621      IF (NSW.NE.2.AND.NSW.NE.4.AND.NSW.NE.6) THEN
    26172622        WRITE(lunout,*) ' ERROR iflag_rrtm=1 and NSW<>2,4,6 not possible'

  • LMDZ6/trunk/libf/phylmd/physiq_mod.F90

    r3900 r3908  
    39673967
    39683968       IF (ok_newmicro) then
    3969           IF (iflag_rrtm.NE.0) THEN
     3969! AI          IF (iflag_rrtm.NE.0) THEN
     3970          IF (iflag_rrtm.EQ.1) THEN
    39703971#ifdef CPP_RRTM
    39713972             IF (ok_cdnc.AND.NRADLP.NE.3) THEN

  • LMDZ6/trunk/libf/phylmd/radlwsw_m.F90

    r3756 r3908  
    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 IF (lldebug) THEN
     452  print*,'Entree dans radlwsw '
     453  print*,'************* INITIALISATIONS *****************************'
     454  print*,'klon, kdlon, klev, kflev =',klon, kdlon, klev, kflev
     455 ENDIF
     456
    403457  CALL assert(size(wo, 1) == klon, size(wo, 2) == klev, "radlwsw wo")
    404   ! initialisation
     458 
    405459  ist=1
    406460  iend=klon
    407461  ktdia=1
    408462  kmode=ist
     463! Aeros
    409464  tauaero(:,:,:,:)=0.
    410465  pizaero(:,:,:,:)=0.
    … …  
    462517  ENDIF
    463518
     519 IF (lldebug) THEN
     520  print*,'************** Debut boucle de 1 a ', nb_gr
     521 ENDIF
     522
    464523  DO j = 1, nb_gr
    465524    iof = kdlon*(j-1)
    466525    DO i = 1, kdlon
    467526      zfract(i) = fract(iof+i)
    468 !     zfract(i) = 1.     !!!!!!  essai MPL 19052010
    469527      zrmu0(i) = rmu0(iof+i)
    470528
    471529
    472 !albedo SB >>>
    473 !
    474530      IF (iflag_rrtm==0) THEN
    475 !
     531!     Albedo
    476532        PALBD(i,1)=alb_dif(iof+i,1)
    477533        PALBD(i,2)=alb_dif(iof+i,2)
    478534        PALBP(i,1)=alb_dir(iof+i,1)
    479535        PALBP(i,2)=alb_dir(iof+i,2)
    480 !
    481       ELSEIF (iflag_rrtm==1) THEn
    482 !
     536! AI 02.2021 cas iflag_rrtm=1 et 2
     537       ELSEIF (iflag_rrtm==1.OR.iflag_rrtm==2) THEN
    483538        DO kk=1,NSW
    484539          PALBD_NEW(i,kk)=alb_dif(iof+i,kk)
    … …  
    488543      ENDIF
    489544!albedo SB <<<
    490 
    491545
    492546      PEMIS(i) = 1.0    !!!!! A REVOIR (MPL)
    … …  
    513567        POZON(i,k, :) = wo(iof+i, k, :) * RG * dobson_u * 1e3 &
    514568             / (paprs(iof+i, k) - paprs(iof+i, k+1))
     569        ZO3_DP(i,k) = wo(iof+i, k, 1) * RG * dobson_u * 1e3
    515570!       A activer pour CCMVAL on prend l'ozone impose (MPL 07042010)
    516571!       POZON(i,k,:) = wo(i,k,:) 
    … …  
    569624      ENDDO
    570625    ENDDO
     626!
     627! AI 02.2021
     628#ifdef CPP_ECRAD
     629  ZEMIS = 1.0
     630  ZEMISW = 1.0
     631  ZGELAM = longitude
     632  ZGEMU = sin(latitude)
     633  ZCO2 = RCO2
     634  ZCH4 = RCH4
     635  ZN2O = RN2O
     636  ZNO2 = 0.0
     637  ZCFC11 = RCFC11
     638  ZCFC12 = RCFC12
     639  ZHCFC22 = 0.0
     640  ZCCL4 = 0.0
     641  ZQ_RAIN = 0.0
     642  ZQ_SNOW = 0.0
     643  ZAEROSOL_OLD = 0.0
     644  ZAEROSOL = 0.0
     645#endif
    571646!
    572647!===== iflag_rrtm ================================================
    … …  
    693768       ENDDO 
    694769!
    695     ELSE
     770    ELSE IF (iflag_rrtm == 1) then
    696771#ifdef CPP_RRTM
    697772!      if (prt_level.gt.10)write(lunout,*)'CPP_RRTM=.T.'
    … …  
    804879            ENDDO
    805880         ENDDO
     881
    806882!       print *,'RADLWSW: avant RECMWFL, RI0,rmu0=',solaire,rmu0
    807883
    … …  
    850926! Nouvel appel a RECMWF (celui du cy32t0)
    851927         CALL RECMWF_AERO (ist , iend, klon , ktdia  , klev   , kmode ,&
     928!                          KST,  KEND, KPROMA, KTDIA , KLEV,  KMODE
    852929         PALBD_NEW,PALBP_NEW, paprs_i , pplay_i , RCO2   , cldfra_i,&
     930!        PALBD ,   PALBP ,    PAPRS ,   PAPRSF ,  PCCO2 , PCLFR
    853931         POZON_i  , PAER_i  , PDP_i   , PEMIS   , rmu0   ,&
     932!        PQO3  ,    PAER  ,   PDP   ,   PEMIS  , PMU0
    854933          q_i     , qsat_i  , fiwc_i  , flwc_i  , zmasq  , t_i  ,tsol,&
     934!         PQ       , PQS   , PQIWP ,     PQLWP , PSLM   , PT    , PTS,
    855935         ref_liq_i, ref_ice_i, &
     936!        PREF_LIQ, PREF_ICE
    856937         ref_liq_pi_i, ref_ice_pi_i, &   ! rajoute par OB pour diagnostiquer effet indirect
     938!        PREF_LIQ_PI, PREF_ICE_PI
    857939         ZEMTD_i  , ZEMTU_i , ZTRSO_i ,&
     940!        PEMTD , PEMTU , PTRSO
    858941         ZTH_i    , ZCTRSO  , ZCEMTR  , ZTRSOD  ,&
     942!        PTH   ,     PCTRSO,   PCEMTR, PTRSOD
    859943         ZLWFC    , ZLWFT_i , ZSWFC   , ZSWFT_i ,&
     944!        PLWFC,     PLWFT,    PSWFC,    PSWFT,
    860945         PSFSWDIR , PSFSWDIF, PFSDNN  , PFSDNV  ,&
     946!        PSFSWDIR,  PSFSWDIF, PFSDNN,    PFSDNV
    861947         PPIZA_TOT, PCGA_TOT,PTAU_TOT,&
     948!        PPIZA_TOT, PCGA_TOT,PTAU_TOT
    862949         PPIZA_NAT, PCGA_NAT,PTAU_NAT,           &  ! rajoute par OB pour diagnostiquer effet direct
     950!        PPIZA_NAT,PCGA_NAT,PTAU_NAT
    863951         PTAU_LW_TOT, PTAU_LW_NAT,               &  ! rajoute par C. Kleinschmitt
     952!        PTAU_LW_TOT, PTAU_LW_NAT,
    864953         ZFLUX_i  , ZFLUC_i ,&
     954!        PFLUX,      PFLUC
    865955         ZFSDWN_i , ZFSUP_i , ZFCDWN_i, ZFCUP_i, ZFCCDWN_i, ZFCCUP_i, ZFLCCDWN_i, ZFLCCUP_i, &
     956!        PFSDN ,    PFSUP ,   PFSCDN , PFSCUP,   PFSCCDN,   PFSCCUP,  PFLCCDN,    PFLCCUP
    866957         ZTOPSWADAERO,ZSOLSWADAERO,&  ! rajoute par OB pour diagnostics
     958!        PTOPSWADAERO,PSOLSWADAERO
    867959         ZTOPSWAD0AERO,ZSOLSWAD0AERO,&
    868960         ZTOPSWAIAERO,ZSOLSWAIAERO, &
    … …  
    874966         ZLWADAERO, & !--NL
    875967         ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat, flag_aer_feedback) ! flags aerosols
     968
     969!--OB diagnostics
     970! & PTOPSWAIAERO,PSOLSWAIAERO,&
     971! & PTOPSWCFAERO,PSOLSWCFAERO,&
     972! & PSWADAERO,& !--NL
     973!!--LW diagnostics CK
     974! & PTOPLWADAERO,PSOLLWADAERO,&
     975! & PTOPLWAD0AERO,PSOLLWAD0AERO,&
     976! & PTOPLWAIAERO,PSOLLWAIAERO,&
     977! & PLWADAERO,& !--NL
     978!!..end
     979! & ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat,&
     980! & flag_aer_feedback)
     981
    876982           
    877983!        print *,'RADLWSW: apres RECMWF'
    … …  
    10611167    call abort_physic(modname, abort_message, 1)
    10621168#endif
    1063     ENDIF ! iflag_rrtm
     1169!======================================================================
     1170! AI fev 2021
     1171    ELSE IF(iflag_rrtm == 2) THEN
     1172    print*,'Traitement cas iflag_rrtm = ',iflag_rrtm
     1173!    print*,'Mise a zero des flux '
     1174#ifdef CPP_ECRAD
     1175      DO k = 1, kflev+1
     1176      DO i = 1, kdlon
     1177        ZEMTD_i(i,k)=0.
     1178        ZEMTU_i(i,k)=0.
     1179        ZTRSO_i(i,k)=0.
     1180        ZTH_i(i,k)=0.
     1181        ZLWFT_i(i,k)=0.
     1182        ZSWFT_i(i,k)=0.
     1183        ZFLUX_i(i,1,k)=0.
     1184        ZFLUX_i(i,2,k)=0.
     1185        ZFLUC_i(i,1,k)=0.
     1186        ZFLUC_i(i,2,k)=0.
     1187        ZFSDWN_i(i,k)=0.
     1188        ZFCDWN_i(i,k)=0.
     1189        ZFCCDWN_i(i,k)=0.
     1190        ZFSUP_i(i,k)=0.
     1191        ZFCUP_i(i,k)=0.
     1192        ZFCCUP_i(i,k)=0.
     1193        ZFLCCDWN_i(i,k)=0.
     1194        ZFLCCUP_i(i,k)=0.
     1195      ENDDO
     1196      ENDDO
     1197!
     1198! Aerosols A REVOIR
     1199!      DO i = 1, kdlon
     1200!      DO k = 1, kflev
     1201!      DO kk=1, NSW
     1202!
     1203!      PTAU_TOT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,2,kk)
     1204!      PPIZA_TOT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,2,kk)
     1205!      PCGA_TOT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,2,kk)
     1206!
     1207!      PTAU_NAT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,1,kk)
     1208!      PPIZA_NAT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,1,kk)
     1209!      PCGA_NAT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,1,kk)
     1210!
     1211!      ENDDO
     1212!      ENDDO
     1213!      ENDDO
     1214!-end OB
     1215!
     1216!      DO i = 1, kdlon
     1217!      DO k = 1, kflev
     1218!      DO kk=1, NLW
     1219!
     1220!      PTAU_LW_TOT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,2,kk)
     1221!      PTAU_LW_NAT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,1,kk)
     1222!
     1223!      ENDDO
     1224!      ENDDO
     1225!      ENDDO
     1226!-end C. Kleinschmitt
     1227!     
     1228      DO i = 1, kdlon
     1229      ZCTRSO(i,1)=0.
     1230      ZCTRSO(i,2)=0.
     1231      ZCEMTR(i,1)=0.
     1232      ZCEMTR(i,2)=0.
     1233      ZTRSOD(i)=0.
     1234      ZLWFC(i,1)=0.
     1235      ZLWFC(i,2)=0.
     1236      ZSWFC(i,1)=0.
     1237      ZSWFC(i,2)=0.
     1238      PFSDNN(i)=0.
     1239      PFSDNV(i)=0.
     1240      DO kk = 1, NSW
     1241        PSFSWDIR(i,kk)=0.
     1242        PSFSWDIF(i,kk)=0.
     1243      ENDDO
     1244      ENDDO
     1245!----- Fin des mises a zero des tableaux output -------------------             
     1246
     1247! On met les donnees dans l'ordre des niveaux ecrad
     1248!         print*,'On inverse sur la verticale '
     1249         paprs_i(:,1)=paprs(:,klev+1)
     1250         DO k=1,klev
     1251            paprs_i(1:klon,k+1) =paprs(1:klon,klev+1-k)
     1252            pplay_i(1:klon,k)   =pplay(1:klon,klev+1-k)
     1253            cldfra_i(1:klon,k)  =cldfra(1:klon,klev+1-k)
     1254            PDP_i(1:klon,k)     =PDP(1:klon,klev+1-k)
     1255            t_i(1:klon,k)       =t(1:klon,klev+1-k)
     1256            q_i(1:klon,k)       =q(1:klon,klev+1-k)
     1257            qsat_i(1:klon,k)    =qsat(1:klon,klev+1-k)
     1258            flwc_i(1:klon,k)    =flwc(1:klon,klev+1-k)
     1259            fiwc_i(1:klon,k)    =fiwc(1:klon,klev+1-k)
     1260            ref_liq_i(1:klon,k) =ref_liq(1:klon,klev+1-k)
     1261            ref_ice_i(1:klon,k) =ref_ice(1:klon,klev+1-k)
     1262!-OB
     1263            ref_liq_pi_i(1:klon,k) =ref_liq_pi(1:klon,klev+1-k)
     1264            ref_ice_pi_i(1:klon,k) =ref_ice_pi(1:klon,klev+1-k)
     1265         ENDDO
     1266         DO k=1,kflev
     1267            ZO3_DP_i(1:klon,k)=ZO3_DP(1:klon,kflev+1-k)
     1268            DO i=1,6
     1269            PAER_i(1:klon,k,i)=PAER(1:klon,kflev+1-k,i)
     1270            ENDDO
     1271         ENDDO
     1272! AI 02.2021
     1273! Calcul of ZTH_i
     1274         DO K=2,KLEV
     1275            ZTH_i(:,K)=&
     1276              & (t_i(:,K-1)*pplay_i(:,K-1)*(pplay_i(:,K)-paprs_i(:,K))&
     1277              & +t_i(:,K)*pplay_i(:,K)*(paprs_i(:,K)-pplay_i(:,K-1)))&
     1278              & *(1.0/(paprs_i(:,K)*(pplay_i(:,K)-pplay_i(:,K-1))))
     1279         ENDDO
     1280            ZTH_i(:,KLEV+1)=tsol(:)
     1281            ZTH_i(:,1)=t_i(:,1)-pplay_i(:,1)*(t_i(:,1)-ZTH_i(:,2))&
     1282                      & /(pplay_i(:,1)-paprs_i(:,2))
     1283
     1284      print *,'RADLWSW: avant RADIATION_SCHEME '
     1285      IF (lldebug) THEN
     1286        CALL writefield_phy('rmu0',rmu0,1)
     1287        CALL writefield_phy('tsol',tsol,1)
     1288        CALL writefield_phy('emissiv_out',ZEMIS,1)
     1289        CALL writefield_phy('emissiv_in',ZEMISW,1)
     1290        CALL writefield_phy('pctsrf_ter',pctsrf(:,is_ter),1)
     1291        CALL writefield_phy('pctsrf_oce',pctsrf(:,is_oce),1)
     1292        CALL writefield_phy('ZGELAM',ZGELAM,1)
     1293        CALL writefield_phy('ZGEMU',ZGEMU,1)
     1294        CALL writefield_phy('zmasq',zmasq,1)
     1295        CALL writefield_phy('paprs_i',paprs_i,klev+1)
     1296        CALL writefield_phy('pplay_i',pplay_i,klev)
     1297        CALL writefield_phy('t_i',t_i,klev)
     1298        CALL writefield_phy('ZTH_i',ZTH_i,klev+1)
     1299        CALL writefield_phy('cldfra_i',cldfra_i,klev)
     1300        CALL writefield_phy('paer_i',PAER_i,klev)
     1301        CALL writefield_phy('q_i',q_i,klev)
     1302        CALL writefield_phy('fiwc_i',fiwc_i,klev)
     1303        CALL writefield_phy('flwc_i',flwc_i,klev)
     1304        CALL writefield_phy('palbd_new',PALBD_NEW,NSW)
     1305        CALL writefield_phy('palbp_new',PALBP_NEW,NSW)
     1306        CALL writefield_phy('ZO3_DP',ZO3_DP,klev)
     1307      ENDIF
     1308
     1309      CALL RADIATION_SCHEME &
     1310      & (ist, iend, klon, klev, naero_tot, NSW, &
     1311! ??? naero_tot
     1312      & day_cur, current_time, &
     1313      & solaire, &
     1314      & rmu0, tsol, PALBD_NEW,PALBP_NEW, &   
     1315!       PEMIS_WINDOW (???), &
     1316      &  ZEMIS, ZEMISW, &
     1317!       PCCN_LAND, PCCN_SEA, & ???
     1318      & pctsrf(:,is_ter), pctsrf(:,is_oce), &
     1319!       longitude(rad), sin(latitude), PMASQ_ ???
     1320      & ZGELAM, ZGEMU, zmasq, &
     1321      & pplay_i, t_i, &
     1322!       PTEMPERATURE_H ?, &
     1323      & paprs_i, ZTH_i, q_i, qsat_i, &
     1324      & ZCO2, ZCH4, ZN2O, ZNO2, ZCFC11, ZCFC12, ZHCFC22, ZCCL4, ZO3_DP_i, &
     1325      & cldfra_i, flwc_i, fiwc_i, ZQ_RAIN, ZQ_SNOW, & 
     1326      & ref_liq_i, ref_ice_i, &
     1327!       aerosols
     1328      & ZAEROSOL_OLD, ZAEROSOL, &
     1329! Outputs
     1330!       Net flux
     1331      & ZSWFT_i, ZLWFT_i, ZSWFT0_i, ZLWFT0_i, &
     1332!       DN flux
     1333      & ZFSDWN_i, ZFLUX_i(:,2,:), ZFCDWN_i, ZFLUC_i(:,2,:), &
     1334!       UP flux
     1335      & ZFSUP_i, ZFLUX_i(:,1,:), ZFCUP_i, ZFLUC_i(:,1,:), &
     1336!       Surf Direct flux
     1337      & ZFLUX_DIR, ZFLUX_DIR_CLEAR, ZFLUX_DIR_INTO_SUN, &
     1338!       UV and para flux
     1339      & ZFLUX_UV, ZFLUX_PAR, ZFLUX_PAR_CLEAR, &
     1340!      & ZFLUX_SW_DN_TOA,
     1341      & ZEMIS_OUT, ZLWDERIVATIVE, &
     1342      & PSFSWDIF, PSFSWDIR)
     1343
     1344      print *,'========= RADLWSW: apres RADIATION_SCHEME ==================== '
     1345  IF (lldebug) THEN
     1346    if (0.eq.0) then
     1347      print *,' Net Flux '
     1348      print *,'ZSWFT_i =', ZSWFT_i
     1349      print *,'ZLWFT_i =', ZLWFT_i
     1350      print *,'ZSWFT0_i =', ZSWFT0_i
     1351      print*,'ZLWFT0_i =', ZLWFT0_i
     1352
     1353      print*,'DN Flux '
     1354      print*,'ZFSDWN_i =', ZFSDWN_i
     1355      print*,'ZFLUX_i(:,2,:)', ZFLUX_i(:,2,:)
     1356      print*,'ZFCDWN_i =', ZFCDWN_i
     1357      print*,'ZFLUC_i(:,2,:) =', ZFLUC_i(:,2,:)
     1358
     1359      print*,'UP Flux '
     1360      print*,'ZFSUP_i =', ZFSUP_i
     1361      print*,'ZFLUX_i(:,1,:) =', ZFLUX_i(:,1,:)
     1362      print*,'ZFCUP_i =', ZFCUP_i
     1363      print*,'ZFLUC_i(:,1,:) =', ZFLUC_i(:,1,:)
     1364
     1365      print*,'UV and para flux '
     1366      print*,'ZFLUX_DIR =', ZFLUX_DIR
     1367      print*,'ZFLUX_DIR_CLEAR', ZFLUX_DIR_CLEAR
     1368      print*,'ZFLUX_DIR_INTO_SUN', ZFLUX_DIR_INTO_SUN
     1369    endif
     1370   ENDIF
     1371
     1372      IF (lldebug) THEN
     1373        CALL writefield_phy('zlwft_i',ZLWFT_i,klev+1)
     1374        CALL writefield_phy('zlwft0_i',ZLWFT0_i,klev+1)
     1375        CALL writefield_phy('zswft_i',ZSWFT_i,klev+1)
     1376        CALL writefield_phy('zswft0_i',ZSWFT0_i,klev+1)
     1377        CALL writefield_phy('psfswdir',PSFSWDIR,6)
     1378        CALL writefield_phy('psfswdif',PSFSWDIF,6)
     1379        CALL writefield_phy('zflux_i',ZFLUX_i,klev+1)
     1380        CALL writefield_phy('zfluc_i',ZFLUC_i,klev+1)
     1381        CALL writefield_phy('zfsdwn_i',ZFSDWN_i,klev+1)
     1382        CALL writefield_phy('zfsup_i',ZFSUP_i,klev+1)
     1383        CALL writefield_phy('zfcdwn_i',ZFCDWN_i,klev+1)
     1384        CALL writefield_phy('zfcup_i',ZFCUP_i,klev+1)
     1385      ENDIF
     1386! ---------
     1387! On retablit l'ordre des niveaux lmd pour les tableaux de sortie
     1388! D autre part, on multiplie les resultats SW par fract pour etre coherent
     1389! avec l ancien rayonnement AR4. Si nuit, fract=0 donc pas de
     1390! rayonnement SW. (MPL 260609)
     1391      print*,'On retablit ordre des niveaux verticaux'
     1392      print*,'On multiplie les flux SW par fract ?'
     1393      DO k=0,klev
     1394         DO i=1,klon
     1395         ZEMTD(i,k+1)  = ZEMTD_i(i,k+1)
     1396         ZEMTU(i,k+1)  = ZEMTU_i(i,k+1)
     1397         ZTRSO(i,k+1)  = ZTRSO_i(i,k+1)
     1398         ZTH(i,k+1)    = ZTH_i(i,k+1)
     1399         ZFLUP(i,k+1)  = ZFLUX_i(i,1,k+1)
     1400         ZFLDN(i,k+1)  = ZFLUX_i(i,2,k+1)
     1401         ZFLUP0(i,k+1) = ZFLUC_i(i,1,k+1)
     1402         ZFLDN0(i,k+1) = ZFLUC_i(i,2,k+1)
     1403         ZFSDN(i,k+1)  = ZFSDWN_i(i,k+1)*fract(i)
     1404         ZFSDN0(i,k+1) = ZFCDWN_i(i,k+1)*fract(i)
     1405         ZFSDNC0(i,k+1)= ZFCCDWN_i(i,k+1)*fract(i)
     1406         ZFSUP (i,k+1) = ZFSUP_i(i,k+1)*fract(i)
     1407         ZFSUP0(i,k+1) = ZFCUP_i(i,k+1)*fract(i)
     1408         ZFSUPC0(i,k+1)= ZFCCUP_i(i,k+1)*fract(i)
     1409         ZFLDNC0(i,k+1)= ZFLCCDWN_i(i,k+1)
     1410         ZFLUPC0(i,k+1)= ZFLCCUP_i(i,k+1)
     1411         IF (ok_volcan) THEN
     1412            ZSWADAERO(i,k+1)=ZSWADAERO(i,k+1)*fract(i) !--NL
     1413         ENDIF
     1414         
     1415!   Nouveau calcul car visiblement ZSWFT et ZSWFC sont nuls dans RRTM cy32
     1416!   en sortie de radlsw.F90 - MPL 7.01.09
     1417         ZSWFT(i,k+1)  = (ZFSDWN_i(i,k+1)-ZFSUP_i(i,k+1))*fract(i)
     1418         ZSWFT0_i(i,k+1) = (ZFCDWN_i(i,k+1)-ZFCUP_i(i,k+1))*fract(i)
     1419         ZLWFT(i,k+1) =-ZFLUX_i(i,2,k+1)-ZFLUX_i(i,1,k+1)
     1420         ZLWFT0_i(i,k+1)=-ZFLUC_i(i,2,k+1)-ZFLUC_i(i,1,k+1)
     1421         ENDDO
     1422      ENDDO
     1423
     1424!--ajout OB
     1425      ZTOPSWADAERO(:) =ZTOPSWADAERO(:) *fract(:)
     1426      ZSOLSWADAERO(:) =ZSOLSWADAERO(:) *fract(:)
     1427      ZTOPSWAD0AERO(:)=ZTOPSWAD0AERO(:)*fract(:)
     1428      ZSOLSWAD0AERO(:)=ZSOLSWAD0AERO(:)*fract(:)
     1429      ZTOPSWAIAERO(:) =ZTOPSWAIAERO(:) *fract(:)
     1430      ZSOLSWAIAERO(:) =ZSOLSWAIAERO(:) *fract(:)
     1431      ZTOPSWCF_AERO(:,1)=ZTOPSWCF_AERO(:,1)*fract(:)
     1432      ZTOPSWCF_AERO(:,2)=ZTOPSWCF_AERO(:,2)*fract(:)
     1433      ZTOPSWCF_AERO(:,3)=ZTOPSWCF_AERO(:,3)*fract(:)
     1434      ZSOLSWCF_AERO(:,1)=ZSOLSWCF_AERO(:,1)*fract(:)
     1435      ZSOLSWCF_AERO(:,2)=ZSOLSWCF_AERO(:,2)*fract(:)
     1436      ZSOLSWCF_AERO(:,3)=ZSOLSWCF_AERO(:,3)*fract(:)
     1437
     1438! ---------
     1439! On renseigne les champs LMDz, pour avoir la meme chose qu'en sortie de
     1440! LW_LMDAR4 et SW_LMDAR4
     1441
     1442      !--fraction of diffuse radiation in surface SW downward radiation
     1443      DO i = 1, kdlon
     1444       IF (fract(i).GT.0.0) THEN
     1445         zdir=SUM(PSFSWDIR(i,:))
     1446         zdif=SUM(PSFSWDIF(i,:))
     1447         zsolswfdiff(i) = zdif/(zdir+zdif)
     1448       ELSE  !--night
     1449         zsolswfdiff(i) = 1.0
     1450       ENDIF
     1451      ENDDO
     1452!
     1453      DO i = 1, kdlon
     1454         zsolsw(i)    = ZSWFT(i,1)
     1455         zsolsw0(i)   = ZSWFT0_i(i,1)
     1456         ztopsw(i)    = ZSWFT(i,klev+1)
     1457         ztopsw0(i)   = ZSWFT0_i(i,klev+1)
     1458         zsollw(i)    = ZLWFT(i,1)
     1459         zsollw0(i)   = ZLWFT0_i(i,1)
     1460         ztoplw(i)    = ZLWFT(i,klev+1)*(-1)
     1461         ztoplw0(i)   = ZLWFT0_i(i,klev+1)*(-1)
     1462!         
     1463         zsollwdown(i)= -1.*ZFLDN(i,1)
     1464      ENDDO
     1465
     1466      DO k=1,kflev
     1467         DO i=1,kdlon
     1468           zheat(i,k)=(ZSWFT(i,k+1)-ZSWFT(i,k))*RDAY*RG/RCPD/PDP(i,k)
     1469           zheat0(i,k)=(ZSWFT0_i(i,k+1)-ZSWFT0_i(i,k))*RDAY*RG/RCPD/PDP(i,k)
     1470           zcool(i,k)=(ZLWFT(i,k)-ZLWFT(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
     1471           zcool0(i,k)=(ZLWFT0_i(i,k)-ZLWFT0_i(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
     1472           IF (ok_volcan) THEN
     1473              zheat_volc(i,k)=(ZSWADAERO(i,k+1)-ZSWADAERO(i,k))*RG/RCPD/PDP(i,k) !NL
     1474              zcool_volc(i,k)=(ZLWADAERO(i,k)-ZLWADAERO(i,k+1))*RG/RCPD/PDP(i,k) !NL
     1475           ENDIF
     1476         ENDDO
     1477      ENDDO
     1478#endif 
     1479  print*,'Fin traitement ECRAD'
     1480! Fin ECRAD
     1481  ENDIF        ! iflag_rrtm
     1482! ecrad
    10641483!======================================================================
    10651484
    … …  
    11021521          solswad_aero(iof+i) = zsolswadaero(i)
    11031522          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,:)
    11091523          topsw_aero(iof+i,:) = ztopsw_aero(i,:)
    11101524          topsw0_aero(iof+i,:) = ztopsw0_aero(i,:)
Note: See TracChangeset for help on using the changeset viewer.