Ignore:
Timestamp:
Apr 1, 2026, 2:42:39 PM (12 days ago)
Author:
evignon
Message:

corrections pour retrouver la convergence en NP et la compilation en debug AVEC aerosols

File:
1 edited

Legend:

Unmodified
Added
Removed
  • LMDZ6/trunk/libf/phylmd/physiq_mod.F90

    r6146 r6148  
    37563756    call surf_wind(klon,nsurfwind,zu10m,zv10m,wake_s,wake_Cstar,zustar,ale_bl,surf_wind_value,surf_wind_proba)
    37573757
     3758
     3759
     3760
    37583761    !===========================================================================
    37593762    ! Large scale condensation and precipitation
     
    37973800     dqsmelt, dqsfreez, zx_rh, zx_rhl, zx_rhi)
    37983801
    3799     !===============================================================================
    3800     ! from Clouds to Radiation (cloud2rad)
    3801     ! computation of cloud variables that are useful for the radiative transfer
    3802 
    3803 
    3804     DO k = 1, klev
    3805        DO i = 1, klon
    3806           cldfra(i,k) = rneb(i,k)
    3807           ! keep only liquid droplets in radocond if not liqice_in_radocond
    3808           IF (.NOT.liqice_in_radocond) radocond(i,k) = ql_seri(i,k)
    3809        ENDDO
    3810     ENDDO
    3811 
    3812 
    3813     ! Option to activate the radiative effect of blowing snow (ok_rad_bs)
    3814     ! makes sense only if the new large scale condensation scheme is active
    3815     ! with the ok_icefra_lscp flag active as well
    3816 
    3817     IF (ok_bs .AND. ok_rad_bs) THEN
    3818        !IF (ok_icefra_lscp) THEN
    3819           DO k=1,klev
    3820              DO i=1,klon
    3821                 radocond(i,k)=radocond(i,k)+qbs_seri(i,k)
    3822                 radicefrac(i,k)=(radocond(i,k)*radicefrac(i,k)+qbs_seri(i,k))/(radocond(i,k))
    3823                 qbsfra=min(qbs_seri(i,k)/qbst_bs,1.0)
    3824                 cldfra(i,k)=max(cldfra(i,k),qbsfra)
    3825              ENDDO
    3826           ENDDO
    3827        !ELSE
    3828        !   WRITE(lunout,*)"PAY ATTENTION, you try to activate the radiative effect of blowing snow"
    3829        !   WRITE(lunout,*)"with ok_icefra_lscp=false. You must use ok_icefra_lscp=y and ok_new_lscp=y."
    3830        !   abort_message='inconsistency in cloud phase for blowing snow'
    3831        !   CALL abort_physic(modname,abort_message,1)
    3832        !ENDIF
    3833 
    3834     ENDIF
    3835 
    3836     IF (mydebug) THEN
    3837        CALL writefield_phy('u_seri',u_seri,nbp_lev)
    3838        CALL writefield_phy('v_seri',v_seri,nbp_lev)
    3839        CALL writefield_phy('t_seri',t_seri,nbp_lev)
    3840        CALL writefield_phy('q_seri',q_seri,nbp_lev)
    3841     ENDIF
    3842 
    3843 
    3844     ! 1. NUAGES CONVECTIFS
    3845     !     IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke
    3846     IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke
    3847        snow_tiedtke=0.
    3848        IF (iflag_cld_th.eq.-1) THEN
    3849           rain_tiedtke=rain_con
    3850        ELSE
    3851           rain_tiedtke=0.
    3852           DO k=1,klev
    3853              DO i=1,klon
    3854                 IF (d_q_con(i,k).lt.0.) THEN
    3855                    rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
    3856                         *(paprs(i,k)-paprs(i,k+1))/rg
    3857                 ENDIF
    3858              ENDDO
    3859           ENDDO
    3860        ENDIF
    3861 
    3862        ! Nuages diagnostiques pour Tiedtke
    3863        CALL diagcld1(paprs,pplay, &
    3864             rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
    3865             diafra,dialiq)
    3866        DO k = 1, klev
    3867           DO i = 1, klon
    3868              IF (diafra(i,k).GT.cldfra(i,k)) THEN
    3869                 radocond(i,k) = dialiq(i,k)
    3870                 cldfra(i,k) = diafra(i,k)
    3871              ENDIF
    3872           ENDDO
    3873        ENDDO
    3874 
    3875     ELSE IF (iflag_cld_th.ge.3) THEN
    3876        !  On prend pour les nuages convectifs le max du calcul de la
    3877        !  convection et du calcul du pas de temps precedent diminue d'un facteur
    3878        !  facttemps
    3879        facteur = pdtphys *facttemps
    3880        DO k=1,klev
    3881           DO i=1,klon
    3882              rnebcon(i,k)=rnebcon(i,k)*facteur
    3883              IF (rnebcon0(i,k)*clwcon0(i,k).GT.rnebcon(i,k)*clwcon(i,k)) THEN
    3884                 rnebcon(i,k)=rnebcon0(i,k)
    3885                 clwcon(i,k)=clwcon0(i,k)
    3886              ENDIF
    3887           ENDDO
    3888        ENDDO
    3889 
    3890        !   On prend la somme des fractions nuageuses et des contenus en eau
    3891 
    3892        IF (iflag_cld_th>=5) THEN
    3893 
    3894           DO k=1,klev
    3895              ptconvth(:,k)=fm_therm(:,k+1)>0.
    3896           ENDDO
    3897 
    3898           IF (iflag_coupl==4) THEN
    3899 
    3900              ! Dans le cas iflag_coupl==4, on prend la somme des convertures
    3901              ! convectives et lsc dans la partie des thermiques
    3902              ! Le controle par iflag_coupl est peut etre provisoire.
    3903              DO k=1,klev
    3904                 DO i=1,klon
    3905                    IF (ptconv(i,k).AND.ptconvth(i,k)) THEN
    3906                       radocond(i,k)=radocond(i,k)+rnebcon(i,k)*clwcon(i,k)
    3907                       cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
    3908                    ELSE IF (ptconv(i,k)) THEN
    3909                       cldfra(i,k)=rnebcon(i,k)
    3910                       radocond(i,k)=rnebcon(i,k)*clwcon(i,k)
    3911                    ENDIF
    3912                 ENDDO
    3913              ENDDO
    3914 
    3915           ELSE IF (iflag_coupl==5) THEN
    3916              DO k=1,klev
    3917                 DO i=1,klon
    3918                    cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
    3919                    radocond(i,k)=radocond(i,k)+rnebcon(i,k)*clwcon(i,k)
    3920                 ENDDO
    3921              ENDDO
    3922 
    3923           ELSE
    3924 
    3925              ! Si on est sur un point touche par la convection
    3926              ! profonde et pas par les thermiques, on prend la
    3927              ! couverture nuageuse et l'eau nuageuse de la convection
    3928              ! profonde.
    3929 
    3930 
    3931              DO k=1,klev
    3932                 DO i=1,klon
    3933                    IF (ptconv(i,k).AND. .NOT.ptconvth(i,k)) THEN
    3934                       cldfra(i,k)=rnebcon(i,k)
    3935                       radocond(i,k)=rnebcon(i,k)*clwcon(i,k)
    3936                    ENDIF
    3937                 ENDDO
    3938              ENDDO
    3939 
    3940           ENDIF
    3941 
    3942        ELSE
    3943 
    3944           ! Ancienne version
    3945           cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
    3946           radocond(:,:)=radocond(:,:)+rnebcon(:,:)*clwcon(:,:)
    3947        ENDIF
    3948 
    3949     ENDIF
    3950 
    3951     ! 2. NUAGES STARTIFORMES
    3952     !
    3953     IF (ok_stratus) THEN
    3954        CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
    3955        DO k = 1, klev
    3956           DO i = 1, klon
    3957              IF (diafra(i,k).GT.cldfra(i,k)) THEN
    3958                 radocond(i,k) = dialiq(i,k)
    3959                 cldfra(i,k) = diafra(i,k)
    3960              ENDIF
    3961           ENDDO
    3962        ENDDO
    3963     ENDIF
    3964 
    3965 
    3966        ! Calculer les parametres optiques des nuages et quelques
    3967        ! parametres pour diagnostiques:
    3968        !
    3969        IF (aerosol_couple.AND.config_inca=='aero') THEN
    3970           mass_solu_aero(:,:)    = ccm(:,:,1)
    3971           mass_solu_aero_pi(:,:) = ccm(:,:,2)
    3972        ENDIF
    3973 
    3974        !Rajout appel a interface calcul proprietes optiques des nuages
    3975        CALL call_cloud_optics_prop(klon, klev, &
    3976             paprs, pplay, t_seri, radocond, radicefrac, cldfra, &
    3977             cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
    3978             flwp, fiwp, flwc, fiwc, ok_aie, &
    3979             mass_solu_aero, mass_solu_aero_pi, &
    3980             cldtaupi, distcltop, temp_cltop, re, fl, ref_liq, ref_ice, &
    3981             ref_liq_pi, ref_ice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, &
    3982             reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra,  &
    3983             zfice, dNovrN, ptconv, rnebcon, clwcon)
    3984        CALL call_cloud_optics_prop_post()
    3985 
    3986        
    3987        !IM betaCRF
    3988        !
    3989        cldtaurad   = cldtau
    3990        cldtaupirad = cldtaupi
    3991        cldemirad   = cldemi
    3992        cldfrarad   = cldfra
    3993 
    3994        !
    3995        IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
    3996             lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
    3997           !
    3998           ! global
    3999           !
    4000           DO k=1, klev
    4001              DO i=1, klon
    4002                 IF (pplay(i,k).GE.pfree) THEN
    4003                    beta(i,k) = beta_pbl
    4004                 ELSE
    4005                    beta(i,k) = beta_free
    4006                 ENDIF
    4007                 IF (mskocean_beta) THEN
    4008                    beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
    4009                 ENDIF
    4010                 cldtaurad(i,k)   = cldtau(i,k) * beta(i,k)
    4011                 cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
    4012                 cldemirad(i,k)   = cldemi(i,k) * beta(i,k)
    4013                 cldfrarad(i,k)   = cldfra(i,k) * beta(i,k)
    4014              ENDDO
    4015           ENDDO
    4016           !
    4017        ELSE
    4018           !
    4019           ! regional
    4020           !
    4021           DO k=1, klev
    4022              DO i=1,klon
    4023                 !
    4024                 IF (longitude_deg(i).ge.lon1_beta.AND. &
    4025                      longitude_deg(i).le.lon2_beta.AND. &
    4026                      latitude_deg(i).le.lat1_beta.AND.  &
    4027                      latitude_deg(i).ge.lat2_beta) THEN
    4028                    IF (pplay(i,k).GE.pfree) THEN
    4029                       beta(i,k) = beta_pbl
    4030                    ELSE
    4031                       beta(i,k) = beta_free
    4032                    ENDIF
    4033                    IF (mskocean_beta) THEN
    4034                       beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
    4035                    ENDIF
    4036                    cldtaurad(i,k)   = cldtau(i,k) * beta(i,k)
    4037                    cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
    4038                    cldemirad(i,k)   = cldemi(i,k) * beta(i,k)
    4039                    cldfrarad(i,k)   = cldfra(i,k) * beta(i,k)
    4040                 ENDIF
    4041                 !
    4042              ENDDO
    4043           ENDDO
    4044           !
    4045        ENDIF
    4046 
    40473802
    40483803
     
    40583813
    40593814          CALL chemtime(itap+itau_phy-1, date0, phys_tstep, itap)
    4060           CALL AEROSOL_METEO_CALC( &
     3815          CALL aerosol_meteo_calc( &
    40613816               calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
    40623817               prfl,psfl,pctsrf,cell_area, &
     
    41233878    ENDIF
    41243879
    4125    
    4126     !===============================================================================
    4127     ! Radiative scheme and associated aerosols
     3880    !===========================================================================
     3881    ! Aerosols
    41283882    !
     3883    ! Read aerosol forcing files to account for direct and first indirect
     3884    ! aerosols effects. Johannes Quaas, 27/11/2003
     3885    ! the radiative properties of aerosols depend on the radiative scheme used
     3886
    41293887    ! Note that the following routines are called every radpas time steps
    4130 
    41313888    IF (MOD(itaprad,radpas).EQ.0) THEN
    41323889
    4133        !
    4134        !jq - introduce the aerosol direct and first indirect radiative forcings
    4135        !jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
     3890       ! Tropospheric aerosols
     3891       !======================     
    41363892       IF (flag_aerosol .GT. 0) THEN
    4137           IF (iflag_rrtm .EQ. 0) THEN !--old radiation
    4138              IF (.NOT. aerosol_couple) THEN
    4139                 !
    4140                 CALL readaerosol_optic( &
     3893         
     3894          IF (iflag_rrtm .EQ. 0) THEN !--old radiative scheme
     3895
     3896             IF (.NOT. aerosol_couple) THEN !--
     3897               
     3898                     CALL readaerosol_optic( &
    41413899                     debut, flag_aerosol, itap, jD_cur-jD_ref, &
    41423900                     pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
     
    41443902                     tau_aero, piz_aero, cg_aero,  &
    41453903                     tausum_aero, tau3d_aero)
     3904
    41463905             ENDIF
    4147           ELSE IF (iflag_rrtm .EQ.1) THEN  ! RRTM radiation
     3906
     3907          ELSE IF (iflag_rrtm .EQ.1) THEN !--RRTM radive scheme
     3908
    41483909             IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
    41493910                abort_message='config_inca=aero et rrtm=1 impossible'
    41503911                CALL abort_physic(modname,abort_message,1)
     3912
    41513913             ELSE
    4152                 !
    41533914#ifdef CPP_RRTM
    41543915                IF (NSW.EQ.6) THEN
    4155                    !--new aerosol properties SW and LW
    4156                    !
     3916                   
    41573917                   IF (CPPKEY_DUST) THEN
    41583918                      !--SPL aerosol model
     
    41623922                           tausum_aero, tau3d_aero)
    41633923                   ELSE
    4164                       !--climatologies or INCA aerosols
     3924                      !--climatologies or inca
    41653925                      CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, ok_volcan, &
    41663926                           flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
     
    41723932
    41733933                   IF (flag_aerosol .EQ. 7) THEN
     3934
    41743935                      CALL macv2sp(pphis,pplay,paprs,longitude_deg,latitude_deg,  &
    41753936                           tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm,dNovrN)
     3937                   
    41763938                   ENDIF
    41773939
    4178                    !
    41793940                ELSE IF (NSW.EQ.2) THEN
    41803941                   !--for now we use the old aerosol properties
    4181                    !
    41823942                   CALL readaerosol_optic( &
    41833943                        debut, flag_aerosol, itap, jD_cur-jD_ref, &
     
    41863946                        tau_aero, piz_aero, cg_aero,  &
    41873947                        tausum_aero, tau3d_aero)
    4188                    !
    41893948                   !--natural aerosols
    41903949                   tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
     
    41953954                   piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
    41963955                   cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
    4197                    !
    41983956                   !--no LW optics
    41993957                   tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
    4200                    !
     3958               
    42013959                ELSE
    42023960                   abort_message='Only NSW=2 or 6 are possible with ' &
     
    42093967                CALL abort_physic(modname,abort_message,1)
    42103968#endif
    4211                 !
    42123969             ENDIF
    4213           ELSE IF (iflag_rrtm .EQ.2) THEN    ! ecrad RADIATION
     3970
     3971          ELSE IF (iflag_rrtm .EQ.2) THEN !--ECRAD radiative scheme
     3972
     3973             IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
     3974                abort_message='config_inca=aero et rrtm=2 impossible'
     3975                CALL abort_physic(modname,abort_message,1)
     3976             END IF                 
    42143977#ifdef CPP_ECRAD
    4215              !--climatologies or INCA aerosols
     3978             ! read climatologies or inca
    42163979             CALL readaerosol_optic_ecrad( debut, aerosol_couple, ok_alw, ok_volcan, &
    42173980                  flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
     
    42253988
    42263989       ELSE   !--flag_aerosol = 0
    4227           tausum_aero(:,:,:) = 0.
     3990          tausum_aero(:,:,:)  = 0.
    42283991          drytausum_aero(:,:) = 0.
    42293992          mass_solu_aero(:,:) = 0.
     
    42404003          ENDIF
    42414004       ENDIF
    4242        !
    4243        !--WMO criterion to determine tropopause
     4005       
     4006       IF (aerosol_couple.AND.config_inca=='aero') THEN
     4007          mass_solu_aero(:,:)    = ccm(:,:,1)
     4008          mass_solu_aero_pi(:,:) = ccm(:,:,2)
     4009       ENDIF
     4010
     4011
     4012
     4013       ! Stratospheric aerosols
     4014       !========================
     4015
     4016       !--WMO criterion to determine tropopause for masking the stratospheric aerosols fields
    42444017       CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
    42454018       !
     
    42484021       IF (flag_aerosol_strat.GT.0) THEN
    42494022          IF (prt_level .GE.10) THEN
    4250              PRINT *,'appel a readaerosolstrat', mth_cur
     4023             PRINT *,'call to readaerosolstrat', mth_cur
    42514024          ENDIF
    42524025          IF (iflag_rrtm.EQ.0) THEN
     
    42804053          tausum_aero(:,:,id_STRAT_phy) = 0.
    42814054       ENDIF
    4282        !
     4055
    42834056#ifdef CPP_RRTM
    42844057       IF (CPPKEY_STRATAER) THEN
     
    42914064       !--fin STRAT AEROSOL
    42924065       !
     4066   ENDIF !IF (MOD(itaprad,radpas).EQ.0)
     4067
     4068
     4069    !===============================================================================
     4070    ! from Clouds to Radiation (cloud2rad)
     4071    ! computation of cloud variables that are useful for the radiative transfer
     4072
     4073
     4074    DO k = 1, klev
     4075       DO i = 1, klon
     4076          cldfra(i,k) = rneb(i,k)
     4077          ! keep only liquid droplets in radocond if not liqice_in_radocond
     4078          IF (.NOT.liqice_in_radocond) radocond(i,k) = ql_seri(i,k)
     4079       ENDDO
     4080    ENDDO
     4081
     4082
     4083    ! Option to activate the radiative effect of blowing snow (ok_rad_bs)
     4084    ! makes sense only if the new large scale condensation scheme is active
     4085    ! with the ok_icefra_lscp flag active as well
     4086
     4087    IF (ok_bs .AND. ok_rad_bs) THEN
     4088       !IF (ok_icefra_lscp) THEN
     4089          DO k=1,klev
     4090             DO i=1,klon
     4091                radocond(i,k)=radocond(i,k)+qbs_seri(i,k)
     4092                radicefrac(i,k)=(radocond(i,k)*radicefrac(i,k)+qbs_seri(i,k))/(radocond(i,k))
     4093                qbsfra=min(qbs_seri(i,k)/qbst_bs,1.0)
     4094                cldfra(i,k)=max(cldfra(i,k),qbsfra)
     4095             ENDDO
     4096          ENDDO
     4097       !ELSE
     4098       !   WRITE(lunout,*)"PAY ATTENTION, you try to activate the radiative effect of blowing snow"
     4099       !   WRITE(lunout,*)"with ok_icefra_lscp=false. You must use ok_icefra_lscp=y and ok_new_lscp=y."
     4100       !   abort_message='inconsistency in cloud phase for blowing snow'
     4101       !   CALL abort_physic(modname,abort_message,1)
     4102       !ENDIF
     4103
     4104    ENDIF
     4105
     4106    IF (mydebug) THEN
     4107       CALL writefield_phy('u_seri',u_seri,nbp_lev)
     4108       CALL writefield_phy('v_seri',v_seri,nbp_lev)
     4109       CALL writefield_phy('t_seri',t_seri,nbp_lev)
     4110       CALL writefield_phy('q_seri',q_seri,nbp_lev)
     4111    ENDIF
     4112
     4113
     4114    ! 1. NUAGES CONVECTIFS
     4115    !     IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke
     4116    IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke
     4117       snow_tiedtke=0.
     4118       IF (iflag_cld_th.eq.-1) THEN
     4119          rain_tiedtke=rain_con
     4120       ELSE
     4121          rain_tiedtke=0.
     4122          DO k=1,klev
     4123             DO i=1,klon
     4124                IF (d_q_con(i,k).lt.0.) THEN
     4125                   rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
     4126                        *(paprs(i,k)-paprs(i,k+1))/rg
     4127                ENDIF
     4128             ENDDO
     4129          ENDDO
     4130       ENDIF
     4131
     4132       ! Nuages diagnostiques pour Tiedtke
     4133       CALL diagcld1(paprs,pplay, &
     4134            rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
     4135            diafra,dialiq)
     4136       DO k = 1, klev
     4137          DO i = 1, klon
     4138             IF (diafra(i,k).GT.cldfra(i,k)) THEN
     4139                radocond(i,k) = dialiq(i,k)
     4140                cldfra(i,k) = diafra(i,k)
     4141             ENDIF
     4142          ENDDO
     4143       ENDDO
     4144
     4145    ELSE IF (iflag_cld_th.ge.3) THEN
     4146       !  On prend pour les nuages convectifs le max du calcul de la
     4147       !  convection et du calcul du pas de temps precedent diminue d'un facteur
     4148       !  facttemps
     4149       facteur = pdtphys *facttemps
     4150       DO k=1,klev
     4151          DO i=1,klon
     4152             rnebcon(i,k)=rnebcon(i,k)*facteur
     4153             IF (rnebcon0(i,k)*clwcon0(i,k).GT.rnebcon(i,k)*clwcon(i,k)) THEN
     4154                rnebcon(i,k)=rnebcon0(i,k)
     4155                clwcon(i,k)=clwcon0(i,k)
     4156             ENDIF
     4157          ENDDO
     4158       ENDDO
     4159
     4160       !   On prend la somme des fractions nuageuses et des contenus en eau
     4161
     4162       IF (iflag_cld_th>=5) THEN
     4163
     4164          DO k=1,klev
     4165             ptconvth(:,k)=fm_therm(:,k+1)>0.
     4166          ENDDO
     4167
     4168          IF (iflag_coupl==4) THEN
     4169
     4170             ! Dans le cas iflag_coupl==4, on prend la somme des convertures
     4171             ! convectives et lsc dans la partie des thermiques
     4172             ! Le controle par iflag_coupl est peut etre provisoire.
     4173             DO k=1,klev
     4174                DO i=1,klon
     4175                   IF (ptconv(i,k).AND.ptconvth(i,k)) THEN
     4176                      radocond(i,k)=radocond(i,k)+rnebcon(i,k)*clwcon(i,k)
     4177                      cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
     4178                   ELSE IF (ptconv(i,k)) THEN
     4179                      cldfra(i,k)=rnebcon(i,k)
     4180                      radocond(i,k)=rnebcon(i,k)*clwcon(i,k)
     4181                   ENDIF
     4182                ENDDO
     4183             ENDDO
     4184
     4185          ELSE IF (iflag_coupl==5) THEN
     4186             DO k=1,klev
     4187                DO i=1,klon
     4188                   cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
     4189                   radocond(i,k)=radocond(i,k)+rnebcon(i,k)*clwcon(i,k)
     4190                ENDDO
     4191             ENDDO
     4192
     4193          ELSE
     4194
     4195             ! Si on est sur un point touche par la convection
     4196             ! profonde et pas par les thermiques, on prend la
     4197             ! couverture nuageuse et l'eau nuageuse de la convection
     4198             ! profonde.
     4199
     4200
     4201             DO k=1,klev
     4202                DO i=1,klon
     4203                   IF (ptconv(i,k).AND. .NOT.ptconvth(i,k)) THEN
     4204                      cldfra(i,k)=rnebcon(i,k)
     4205                      radocond(i,k)=rnebcon(i,k)*clwcon(i,k)
     4206                   ENDIF
     4207                ENDDO
     4208             ENDDO
     4209
     4210          ENDIF
     4211
     4212       ELSE
     4213
     4214          ! Ancienne version
     4215          cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
     4216          radocond(:,:)=radocond(:,:)+rnebcon(:,:)*clwcon(:,:)
     4217       ENDIF
     4218
     4219    ENDIF
     4220
     4221    ! 2. NUAGES STARTIFORMES
     4222    !
     4223    IF (ok_stratus) THEN
     4224       CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
     4225       DO k = 1, klev
     4226          DO i = 1, klon
     4227             IF (diafra(i,k).GT.cldfra(i,k)) THEN
     4228                radocond(i,k) = dialiq(i,k)
     4229                cldfra(i,k) = diafra(i,k)
     4230             ENDIF
     4231          ENDDO
     4232       ENDDO
     4233    ENDIF
     4234
     4235
     4236   ! Cloud optical parameters and some diagnostics
     4237
     4238   ! Note that the following routines are called every radpas time steps
     4239    IF (MOD(itaprad,radpas).EQ.0) THEN
     4240
     4241       CALL call_cloud_optics_prop(klon, klev, &
     4242            paprs, pplay, t_seri, radocond, radicefrac, cldfra, &
     4243            cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
     4244            flwp, fiwp, flwc, fiwc, ok_aie, &
     4245            mass_solu_aero, mass_solu_aero_pi, &
     4246            cldtaupi, distcltop, temp_cltop, re, fl, ref_liq, ref_ice, &
     4247            ref_liq_pi, ref_ice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, &
     4248            reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra,  &
     4249            zfice, dNovrN, ptconv, rnebcon, clwcon)
     4250       CALL call_cloud_optics_prop_post()
     4251
     4252       
     4253       !IM betaCRF
     4254       cldtaurad   = cldtau
     4255       cldtaupirad = cldtaupi
     4256       cldemirad   = cldemi
     4257       cldfrarad   = cldfra
     4258
     4259       IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
     4260            lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
     4261          !
     4262          ! global
     4263          !
     4264          DO k=1, klev
     4265             DO i=1, klon
     4266                IF (pplay(i,k).GE.pfree) THEN
     4267                   beta(i,k) = beta_pbl
     4268                ELSE
     4269                   beta(i,k) = beta_free
     4270                ENDIF
     4271                IF (mskocean_beta) THEN
     4272                   beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
     4273                ENDIF
     4274                cldtaurad(i,k)   = cldtau(i,k) * beta(i,k)
     4275                cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
     4276                cldemirad(i,k)   = cldemi(i,k) * beta(i,k)
     4277                cldfrarad(i,k)   = cldfra(i,k) * beta(i,k)
     4278             ENDDO
     4279          ENDDO
     4280       ELSE
     4281          !
     4282          ! regional
     4283          !
     4284          DO k=1, klev
     4285             DO i=1,klon
     4286                !
     4287                IF (longitude_deg(i).ge.lon1_beta.AND. &
     4288                     longitude_deg(i).le.lon2_beta.AND. &
     4289                     latitude_deg(i).le.lat1_beta.AND.  &
     4290                     latitude_deg(i).ge.lat2_beta) THEN
     4291                   IF (pplay(i,k).GE.pfree) THEN
     4292                      beta(i,k) = beta_pbl
     4293                   ELSE
     4294                      beta(i,k) = beta_free
     4295                   ENDIF
     4296                   IF (mskocean_beta) THEN
     4297                      beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
     4298                   ENDIF
     4299                   cldtaurad(i,k)   = cldtau(i,k) * beta(i,k)
     4300                   cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
     4301                   cldemirad(i,k)   = cldemi(i,k) * beta(i,k)
     4302                   cldfrarad(i,k)   = cldfra(i,k) * beta(i,k)
     4303                ENDIF
     4304             ENDDO
     4305          ENDDO
     4306       ENDIF
     4307   
     4308    END IF ! radpas time steps
     4309   
     4310    !===============================================================================
     4311    ! Radiation scheme
     4312    !
     4313    ! Note that the following routines are called every radpas time steps
     4314
     4315    IF (MOD(itaprad,radpas).EQ.0) THEN
    42934316
    42944317       !lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek
Note: See TracChangeset for help on using the changeset viewer.