Index: LMDZ5/trunk/libf/phylmd/physiq_mod.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/physiq_mod.F90 (revision 2617) +++ LMDZ5/trunk/libf/phylmd/physiq_mod.F90 (revision 2618) @@ -3382,119 +3382,230 @@ #endif END IF !type_trac = inca - ! - ! 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) - END IF - - if (ok_newmicro) then - IF (iflag_rrtm.NE.0) THEN + + + ! + ! 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, new_aod, 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 ! 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 (ok_cdnc.AND.NRADLP.NE.3) THEN + IF (NSW.EQ.6) THEN + !--new aerosol properties + ! + CALL readaerosol_optic_rrtm( debut, aerosol_couple, & + new_aod, flag_aerosol, 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, tau3d_aero) + + ELSE IF (NSW.EQ.2) THEN + !--for now we use the old aerosol properties + ! + CALL readaerosol_optic( & + debut, new_aod, 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,:) + ELSE + abort_message='Only NSW=2 or 6 are possible with ' & + // 'aerosols and iflag_rrtm=1' + call abort_physic(modname,abort_message,1) + ENDIF + + !--call LW optical properties for tropospheric aerosols + !--only works for INCA aerosol (aerosol_couple = TRUE) + CALL aeropt_lw_rrtm(aerosol_couple,paprs,tr_seri) + ! +#else + abort_message='You should compile with -rrtm if running ' & + // 'with iflag_rrtm=1' + call abort_physic(modname,abort_message,1) +#endif + ! + ENDIF + ENDIF + ELSE + tausum_aero(:,:,:) = 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 + ! + !--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 (flag_aerosol_strat.EQ.1) THEN + CALL readaerosolstrato1_rrtm(debut) + ELSEIF (flag_aerosol_strat.EQ.2) THEN + CALL stratosphere_mask(t_seri, pplay, latitude_deg) + CALL readaerosolstrato2_rrtm(debut) + ELSE + abort_message='flag_aerosol_strat must equal 1 or 2 for 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 + ENDIF + !--fin STRAT AEROSOL + ! + + ! 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) + END IF + + IF (ok_newmicro) then + IF (iflag_rrtm.NE.0) THEN +#ifdef CPP_RRTM + IF (ok_cdnc.AND.NRADLP.NE.3) THEN abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 ' & // 'pour ok_cdnc' - call abort_physic(modname,abort_message,1) - endif + 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) + abort_message='You should compile with -rrtm if running with '//'iflag_rrtm=1' + CALL abort_physic(modname,abort_message,1) #endif + ENDIF + CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, & + paprs, pplay, t_seri, cldliq, cldfra, & + cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, & + flwp, fiwp, flwc, fiwc, & + mass_solu_aero, mass_solu_aero_pi, & + cldtaupi, re, fl, ref_liq, ref_ice, & + ref_liq_pi, ref_ice_pi) + ELSE + CALL nuage (paprs, pplay, & + t_seri, cldliq, cldfra, cldtau, cldemi, & + cldh, cldl, cldm, cldt, cldq, & + ok_aie, & + mass_solu_aero, mass_solu_aero_pi, & + bl95_b0, bl95_b1, & + cldtaupi, re, fl) ENDIF - CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, & - paprs, pplay, t_seri, cldliq, cldfra, & - cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, & - flwp, fiwp, flwc, fiwc, & - mass_solu_aero, mass_solu_aero_pi, & - cldtaupi, re, fl, ref_liq, ref_ice, & - ref_liq_pi, ref_ice_pi) - else - CALL nuage (paprs, pplay, & - t_seri, cldliq, cldfra, cldtau, cldemi, & - cldh, cldl, cldm, cldt, cldq, & - ok_aie, & - mass_solu_aero, mass_solu_aero_pi, & - bl95_b0, bl95_b1, & - cldtaupi, re, fl) - endif - ! - !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 + ! + !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 + ELSE beta(i,k) = beta_free - endif - if (mskocean_beta) THEN + ENDIF + IF (mskocean_beta) THEN beta(i,k) = beta(i,k) * pctsrf(i,is_oce) - endif + 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 + ! + 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 - ENDDO - ! - endif - ! - ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol. - ! - IF (MOD(itaprad,radpas).EQ.0) THEN - - !albedo SB >>> - if(ok_chlorophyll)then + ! + ENDIF + + !lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek + IF (ok_chlorophyll) THEN print*,"-- reading chlorophyll" - call readchlorophyll(debut) - endif - !do i=1,klon - !if(chl_con(i)>1.) print*,i,chl_con(i),pctsrf(i,is_ter) - !enddo - !albedo SB <<< + CALL readchlorophyll(debut) + ENDIF ! @@ -3659,5 +3770,4 @@ cldtaupirad, & topswai_aero, solswai_aero) - #endif ELSE @@ -3674,5 +3784,4 @@ print *,' ->radlwsw, number 1 ' ENDIF - ! CALL radlwsw & @@ -3814,4 +3923,5 @@ ! radsol=solsw*swradcorr+sollw + if (ok_4xCO2atm) then radsolp=solswp*swradcorr+sollwp