Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 6144) +++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 6145) @@ -3802,5 +3802,4 @@ - DO k = 1, klev DO i = 1, klon @@ -3929,6 +3928,4 @@ ! profonde. - !IM/FH: 2011/02/23 - ! definition des points sur lesquels ls thermiques sont actifs DO k=1,klev @@ -3966,248 +3963,4 @@ ENDIF - !=============================================================================== - ! Interactive chemistry through coupling with INCA or REPROBUS chemistry models - ! - - IF (ANY(type_trac == ['inca','inco'])) THEN ! ModThL - IF (CPPKEY_INCA) THEN - CALL VTe(VTphysiq) - CALL VTb(VTinca) - calday = REAL(days_elapsed + 1) + jH_cur - - CALL chemtime(itap+itau_phy-1, date0, phys_tstep, itap) - CALL AEROSOL_METEO_CALC( & - calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, & - prfl,psfl,pctsrf,cell_area, & - latitude_deg,longitude_deg,u10m,v10m) - - zxsnow_dummy(:) = 0.0 - ! INCA needs a cloud fraction that is not necessarily that - ! for radiation. Here we provide a cloud fraction calculated - ! the same manner as that in LMDZ5, LMDZ6 and LMDZ7 - DO k=1,klev - DO i=1,klon - cldfra_inca(i,k)=min(rneb(i,k)+rnebcon(i,k),1.) - ENDDO - ENDDO - - - CALL chemhook_begin (calday, & - days_elapsed+1, & - jH_cur, & - pctsrf(1,1), & - latitude_deg, & - longitude_deg, & - cell_area, & - paprs, & - pplay, & - coefh(1:klon,1:klev,is_ave), & - pphi, & - t_seri, & - u, & - v, & - rot, & - wo(:, :, 1), & - q_seri, & - zxtsol, & - zt2m, & - zxsnow_dummy, & - solsw, & - albsol1, & - rain_fall, & - snow_fall, & - itop_con, & - ibas_con, & - cldfra_inca, & - nbp_lon, & - nbp_lat-1, & - tr_seri(:,:,1+nqCO2:nbtr), & - ftsol, & - paprs, & - cdragh, & - cdragm, & - pctsrf, & - pdtphys, & - itap) - - CALL VTe(VTinca) - CALL VTb(VTphysiq) - END IF - ENDIF !type_trac = inca or inco - - IF (type_trac == 'repr') THEN - IF (CPPKEY_REPROBUS) THEN - CALL chemtime_rep(itap+itau_phy-1, date0, phys_tstep, itap) - END IF - ENDIF - - - !=============================================================================== - ! Radiative scheme and associated aerosols - ! - ! 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) - IF (flag_aerosol .GT. 0) THEN - IF (iflag_rrtm .EQ. 0) THEN !--old radiation - IF (.NOT. aerosol_couple) THEN - ! - CALL readaerosol_optic( & - debut, flag_aerosol, itap, jD_cur-jD_ref, & - pdtphys, pplay, paprs, t_seri, rhcl, presnivs, & - mass_solu_aero, mass_solu_aero_pi, & - tau_aero, piz_aero, cg_aero, & - tausum_aero, tau3d_aero) - ENDIF - ELSE IF (iflag_rrtm .EQ.1) THEN ! RRTM radiation - 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 - CALL splaerosol_optic_rrtm( ok_alw, pplay, paprs, t_seri, rhcl, & - tr_seri, mass_solu_aero, mass_solu_aero_pi, & - tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, & - tausum_aero, tau3d_aero) - ELSE - !--climatologies or INCA aerosols - CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, ok_volcan, & - flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, & - pdtphys, pplay, paprs, t_seri, rhcl, presnivs, & - tr_seri, mass_solu_aero, mass_solu_aero_pi, & - tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, & - tausum_aero, drytausum_aero, tau3d_aero) - END IF - - 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, & - pdtphys, pplay, paprs, t_seri, rhcl, presnivs, & - mass_solu_aero, mass_solu_aero_pi, & - tau_aero, piz_aero, cg_aero, & - tausum_aero, tau3d_aero) - ! - !--natural aerosols - tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:) - piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:) - cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:) - !--all aerosols - tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:) - 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 ' & - // 'aerosols and iflag_rrtm=1' - CALL abort_physic(modname,abort_message,1) - ENDIF -#else - abort_message='You should compile with -rrtm if running ' & - // 'with iflag_rrtm=1' - CALL abort_physic(modname,abort_message,1) -#endif - ! - ENDIF - ELSE IF (iflag_rrtm .EQ.2) THEN ! ecrad RADIATION -#ifdef CPP_ECRAD - !--climatologies or INCA aerosols - CALL readaerosol_optic_ecrad( debut, aerosol_couple, ok_alw, ok_volcan, & - flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, & - pdtphys, pplay, paprs, t_seri, rhcl, presnivs, & - tr_seri, mass_solu_aero, mass_solu_aero_pi, m_allaer) -#else - abort_message='You should compile with -rad ecrad if running with iflag_rrtm=2' - CALL abort_physic(modname,abort_message,1) -#endif - ENDIF - - ELSE !--flag_aerosol = 0 - tausum_aero(:,:,:) = 0. - drytausum_aero(:,:) = 0. - mass_solu_aero(:,:) = 0. - mass_solu_aero_pi(:,:) = 0. - IF (iflag_rrtm .EQ. 0) THEN !--old radiation - tau_aero(:,:,:,:) = 1.e-15 - piz_aero(:,:,:,:) = 1. - cg_aero(:,:,:,:) = 0. - ELSE - tau_aero_sw_rrtm(:,:,:,:) = 1.e-15 - tau_aero_lw_rrtm(:,:,:,:) = 1.e-15 - piz_aero_sw_rrtm(:,:,:,:) = 1.0 - cg_aero_sw_rrtm(:,:,:,:) = 0.0 - ENDIF - ENDIF - ! - !--WMO criterion to determine tropopause - CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg) - ! - !--STRAT AEROSOL - !--updates tausum_aero,tau_aero,piz_aero,cg_aero - IF (flag_aerosol_strat.GT.0) THEN - IF (prt_level .GE.10) THEN - PRINT *,'appel a readaerosolstrat', mth_cur - ENDIF - IF (iflag_rrtm.EQ.0) THEN - IF (flag_aerosol_strat.EQ.1) THEN - CALL readaerosolstrato(debut) - ELSE - abort_message='flag_aerosol_strat must equal 1 for rrtm=0' - CALL abort_physic(modname,abort_message,1) - ENDIF - ELSE -#ifdef CPP_RRTM - IF (.NOT. CPPKEY_STRATAER) THEN - !--prescribed strat aerosols - !--only in the case of non-interactive strat aerosols - IF (flag_aerosol_strat.EQ.1) THEN - CALL readaerosolstrato1_rrtm(debut) - ELSEIF (flag_aerosol_strat.EQ.2) THEN - CALL readaerosolstrato2_rrtm(debut, ok_volcan) - ELSE - abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1' - CALL abort_physic(modname,abort_message,1) - ENDIF - END IF -#else - abort_message='You should compile with -rrtm if running ' & - // 'with iflag_rrtm=1' - CALL abort_physic(modname,abort_message,1) -#endif - ENDIF - ELSE - tausum_aero(:,:,id_STRAT_phy) = 0. - ENDIF - ! -#ifdef CPP_RRTM - IF (CPPKEY_STRATAER) THEN - !--compute stratospheric mask - CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg) - !--interactive strat aerosols - CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut) - END IF -#endif - !--fin STRAT AEROSOL - ! ! Calculer les parametres optiques des nuages et quelques @@ -4231,5 +3984,5 @@ CALL call_cloud_optics_prop_post() - ! + !IM betaCRF ! @@ -4291,4 +4044,251 @@ ! ENDIF + + + + !=============================================================================== + ! Interactive chemistry through coupling with INCA or REPROBUS chemistry models + ! + + IF (ANY(type_trac == ['inca','inco'])) THEN ! ModThL + IF (CPPKEY_INCA) THEN + CALL VTe(VTphysiq) + CALL VTb(VTinca) + calday = REAL(days_elapsed + 1) + jH_cur + + CALL chemtime(itap+itau_phy-1, date0, phys_tstep, itap) + CALL AEROSOL_METEO_CALC( & + calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, & + prfl,psfl,pctsrf,cell_area, & + latitude_deg,longitude_deg,u10m,v10m) + + zxsnow_dummy(:) = 0.0 + ! INCA needs a cloud fraction that is not necessarily that + ! for radiation. Here we provide a cloud fraction calculated + ! the same manner as that in LMDZ5, LMDZ6 and LMDZ7 + DO k=1,klev + DO i=1,klon + cldfra_inca(i,k)=min(rneb(i,k)+rnebcon(i,k),1.) + ENDDO + ENDDO + + + CALL chemhook_begin (calday, & + days_elapsed+1, & + jH_cur, & + pctsrf(1,1), & + latitude_deg, & + longitude_deg, & + cell_area, & + paprs, & + pplay, & + coefh(1:klon,1:klev,is_ave), & + pphi, & + t_seri, & + u, & + v, & + rot, & + wo(:, :, 1), & + q_seri, & + zxtsol, & + zt2m, & + zxsnow_dummy, & + solsw, & + albsol1, & + rain_fall, & + snow_fall, & + itop_con, & + ibas_con, & + cldfra_inca, & + nbp_lon, & + nbp_lat-1, & + tr_seri(:,:,1+nqCO2:nbtr), & + ftsol, & + paprs, & + cdragh, & + cdragm, & + pctsrf, & + pdtphys, & + itap) + + CALL VTe(VTinca) + CALL VTb(VTphysiq) + END IF + ENDIF !type_trac = inca or inco + + IF (type_trac == 'repr') THEN + IF (CPPKEY_REPROBUS) THEN + CALL chemtime_rep(itap+itau_phy-1, date0, phys_tstep, itap) + END IF + ENDIF + + + !=============================================================================== + ! Radiative scheme and associated aerosols + ! + ! 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) + IF (flag_aerosol .GT. 0) THEN + IF (iflag_rrtm .EQ. 0) THEN !--old radiation + IF (.NOT. aerosol_couple) THEN + ! + CALL readaerosol_optic( & + debut, flag_aerosol, itap, jD_cur-jD_ref, & + pdtphys, pplay, paprs, t_seri, rhcl, presnivs, & + mass_solu_aero, mass_solu_aero_pi, & + tau_aero, piz_aero, cg_aero, & + tausum_aero, tau3d_aero) + ENDIF + ELSE IF (iflag_rrtm .EQ.1) THEN ! RRTM radiation + 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 + CALL splaerosol_optic_rrtm( ok_alw, pplay, paprs, t_seri, rhcl, & + tr_seri, mass_solu_aero, mass_solu_aero_pi, & + tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, & + tausum_aero, tau3d_aero) + ELSE + !--climatologies or INCA aerosols + CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, ok_volcan, & + flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, & + pdtphys, pplay, paprs, t_seri, rhcl, presnivs, & + tr_seri, mass_solu_aero, mass_solu_aero_pi, & + tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, & + tausum_aero, drytausum_aero, tau3d_aero) + END IF + + 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, & + pdtphys, pplay, paprs, t_seri, rhcl, presnivs, & + mass_solu_aero, mass_solu_aero_pi, & + tau_aero, piz_aero, cg_aero, & + tausum_aero, tau3d_aero) + ! + !--natural aerosols + tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:) + piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:) + cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:) + !--all aerosols + tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:) + 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 ' & + // 'aerosols and iflag_rrtm=1' + CALL abort_physic(modname,abort_message,1) + ENDIF +#else + abort_message='You should compile with -rrtm if running ' & + // 'with iflag_rrtm=1' + CALL abort_physic(modname,abort_message,1) +#endif + ! + ENDIF + ELSE IF (iflag_rrtm .EQ.2) THEN ! ecrad RADIATION +#ifdef CPP_ECRAD + !--climatologies or INCA aerosols + CALL readaerosol_optic_ecrad( debut, aerosol_couple, ok_alw, ok_volcan, & + flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, & + pdtphys, pplay, paprs, t_seri, rhcl, presnivs, & + tr_seri, mass_solu_aero, mass_solu_aero_pi, m_allaer) +#else + abort_message='You should compile with -rad ecrad if running with iflag_rrtm=2' + CALL abort_physic(modname,abort_message,1) +#endif + ENDIF + + ELSE !--flag_aerosol = 0 + tausum_aero(:,:,:) = 0. + drytausum_aero(:,:) = 0. + mass_solu_aero(:,:) = 0. + mass_solu_aero_pi(:,:) = 0. + IF (iflag_rrtm .EQ. 0) THEN !--old radiation + tau_aero(:,:,:,:) = 1.e-15 + piz_aero(:,:,:,:) = 1. + cg_aero(:,:,:,:) = 0. + ELSE + tau_aero_sw_rrtm(:,:,:,:) = 1.e-15 + tau_aero_lw_rrtm(:,:,:,:) = 1.e-15 + piz_aero_sw_rrtm(:,:,:,:) = 1.0 + cg_aero_sw_rrtm(:,:,:,:) = 0.0 + ENDIF + ENDIF + ! + !--WMO criterion to determine tropopause + CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg) + ! + !--STRAT AEROSOL + !--updates tausum_aero,tau_aero,piz_aero,cg_aero + IF (flag_aerosol_strat.GT.0) THEN + IF (prt_level .GE.10) THEN + PRINT *,'appel a readaerosolstrat', mth_cur + ENDIF + IF (iflag_rrtm.EQ.0) THEN + IF (flag_aerosol_strat.EQ.1) THEN + CALL readaerosolstrato(debut) + ELSE + abort_message='flag_aerosol_strat must equal 1 for rrtm=0' + CALL abort_physic(modname,abort_message,1) + ENDIF + ELSE +#ifdef CPP_RRTM + IF (.NOT. CPPKEY_STRATAER) THEN + !--prescribed strat aerosols + !--only in the case of non-interactive strat aerosols + IF (flag_aerosol_strat.EQ.1) THEN + CALL readaerosolstrato1_rrtm(debut) + ELSEIF (flag_aerosol_strat.EQ.2) THEN + CALL readaerosolstrato2_rrtm(debut, ok_volcan) + ELSE + abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1' + CALL abort_physic(modname,abort_message,1) + ENDIF + END IF +#else + abort_message='You should compile with -rrtm if running ' & + // 'with iflag_rrtm=1' + CALL abort_physic(modname,abort_message,1) +#endif + ENDIF + ELSE + tausum_aero(:,:,id_STRAT_phy) = 0. + ENDIF + ! +#ifdef CPP_RRTM + IF (CPPKEY_STRATAER) THEN + !--compute stratospheric mask + CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg) + !--interactive strat aerosols + CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut) + END IF +#endif + !--fin STRAT AEROSOL + ! !lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek Index: LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90 (revision 6144) +++ LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90 (revision 6145) @@ -112,4 +112,5 @@ USE calwake_mod, ONLY : calwake, calwake_first USE lmdz_wake_ini, ONLY : wake_ini + USE lmdz_cv_ini, ONLY : cv_ini USE lmdz_cv_ini, ONLY : epmax, coef_epmax_cape, cvl_comp_threshold, cvl_sig2feed USE lmdz_cv_ini, ONLY : iflag_cvl_sigd, iflag_clw, ok_adj_ema @@ -5119,35 +5120,4 @@ !--------------------------------------------------------------------------- - 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) - picefra(i,k)=(radocond(i,k)*picefra(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_new_lscp=false and/or ok_icefra_lscp=false" - ! abort_message='inconsistency in cloud phase for blowing snow' - ! CALL abort_physic(modname,abort_message,1) - ! ENDIF - - ENDIF #ifdef ISO !#ifdef ISOVERIF @@ -5273,141 +5243,4 @@ ! - !------------------------------------------------------------------- - ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT - !------------------------------------------------------------------- - - ! 1. NUAGES CONVECTIFS - ! - !IM cf FH - ! IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke - IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke - snow_tiedtke=0. - ! print*,'avant calcul de la pseudo precip ' - ! print*,'iflag_cld_th',iflag_cld_th - IF (iflag_cld_th.eq.-1) THEN - rain_tiedtke=rain_con - ELSE - ! print*,'calcul de la pseudo precip ' - rain_tiedtke=0. - ! print*,'calcul de la pseudo precip 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 - ! - ! call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ') - ! - - ! Nuages diagnostiques pour Tiedtke - CALL diagcld1(paprs,pplay, & - !IM cf FH. rain_con,snow_con,ibas_con,itop_con, - 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. - - !IM/FH: 2011/02/23 - ! definition des points sur lesquels ls thermiques sont actifs - - 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 - ! ! Precipitation totale ! @@ -5480,4 +5313,5 @@ endif #endif + ! ! Calculer l'humidite relative pour diagnostique @@ -5513,243 +5347,174 @@ ENDDO - IF (ANY(type_trac == ['inca','inco'])) THEN ! ModThL -IF (CPPKEY_INCA) THEN - CALL VTe(VTphysiq) - CALL VTb(VTinca) - calday = REAL(days_elapsed + 1) + jH_cur - - CALL chemtime(itap+itau_phy-1, date0, phys_tstep, itap) - CALL AEROSOL_METEO_CALC( & - calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, & - prfl,psfl,pctsrf,cell_area, & - latitude_deg,longitude_deg,u10m,v10m) - - zxsnow_dummy(:) = 0.0 - ! INCA needs a cloud fraction that is not necessarily that - ! for radiation. Here we provide a cloud fraction calculated - ! the same manner as that in LMDZ5, LMDZ6 and LMDZ7 + + ! + !------------------------------------------------------------------- + ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT + !------------------------------------------------------------------- + + 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) + picefra(i,k)=(radocond(i,k)*picefra(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_new_lscp=false and/or ok_icefra_lscp=false" + ! abort_message='inconsistency in cloud phase for blowing snow' + ! CALL abort_physic(modname,abort_message,1) + ! ENDIF + + ENDIF + + ! 1. NUAGES CONVECTIFS + ! + !IM cf FH + ! IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke + IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke + snow_tiedtke=0. + ! print*,'avant calcul de la pseudo precip ' + ! print*,'iflag_cld_th',iflag_cld_th + IF (iflag_cld_th.eq.-1) THEN + rain_tiedtke=rain_con + ELSE + ! print*,'calcul de la pseudo precip ' + rain_tiedtke=0. + ! print*,'calcul de la pseudo precip 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 + ! + ! call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ') + ! + + ! Nuages diagnostiques pour Tiedtke + CALL diagcld1(paprs,pplay, & + !IM cf FH. rain_con,snow_con,ibas_con,itop_con, + 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 - cldfra_inca(i,k)=min(rneb(i,k)+rnebcon(i,k),1.) + 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. + + !IM/FH: 2011/02/23 + ! definition des points sur lesquels ls thermiques sont actifs + + 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 - - - CALL chemhook_begin (calday, & - days_elapsed+1, & - jH_cur, & - pctsrf(1,1), & - latitude_deg, & - longitude_deg, & - cell_area, & - paprs, & - pplay, & - coefh(1:klon,1:klev,is_ave), & - pphi, & - t_seri, & - u, & - v, & - rot, & - wo(:, :, 1), & - q_seri, & - zxtsol, & - zt2m, & - zxsnow_dummy, & - solsw, & - albsol1, & - rain_fall, & - snow_fall, & - itop_con, & - ibas_con, & - cldfra, & - nbp_lon, & - nbp_lat-1, & - tr_seri(:,:,1+nqCO2:nbtr), & - ftsol, & - paprs, & - cdragh, & - cdragm, & - pctsrf, & - pdtphys, & - itap) - - CALL VTe(VTinca) - CALL VTb(VTphysiq) -END IF - ENDIF !type_trac = inca or inco - IF (type_trac == 'repr') THEN -IF (CPPKEY_REPROBUS) THEN - !CALL chemtime_rep(itap+itau_phy-1, date0, dtime, itap) - CALL chemtime_rep(itap+itau_phy-1, date0, phys_tstep, itap) -END IF + ENDDO ENDIF - - ! - ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol. - ! - 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) - IF (flag_aerosol .GT. 0) THEN - IF (iflag_rrtm .EQ. 0) THEN !--old radiation - IF (.NOT. aerosol_couple) THEN - ! - CALL readaerosol_optic( & - debut, flag_aerosol, itap, jD_cur-jD_ref, & - pdtphys, pplay, paprs, t_seri, rhcl, presnivs, & - mass_solu_aero, mass_solu_aero_pi, & - tau_aero, piz_aero, cg_aero, & - tausum_aero, tau3d_aero) - ENDIF - ELSE IF (iflag_rrtm .EQ.1) THEN ! RRTM radiation - 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 - CALL splaerosol_optic_rrtm( ok_alw, pplay, paprs, t_seri, rhcl, & - tr_seri, mass_solu_aero, mass_solu_aero_pi, & - tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, & - tausum_aero, tau3d_aero) -ELSE - !--climatologies or INCA aerosols - CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, ok_volcan, & - flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, & - pdtphys, pplay, paprs, t_seri, rhcl, presnivs, & - tr_seri, mass_solu_aero, mass_solu_aero_pi, & - tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, & - tausum_aero, drytausum_aero, tau3d_aero) -END IF - - 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, & - pdtphys, pplay, paprs, t_seri, rhcl, presnivs, & - mass_solu_aero, mass_solu_aero_pi, & - tau_aero, piz_aero, cg_aero, & - tausum_aero, tau3d_aero) - ! - !--natural aerosols - tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:) - piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:) - cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:) - !--all aerosols - tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:) - 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 ' & - // 'aerosols and iflag_rrtm=1' - CALL abort_physic(modname,abort_message,1) - ENDIF -#else - abort_message='You should compile with -rrtm if running ' & - // 'with iflag_rrtm=1' - CALL abort_physic(modname,abort_message,1) -#endif - ! - ENDIF - ELSE IF (iflag_rrtm .EQ.2) THEN ! ecrad RADIATION -#ifdef CPP_ECRAD - !--climatologies or INCA aerosols - CALL readaerosol_optic_ecrad( debut, aerosol_couple, ok_alw, ok_volcan, & - flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, & - pdtphys, pplay, paprs, t_seri, rhcl, presnivs, & - tr_seri, mass_solu_aero, mass_solu_aero_pi, m_allaer) -#else - abort_message='You should compile with -rad ecrad if running with iflag_rrtm=2' - CALL abort_physic(modname,abort_message,1) -#endif - ENDIF - ELSE !--flag_aerosol = 0 - tausum_aero(:,:,:) = 0. - drytausum_aero(:,:) = 0. - mass_solu_aero(:,:) = 0. - mass_solu_aero_pi(:,:) = 0. - IF (iflag_rrtm .EQ. 0) THEN !--old radiation - tau_aero(:,:,:,:) = 1.e-15 - piz_aero(:,:,:,:) = 1. - cg_aero(:,:,:,:) = 0. - ELSE - tau_aero_sw_rrtm(:,:,:,:) = 1.e-15 - tau_aero_lw_rrtm(:,:,:,:) = 1.e-15 - piz_aero_sw_rrtm(:,:,:,:) = 1.0 - cg_aero_sw_rrtm(:,:,:,:) = 0.0 - ENDIF - ENDIF - ! - !--WMO criterion to determine tropopause - CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg) - ! - !--STRAT AEROSOL - !--updates tausum_aero,tau_aero,piz_aero,cg_aero - IF (flag_aerosol_strat.GT.0) THEN - IF (prt_level .GE.10) THEN - PRINT *,'appel a readaerosolstrat', mth_cur - ENDIF - IF (iflag_rrtm.EQ.0) THEN - IF (flag_aerosol_strat.EQ.1) THEN - CALL readaerosolstrato(debut) - ELSE - abort_message='flag_aerosol_strat must equal 1 for rrtm=0' - CALL abort_physic(modname,abort_message,1) - ENDIF - ELSE -#ifdef CPP_RRTM -IF (.NOT. CPPKEY_STRATAER) THEN - !--prescribed strat aerosols - !--only in the case of non-interactive strat aerosols - IF (flag_aerosol_strat.EQ.1) THEN - CALL readaerosolstrato1_rrtm(debut) - ELSEIF (flag_aerosol_strat.EQ.2) THEN - CALL readaerosolstrato2_rrtm(debut, ok_volcan) - ELSE - abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1' - CALL abort_physic(modname,abort_message,1) - ENDIF -END IF -#else - abort_message='You should compile with -rrtm if running ' & - // 'with iflag_rrtm=1' - CALL abort_physic(modname,abort_message,1) -#endif - ENDIF - ELSE - tausum_aero(:,:,id_STRAT_phy) = 0. - ENDIF -! -#ifdef CPP_RRTM -IF (CPPKEY_STRATAER) THEN -#ifdef ISO - CALL abort_gcm("physiq_mod", "StratAer isn't ISO-compatible for now, 07/24",1) -#endif - !--compute stratospheric mask - CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg) - !--interactive strat aerosols - CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut) -END IF -#endif - !--fin STRAT AEROSOL - ! ! Calculer les parametres optiques des nuages et quelques @@ -5833,4 +5598,245 @@ ! ENDIF + + + IF (ANY(type_trac == ['inca','inco'])) THEN ! ModThL +IF (CPPKEY_INCA) THEN + CALL VTe(VTphysiq) + CALL VTb(VTinca) + calday = REAL(days_elapsed + 1) + jH_cur + + CALL chemtime(itap+itau_phy-1, date0, phys_tstep, itap) + CALL AEROSOL_METEO_CALC( & + calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, & + prfl,psfl,pctsrf,cell_area, & + latitude_deg,longitude_deg,u10m,v10m) + + zxsnow_dummy(:) = 0.0 + ! INCA needs a cloud fraction that is not necessarily that + ! for radiation. Here we provide a cloud fraction calculated + ! the same manner as that in LMDZ5, LMDZ6 and LMDZ7 + DO k=1,klev + DO i=1,klon + cldfra_inca(i,k)=min(rneb(i,k)+rnebcon(i,k),1.) + ENDDO + ENDDO + + + CALL chemhook_begin (calday, & + days_elapsed+1, & + jH_cur, & + pctsrf(1,1), & + latitude_deg, & + longitude_deg, & + cell_area, & + paprs, & + pplay, & + coefh(1:klon,1:klev,is_ave), & + pphi, & + t_seri, & + u, & + v, & + rot, & + wo(:, :, 1), & + q_seri, & + zxtsol, & + zt2m, & + zxsnow_dummy, & + solsw, & + albsol1, & + rain_fall, & + snow_fall, & + itop_con, & + ibas_con, & + cldfra, & + nbp_lon, & + nbp_lat-1, & + tr_seri(:,:,1+nqCO2:nbtr), & + ftsol, & + paprs, & + cdragh, & + cdragm, & + pctsrf, & + pdtphys, & + itap) + + CALL VTe(VTinca) + CALL VTb(VTphysiq) +END IF + ENDIF !type_trac = inca or inco + IF (type_trac == 'repr') THEN +IF (CPPKEY_REPROBUS) THEN + !CALL chemtime_rep(itap+itau_phy-1, date0, dtime, itap) + CALL chemtime_rep(itap+itau_phy-1, date0, phys_tstep, itap) +END IF + ENDIF + + ! + ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol. + ! + 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) + IF (flag_aerosol .GT. 0) THEN + IF (iflag_rrtm .EQ. 0) THEN !--old radiation + IF (.NOT. aerosol_couple) THEN + ! + CALL readaerosol_optic( & + debut, flag_aerosol, itap, jD_cur-jD_ref, & + pdtphys, pplay, paprs, t_seri, rhcl, presnivs, & + mass_solu_aero, mass_solu_aero_pi, & + tau_aero, piz_aero, cg_aero, & + tausum_aero, tau3d_aero) + ENDIF + ELSE IF (iflag_rrtm .EQ.1) THEN ! RRTM radiation + 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 + CALL splaerosol_optic_rrtm( ok_alw, pplay, paprs, t_seri, rhcl, & + tr_seri, mass_solu_aero, mass_solu_aero_pi, & + tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, & + tausum_aero, tau3d_aero) +ELSE + !--climatologies or INCA aerosols + CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, ok_volcan, & + flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, & + pdtphys, pplay, paprs, t_seri, rhcl, presnivs, & + tr_seri, mass_solu_aero, mass_solu_aero_pi, & + tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, & + tausum_aero, drytausum_aero, tau3d_aero) +END IF + + 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, & + pdtphys, pplay, paprs, t_seri, rhcl, presnivs, & + mass_solu_aero, mass_solu_aero_pi, & + tau_aero, piz_aero, cg_aero, & + tausum_aero, tau3d_aero) + ! + !--natural aerosols + tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:) + piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:) + cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:) + !--all aerosols + tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:) + 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 ' & + // 'aerosols and iflag_rrtm=1' + CALL abort_physic(modname,abort_message,1) + ENDIF +#else + abort_message='You should compile with -rrtm if running ' & + // 'with iflag_rrtm=1' + CALL abort_physic(modname,abort_message,1) +#endif + ! + ENDIF + ELSE IF (iflag_rrtm .EQ.2) THEN ! ecrad RADIATION +#ifdef CPP_ECRAD + !--climatologies or INCA aerosols + CALL readaerosol_optic_ecrad( debut, aerosol_couple, ok_alw, ok_volcan, & + flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, & + pdtphys, pplay, paprs, t_seri, rhcl, presnivs, & + tr_seri, mass_solu_aero, mass_solu_aero_pi, m_allaer) +#else + abort_message='You should compile with -rad ecrad if running with iflag_rrtm=2' + CALL abort_physic(modname,abort_message,1) +#endif + ENDIF + ELSE !--flag_aerosol = 0 + tausum_aero(:,:,:) = 0. + drytausum_aero(:,:) = 0. + mass_solu_aero(:,:) = 0. + mass_solu_aero_pi(:,:) = 0. + IF (iflag_rrtm .EQ. 0) THEN !--old radiation + tau_aero(:,:,:,:) = 1.e-15 + piz_aero(:,:,:,:) = 1. + cg_aero(:,:,:,:) = 0. + ELSE + tau_aero_sw_rrtm(:,:,:,:) = 1.e-15 + tau_aero_lw_rrtm(:,:,:,:) = 1.e-15 + piz_aero_sw_rrtm(:,:,:,:) = 1.0 + cg_aero_sw_rrtm(:,:,:,:) = 0.0 + ENDIF + ENDIF + ! + !--WMO criterion to determine tropopause + CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg) + ! + !--STRAT AEROSOL + !--updates tausum_aero,tau_aero,piz_aero,cg_aero + IF (flag_aerosol_strat.GT.0) THEN + IF (prt_level .GE.10) THEN + PRINT *,'appel a readaerosolstrat', mth_cur + ENDIF + IF (iflag_rrtm.EQ.0) THEN + IF (flag_aerosol_strat.EQ.1) THEN + CALL readaerosolstrato(debut) + ELSE + abort_message='flag_aerosol_strat must equal 1 for rrtm=0' + CALL abort_physic(modname,abort_message,1) + ENDIF + ELSE +#ifdef CPP_RRTM +IF (.NOT. CPPKEY_STRATAER) THEN + !--prescribed strat aerosols + !--only in the case of non-interactive strat aerosols + IF (flag_aerosol_strat.EQ.1) THEN + CALL readaerosolstrato1_rrtm(debut) + ELSEIF (flag_aerosol_strat.EQ.2) THEN + CALL readaerosolstrato2_rrtm(debut, ok_volcan) + ELSE + abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1' + CALL abort_physic(modname,abort_message,1) + ENDIF +END IF +#else + abort_message='You should compile with -rrtm if running ' & + // 'with iflag_rrtm=1' + CALL abort_physic(modname,abort_message,1) +#endif + ENDIF + ELSE + tausum_aero(:,:,id_STRAT_phy) = 0. + ENDIF +! +#ifdef CPP_RRTM +IF (CPPKEY_STRATAER) THEN +#ifdef ISO + CALL abort_gcm("physiq_mod", "StratAer isn't ISO-compatible for now, 07/24",1) +#endif + !--compute stratospheric mask + CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg) + !--interactive strat aerosols + CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut) +END IF +#endif + !--fin STRAT AEROSOL + ! !lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek