Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 6146) +++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 6148) @@ -3756,4 +3756,7 @@ call surf_wind(klon,nsurfwind,zu10m,zv10m,wake_s,wake_Cstar,zustar,ale_bl,surf_wind_value,surf_wind_proba) + + + !=========================================================================== ! Large scale condensation and precipitation @@ -3797,252 +3800,4 @@ dqsmelt, dqsfreez, zx_rh, zx_rhl, zx_rhi) - !=============================================================================== - ! from Clouds to Radiation (cloud2rad) - ! computation of cloud variables that are useful for the radiative transfer - - - DO k = 1, klev - DO i = 1, klon - cldfra(i,k) = rneb(i,k) - ! keep only liquid droplets in radocond if not liqice_in_radocond - IF (.NOT.liqice_in_radocond) radocond(i,k) = ql_seri(i,k) - ENDDO - ENDDO - - - ! Option to activate the radiative effect of blowing snow (ok_rad_bs) - ! makes sense only if the new large scale condensation scheme is active - ! with the ok_icefra_lscp flag active as well - - IF (ok_bs .AND. ok_rad_bs) THEN - !IF (ok_icefra_lscp) THEN - DO k=1,klev - DO i=1,klon - radocond(i,k)=radocond(i,k)+qbs_seri(i,k) - radicefrac(i,k)=(radocond(i,k)*radicefrac(i,k)+qbs_seri(i,k))/(radocond(i,k)) - qbsfra=min(qbs_seri(i,k)/qbst_bs,1.0) - cldfra(i,k)=max(cldfra(i,k),qbsfra) - ENDDO - ENDDO - !ELSE - ! WRITE(lunout,*)"PAY ATTENTION, you try to activate the radiative effect of blowing snow" - ! WRITE(lunout,*)"with ok_icefra_lscp=false. You must use ok_icefra_lscp=y and ok_new_lscp=y." - ! abort_message='inconsistency in cloud phase for blowing snow' - ! CALL abort_physic(modname,abort_message,1) - !ENDIF - - ENDIF - - IF (mydebug) THEN - CALL writefield_phy('u_seri',u_seri,nbp_lev) - CALL writefield_phy('v_seri',v_seri,nbp_lev) - CALL writefield_phy('t_seri',t_seri,nbp_lev) - CALL writefield_phy('q_seri',q_seri,nbp_lev) - ENDIF - - - ! 1. NUAGES CONVECTIFS - ! IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke - IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke - snow_tiedtke=0. - IF (iflag_cld_th.eq.-1) THEN - rain_tiedtke=rain_con - ELSE - rain_tiedtke=0. - DO k=1,klev - DO i=1,klon - IF (d_q_con(i,k).lt.0.) THEN - rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys & - *(paprs(i,k)-paprs(i,k+1))/rg - ENDIF - ENDDO - ENDDO - ENDIF - - ! Nuages diagnostiques pour Tiedtke - CALL diagcld1(paprs,pplay, & - rain_tiedtke,snow_tiedtke,ibas_con,itop_con, & - diafra,dialiq) - DO k = 1, klev - DO i = 1, klon - IF (diafra(i,k).GT.cldfra(i,k)) THEN - radocond(i,k) = dialiq(i,k) - cldfra(i,k) = diafra(i,k) - ENDIF - ENDDO - ENDDO - - ELSE IF (iflag_cld_th.ge.3) THEN - ! On prend pour les nuages convectifs le max du calcul de la - ! convection et du calcul du pas de temps precedent diminue d'un facteur - ! facttemps - facteur = pdtphys *facttemps - DO k=1,klev - DO i=1,klon - rnebcon(i,k)=rnebcon(i,k)*facteur - IF (rnebcon0(i,k)*clwcon0(i,k).GT.rnebcon(i,k)*clwcon(i,k)) THEN - rnebcon(i,k)=rnebcon0(i,k) - clwcon(i,k)=clwcon0(i,k) - ENDIF - ENDDO - ENDDO - - ! On prend la somme des fractions nuageuses et des contenus en eau - - IF (iflag_cld_th>=5) THEN - - DO k=1,klev - ptconvth(:,k)=fm_therm(:,k+1)>0. - ENDDO - - IF (iflag_coupl==4) THEN - - ! Dans le cas iflag_coupl==4, on prend la somme des convertures - ! convectives et lsc dans la partie des thermiques - ! Le controle par iflag_coupl est peut etre provisoire. - DO k=1,klev - DO i=1,klon - IF (ptconv(i,k).AND.ptconvth(i,k)) THEN - radocond(i,k)=radocond(i,k)+rnebcon(i,k)*clwcon(i,k) - cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.) - ELSE IF (ptconv(i,k)) THEN - cldfra(i,k)=rnebcon(i,k) - radocond(i,k)=rnebcon(i,k)*clwcon(i,k) - ENDIF - ENDDO - ENDDO - - ELSE IF (iflag_coupl==5) THEN - DO k=1,klev - DO i=1,klon - cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.) - radocond(i,k)=radocond(i,k)+rnebcon(i,k)*clwcon(i,k) - ENDDO - ENDDO - - ELSE - - ! Si on est sur un point touche par la convection - ! profonde et pas par les thermiques, on prend la - ! couverture nuageuse et l'eau nuageuse de la convection - ! profonde. - - - DO k=1,klev - DO i=1,klon - IF (ptconv(i,k).AND. .NOT.ptconvth(i,k)) THEN - cldfra(i,k)=rnebcon(i,k) - radocond(i,k)=rnebcon(i,k)*clwcon(i,k) - ENDIF - ENDDO - ENDDO - - ENDIF - - ELSE - - ! Ancienne version - cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.) - radocond(:,:)=radocond(:,:)+rnebcon(:,:)*clwcon(:,:) - ENDIF - - ENDIF - - ! 2. NUAGES STARTIFORMES - ! - IF (ok_stratus) THEN - CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq) - DO k = 1, klev - DO i = 1, klon - IF (diafra(i,k).GT.cldfra(i,k)) THEN - radocond(i,k) = dialiq(i,k) - cldfra(i,k) = diafra(i,k) - ENDIF - ENDDO - ENDDO - ENDIF - - - ! Calculer les parametres optiques des nuages et quelques - ! parametres pour diagnostiques: - ! - IF (aerosol_couple.AND.config_inca=='aero') THEN - mass_solu_aero(:,:) = ccm(:,:,1) - mass_solu_aero_pi(:,:) = ccm(:,:,2) - ENDIF - - !Rajout appel a interface calcul proprietes optiques des nuages - CALL call_cloud_optics_prop(klon, klev, & - paprs, pplay, t_seri, radocond, radicefrac, cldfra, & - cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, & - flwp, fiwp, flwc, fiwc, ok_aie, & - mass_solu_aero, mass_solu_aero_pi, & - cldtaupi, distcltop, temp_cltop, re, fl, ref_liq, ref_ice, & - ref_liq_pi, ref_ice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, & - reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra, & - zfice, dNovrN, ptconv, rnebcon, clwcon) - CALL call_cloud_optics_prop_post() - - - !IM betaCRF - ! - cldtaurad = cldtau - cldtaupirad = cldtaupi - cldemirad = cldemi - cldfrarad = cldfra - - ! - IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. & - lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN - ! - ! global - ! - DO k=1, klev - DO i=1, klon - IF (pplay(i,k).GE.pfree) THEN - beta(i,k) = beta_pbl - ELSE - beta(i,k) = beta_free - ENDIF - IF (mskocean_beta) THEN - beta(i,k) = beta(i,k) * pctsrf(i,is_oce) - ENDIF - cldtaurad(i,k) = cldtau(i,k) * beta(i,k) - cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k) - cldemirad(i,k) = cldemi(i,k) * beta(i,k) - cldfrarad(i,k) = cldfra(i,k) * beta(i,k) - ENDDO - ENDDO - ! - ELSE - ! - ! regional - ! - DO k=1, klev - DO i=1,klon - ! - IF (longitude_deg(i).ge.lon1_beta.AND. & - longitude_deg(i).le.lon2_beta.AND. & - latitude_deg(i).le.lat1_beta.AND. & - latitude_deg(i).ge.lat2_beta) THEN - IF (pplay(i,k).GE.pfree) THEN - beta(i,k) = beta_pbl - ELSE - beta(i,k) = beta_free - ENDIF - IF (mskocean_beta) THEN - beta(i,k) = beta(i,k) * pctsrf(i,is_oce) - ENDIF - cldtaurad(i,k) = cldtau(i,k) * beta(i,k) - cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k) - cldemirad(i,k) = cldemi(i,k) * beta(i,k) - cldfrarad(i,k) = cldfra(i,k) * beta(i,k) - ENDIF - ! - ENDDO - ENDDO - ! - ENDIF - @@ -4058,5 +3813,5 @@ CALL chemtime(itap+itau_phy-1, date0, phys_tstep, itap) - CALL AEROSOL_METEO_CALC( & + CALL aerosol_meteo_calc( & calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, & prfl,psfl,pctsrf,cell_area, & @@ -4123,20 +3878,23 @@ ENDIF - - !=============================================================================== - ! Radiative scheme and associated aerosols + !=========================================================================== + ! Aerosols ! + ! Read aerosol forcing files to account for direct and first indirect + ! aerosols effects. Johannes Quaas, 27/11/2003 + ! the radiative properties of aerosols depend on the radiative scheme used + ! Note that the following routines are called every radpas time steps - IF (MOD(itaprad,radpas).EQ.0) THEN - ! - !jq - introduce the aerosol direct and first indirect radiative forcings - !jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr) + ! Tropospheric aerosols + !====================== IF (flag_aerosol .GT. 0) THEN - IF (iflag_rrtm .EQ. 0) THEN !--old radiation - IF (.NOT. aerosol_couple) THEN - ! - CALL readaerosol_optic( & + + IF (iflag_rrtm .EQ. 0) THEN !--old radiative scheme + + IF (.NOT. aerosol_couple) THEN !-- + + CALL readaerosol_optic( & debut, flag_aerosol, itap, jD_cur-jD_ref, & pdtphys, pplay, paprs, t_seri, rhcl, presnivs, & @@ -4144,15 +3902,17 @@ tau_aero, piz_aero, cg_aero, & tausum_aero, tau3d_aero) + ENDIF - ELSE IF (iflag_rrtm .EQ.1) THEN ! RRTM radiation + + ELSE IF (iflag_rrtm .EQ.1) THEN !--RRTM radive scheme + IF (aerosol_couple .AND. config_inca == 'aero' ) THEN abort_message='config_inca=aero et rrtm=1 impossible' CALL abort_physic(modname,abort_message,1) + ELSE - ! #ifdef CPP_RRTM IF (NSW.EQ.6) THEN - !--new aerosol properties SW and LW - ! + IF (CPPKEY_DUST) THEN !--SPL aerosol model @@ -4162,5 +3922,5 @@ tausum_aero, tau3d_aero) ELSE - !--climatologies or INCA aerosols + !--climatologies or inca CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, ok_volcan, & flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, & @@ -4172,12 +3932,12 @@ IF (flag_aerosol .EQ. 7) THEN + CALL macv2sp(pphis,pplay,paprs,longitude_deg,latitude_deg, & tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm,dNovrN) + ENDIF - ! ELSE IF (NSW.EQ.2) THEN !--for now we use the old aerosol properties - ! CALL readaerosol_optic( & debut, flag_aerosol, itap, jD_cur-jD_ref, & @@ -4186,5 +3946,4 @@ tau_aero, piz_aero, cg_aero, & tausum_aero, tau3d_aero) - ! !--natural aerosols tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:) @@ -4195,8 +3954,7 @@ piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:) cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:) - ! !--no LW optics tau_aero_lw_rrtm(:,:,:,:) = 1.e-15 - ! + ELSE abort_message='Only NSW=2 or 6 are possible with ' & @@ -4209,9 +3967,14 @@ CALL abort_physic(modname,abort_message,1) #endif - ! ENDIF - ELSE IF (iflag_rrtm .EQ.2) THEN ! ecrad RADIATION + + ELSE IF (iflag_rrtm .EQ.2) THEN !--ECRAD radiative scheme + + IF (aerosol_couple .AND. config_inca == 'aero' ) THEN + abort_message='config_inca=aero et rrtm=2 impossible' + CALL abort_physic(modname,abort_message,1) + END IF #ifdef CPP_ECRAD - !--climatologies or INCA aerosols + ! read climatologies or inca CALL readaerosol_optic_ecrad( debut, aerosol_couple, ok_alw, ok_volcan, & flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, & @@ -4225,5 +3988,5 @@ ELSE !--flag_aerosol = 0 - tausum_aero(:,:,:) = 0. + tausum_aero(:,:,:) = 0. drytausum_aero(:,:) = 0. mass_solu_aero(:,:) = 0. @@ -4240,6 +4003,16 @@ ENDIF ENDIF - ! - !--WMO criterion to determine tropopause + + IF (aerosol_couple.AND.config_inca=='aero') THEN + mass_solu_aero(:,:) = ccm(:,:,1) + mass_solu_aero_pi(:,:) = ccm(:,:,2) + ENDIF + + + + ! Stratospheric aerosols + !======================== + + !--WMO criterion to determine tropopause for masking the stratospheric aerosols fields CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg) ! @@ -4248,5 +4021,5 @@ IF (flag_aerosol_strat.GT.0) THEN IF (prt_level .GE.10) THEN - PRINT *,'appel a readaerosolstrat', mth_cur + PRINT *,'call to readaerosolstrat', mth_cur ENDIF IF (iflag_rrtm.EQ.0) THEN @@ -4280,5 +4053,5 @@ tausum_aero(:,:,id_STRAT_phy) = 0. ENDIF - ! + #ifdef CPP_RRTM IF (CPPKEY_STRATAER) THEN @@ -4291,4 +4064,254 @@ !--fin STRAT AEROSOL ! + ENDIF !IF (MOD(itaprad,radpas).EQ.0) + + + !=============================================================================== + ! from Clouds to Radiation (cloud2rad) + ! computation of cloud variables that are useful for the radiative transfer + + + DO k = 1, klev + DO i = 1, klon + cldfra(i,k) = rneb(i,k) + ! keep only liquid droplets in radocond if not liqice_in_radocond + IF (.NOT.liqice_in_radocond) radocond(i,k) = ql_seri(i,k) + ENDDO + ENDDO + + + ! Option to activate the radiative effect of blowing snow (ok_rad_bs) + ! makes sense only if the new large scale condensation scheme is active + ! with the ok_icefra_lscp flag active as well + + IF (ok_bs .AND. ok_rad_bs) THEN + !IF (ok_icefra_lscp) THEN + DO k=1,klev + DO i=1,klon + radocond(i,k)=radocond(i,k)+qbs_seri(i,k) + radicefrac(i,k)=(radocond(i,k)*radicefrac(i,k)+qbs_seri(i,k))/(radocond(i,k)) + qbsfra=min(qbs_seri(i,k)/qbst_bs,1.0) + cldfra(i,k)=max(cldfra(i,k),qbsfra) + ENDDO + ENDDO + !ELSE + ! WRITE(lunout,*)"PAY ATTENTION, you try to activate the radiative effect of blowing snow" + ! WRITE(lunout,*)"with ok_icefra_lscp=false. You must use ok_icefra_lscp=y and ok_new_lscp=y." + ! abort_message='inconsistency in cloud phase for blowing snow' + ! CALL abort_physic(modname,abort_message,1) + !ENDIF + + ENDIF + + IF (mydebug) THEN + CALL writefield_phy('u_seri',u_seri,nbp_lev) + CALL writefield_phy('v_seri',v_seri,nbp_lev) + CALL writefield_phy('t_seri',t_seri,nbp_lev) + CALL writefield_phy('q_seri',q_seri,nbp_lev) + ENDIF + + + ! 1. NUAGES CONVECTIFS + ! IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke + IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke + snow_tiedtke=0. + IF (iflag_cld_th.eq.-1) THEN + rain_tiedtke=rain_con + ELSE + rain_tiedtke=0. + DO k=1,klev + DO i=1,klon + IF (d_q_con(i,k).lt.0.) THEN + rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys & + *(paprs(i,k)-paprs(i,k+1))/rg + ENDIF + ENDDO + ENDDO + ENDIF + + ! Nuages diagnostiques pour Tiedtke + CALL diagcld1(paprs,pplay, & + rain_tiedtke,snow_tiedtke,ibas_con,itop_con, & + diafra,dialiq) + DO k = 1, klev + DO i = 1, klon + IF (diafra(i,k).GT.cldfra(i,k)) THEN + radocond(i,k) = dialiq(i,k) + cldfra(i,k) = diafra(i,k) + ENDIF + ENDDO + ENDDO + + ELSE IF (iflag_cld_th.ge.3) THEN + ! On prend pour les nuages convectifs le max du calcul de la + ! convection et du calcul du pas de temps precedent diminue d'un facteur + ! facttemps + facteur = pdtphys *facttemps + DO k=1,klev + DO i=1,klon + rnebcon(i,k)=rnebcon(i,k)*facteur + IF (rnebcon0(i,k)*clwcon0(i,k).GT.rnebcon(i,k)*clwcon(i,k)) THEN + rnebcon(i,k)=rnebcon0(i,k) + clwcon(i,k)=clwcon0(i,k) + ENDIF + ENDDO + ENDDO + + ! On prend la somme des fractions nuageuses et des contenus en eau + + IF (iflag_cld_th>=5) THEN + + DO k=1,klev + ptconvth(:,k)=fm_therm(:,k+1)>0. + ENDDO + + IF (iflag_coupl==4) THEN + + ! Dans le cas iflag_coupl==4, on prend la somme des convertures + ! convectives et lsc dans la partie des thermiques + ! Le controle par iflag_coupl est peut etre provisoire. + DO k=1,klev + DO i=1,klon + IF (ptconv(i,k).AND.ptconvth(i,k)) THEN + radocond(i,k)=radocond(i,k)+rnebcon(i,k)*clwcon(i,k) + cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.) + ELSE IF (ptconv(i,k)) THEN + cldfra(i,k)=rnebcon(i,k) + radocond(i,k)=rnebcon(i,k)*clwcon(i,k) + ENDIF + ENDDO + ENDDO + + ELSE IF (iflag_coupl==5) THEN + DO k=1,klev + DO i=1,klon + cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.) + radocond(i,k)=radocond(i,k)+rnebcon(i,k)*clwcon(i,k) + ENDDO + ENDDO + + ELSE + + ! Si on est sur un point touche par la convection + ! profonde et pas par les thermiques, on prend la + ! couverture nuageuse et l'eau nuageuse de la convection + ! profonde. + + + DO k=1,klev + DO i=1,klon + IF (ptconv(i,k).AND. .NOT.ptconvth(i,k)) THEN + cldfra(i,k)=rnebcon(i,k) + radocond(i,k)=rnebcon(i,k)*clwcon(i,k) + ENDIF + ENDDO + ENDDO + + ENDIF + + ELSE + + ! Ancienne version + cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.) + radocond(:,:)=radocond(:,:)+rnebcon(:,:)*clwcon(:,:) + ENDIF + + ENDIF + + ! 2. NUAGES STARTIFORMES + ! + IF (ok_stratus) THEN + CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq) + DO k = 1, klev + DO i = 1, klon + IF (diafra(i,k).GT.cldfra(i,k)) THEN + radocond(i,k) = dialiq(i,k) + cldfra(i,k) = diafra(i,k) + ENDIF + ENDDO + ENDDO + ENDIF + + + ! Cloud optical parameters and some diagnostics + + ! Note that the following routines are called every radpas time steps + IF (MOD(itaprad,radpas).EQ.0) THEN + + CALL call_cloud_optics_prop(klon, klev, & + paprs, pplay, t_seri, radocond, radicefrac, cldfra, & + cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, & + flwp, fiwp, flwc, fiwc, ok_aie, & + mass_solu_aero, mass_solu_aero_pi, & + cldtaupi, distcltop, temp_cltop, re, fl, ref_liq, ref_ice, & + ref_liq_pi, ref_ice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, & + reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra, & + zfice, dNovrN, ptconv, rnebcon, clwcon) + CALL call_cloud_optics_prop_post() + + + !IM betaCRF + cldtaurad = cldtau + cldtaupirad = cldtaupi + cldemirad = cldemi + cldfrarad = cldfra + + IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. & + lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN + ! + ! global + ! + DO k=1, klev + DO i=1, klon + IF (pplay(i,k).GE.pfree) THEN + beta(i,k) = beta_pbl + ELSE + beta(i,k) = beta_free + ENDIF + IF (mskocean_beta) THEN + beta(i,k) = beta(i,k) * pctsrf(i,is_oce) + ENDIF + cldtaurad(i,k) = cldtau(i,k) * beta(i,k) + cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k) + cldemirad(i,k) = cldemi(i,k) * beta(i,k) + cldfrarad(i,k) = cldfra(i,k) * beta(i,k) + ENDDO + ENDDO + ELSE + ! + ! regional + ! + DO k=1, klev + DO i=1,klon + ! + IF (longitude_deg(i).ge.lon1_beta.AND. & + longitude_deg(i).le.lon2_beta.AND. & + latitude_deg(i).le.lat1_beta.AND. & + latitude_deg(i).ge.lat2_beta) THEN + IF (pplay(i,k).GE.pfree) THEN + beta(i,k) = beta_pbl + ELSE + beta(i,k) = beta_free + ENDIF + IF (mskocean_beta) THEN + beta(i,k) = beta(i,k) * pctsrf(i,is_oce) + ENDIF + cldtaurad(i,k) = cldtau(i,k) * beta(i,k) + cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k) + cldemirad(i,k) = cldemi(i,k) * beta(i,k) + cldfrarad(i,k) = cldfra(i,k) * beta(i,k) + ENDIF + ENDDO + ENDDO + ENDIF + + END IF ! radpas time steps + + !=============================================================================== + ! Radiation scheme + ! + ! Note that the following routines are called every radpas time steps + + IF (MOD(itaprad,radpas).EQ.0) THEN !lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek