| 1 | SUBROUTINE MACv2SP(pphis,pplay,paprs,xlon,xlat,tau_allaer,piz_allaer,cg_allaer) |
|---|
| 2 | ! |
|---|
| 3 | !--routine to read the MACv2SP plume and compute optical properties |
|---|
| 4 | !--requires flag_aerosol = 7 |
|---|
| 5 | !--feeds into aerosol optical properties and newmicro cloud droplet size if ok_cdnc activated |
|---|
| 6 | !--for this one needs to feed natural (pre-industrial) aerosols twice for nat and 1980 files |
|---|
| 7 | !--pre-ind aerosols (index=1) are not changed, present-day aerosols (index=2) are incremented |
|---|
| 8 | !--uses model year so year_cur needs to be correct in the model simulation |
|---|
| 9 | ! |
|---|
| 10 | !--aod_prof = AOD per layer |
|---|
| 11 | !--ssa_prof = SSA |
|---|
| 12 | !--asy_prof = asymetry parameter |
|---|
| 13 | !--dNovrN = enhancement factor for CDNC |
|---|
| 14 | ! |
|---|
| 15 | USE mo_simple_plumes, ONLY: sp_aop_profile |
|---|
| 16 | USE phys_cal_mod, ONLY : year_cur, day_cur, year_len |
|---|
| 17 | USE dimphy |
|---|
| 18 | USE aero_mod |
|---|
| 19 | USE phys_local_var_mod, ONLY: t_seri, od443aer, od550aer, od865aer, ec550aer, dryod550aer, od550lt1aer, dNovrN |
|---|
| 20 | !!USE YOMCST, ONLY : RD, RG |
|---|
| 21 | ! |
|---|
| 22 | USE yomcst_mod_h, ONLY: RPI, RCLUM, RHPLA, RKBOL, RNAVO & |
|---|
| 23 | , RDAY, REA, REPSM, RSIYEA, RSIDAY, ROMEGA & |
|---|
| 24 | , R_ecc, R_peri, R_incl & |
|---|
| 25 | , RA, RG, R1SA & |
|---|
| 26 | , RSIGMA & |
|---|
| 27 | , R, RMD, RMV, RD, RV, RCPD & |
|---|
| 28 | , RMO3, RMCO2, RMC, RMCH4, RMN2O, RMCFC11, RMCFC12 & |
|---|
| 29 | , RCPV, RCVD, RCVV, RKAPPA, RETV, eps_w & |
|---|
| 30 | , RCW, RCS & |
|---|
| 31 | , RLVTT, RLSTT, RLMLT, RTT, RATM & |
|---|
| 32 | , RESTT, RALPW, RBETW, RGAMW, RALPS, RBETS, RGAMS & |
|---|
| 33 | , RALPD, RBETD, RGAMD |
|---|
| 34 | IMPLICIT NONE |
|---|
| 35 | ! |
|---|
| 36 | |
|---|
| 37 | ! |
|---|
| 38 | REAL,DIMENSION(klon),INTENT(IN) :: pphis ! Geopotentiel de surface |
|---|
| 39 | REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) |
|---|
| 40 | REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour les interfaces de chaque couche (en Pa) |
|---|
| 41 | REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point |
|---|
| 42 | REAL,DIMENSION(klon),INTENT(IN) :: xlon ! longitudes pour chaque point |
|---|
| 43 | ! |
|---|
| 44 | REAL, DIMENSION(klon,klev,2,nbands_sw_rrtm), INTENT(OUT) :: tau_allaer ! epaisseur optique aerosol |
|---|
| 45 | REAL, DIMENSION(klon,klev,2,nbands_sw_rrtm), INTENT(OUT) :: piz_allaer ! single scattering albedo aerosol |
|---|
| 46 | REAL, DIMENSION(klon,klev,2,nbands_sw_rrtm), INTENT(OUT) :: cg_allaer ! asymmetry parameter aerosol |
|---|
| 47 | ! |
|---|
| 48 | REAL,DIMENSION(klon,klev) :: aod_prof, ssa_prof, asy_prof |
|---|
| 49 | REAL,DIMENSION(klon,klev) :: z, dz |
|---|
| 50 | REAL,DIMENSION(klon) :: oro, zrho, zt |
|---|
| 51 | ! |
|---|
| 52 | INTEGER, PARAMETER :: nmon = 12 |
|---|
| 53 | ! |
|---|
| 54 | REAL, PARAMETER :: l443 = 443.0, l550 = 550.0, l865 = 865.0 !--wavelengths in nm |
|---|
| 55 | ! |
|---|
| 56 | INTEGER, PARAMETER :: Nwvmax=25 |
|---|
| 57 | REAL, DIMENSION(0:Nwvmax), PARAMETER :: lambda=(/ 240.0, & !--this one is for band 1 |
|---|
| 58 | 280.0, 300.0, 330.0, 360.0, 400.0, & !--these are bounds of Streamer bands |
|---|
| 59 | 440.0, 480.0, 520.0, 570.0, 640.0, & |
|---|
| 60 | 690.0, 750.0, 780.0, 870.0, 1000.0, & |
|---|
| 61 | 1100.0, 1190.0, 1280.0, 1530.0, 1640.0, & |
|---|
| 62 | 2130.0, 2380.0, 2910.0, 3420.0, 4000.0 /) |
|---|
| 63 | ! |
|---|
| 64 | REAL, DIMENSION(1:Nwvmax-1), PARAMETER :: weight =(/ & !--and the weights to be given to the bands |
|---|
| 65 | 0.01, 4.05, 9.51, 15.99, 26.07, 33.10, & !--corresponding to a typical solar spectrum |
|---|
| 66 | 33.07, 39.91, 52.67, 27.89, 43.60, 13.67, & |
|---|
| 67 | 42.22, 40.12, 32.70, 14.44, 19.48, 14.23, & |
|---|
| 68 | 13.43, 16.42, 8.33, 0.95, 0.65, 2.76 /) |
|---|
| 69 | ! |
|---|
| 70 | REAL :: zlambda, zweight |
|---|
| 71 | REAL :: year_fr |
|---|
| 72 | ! |
|---|
| 73 | INTEGER band, i, k, Nwv |
|---|
| 74 | ! |
|---|
| 75 | ! define the height and dheight arrays |
|---|
| 76 | ! |
|---|
| 77 | oro(:) = pphis(:)/RG ! surface height in m |
|---|
| 78 | ! |
|---|
| 79 | DO k = 1, klev |
|---|
| 80 | zrho(:) = pplay(:,k)/t_seri(:,k)/RD ! air density in kg/m3 |
|---|
| 81 | dz(:,k) = (paprs(:,k)-paprs(:,k+1))/zrho(:)/RG ! layer thickness in m |
|---|
| 82 | IF (k==1) THEN |
|---|
| 83 | z(:,1) = oro(:) + (paprs(:,1)-pplay(:,1))/zrho(:)/RG ! altitude middle of first layer in m |
|---|
| 84 | zt(:) = oro(:) + dz(:,1) ! altitude top of first layer in m |
|---|
| 85 | ELSE |
|---|
| 86 | z(:,k) = zt(:) + (paprs(:,k)-pplay(:,k))/zrho(:)/RG ! altitude middle of layer k in m |
|---|
| 87 | zt(:) = zt(:) + dz(:,k) ! altitude top of layer k in m |
|---|
| 88 | ENDIF |
|---|
| 89 | ENDDO |
|---|
| 90 | ! |
|---|
| 91 | !--fractional year |
|---|
| 92 | ! |
|---|
| 93 | year_fr = FLOAT(year_cur) + (FLOAT(day_cur)-0.5) / FLOAT(year_len) |
|---|
| 94 | IF (year_fr.LT.1850.0.OR.year_fr.GE.2017.0) THEN |
|---|
| 95 | CALL abort_physic ('macv2sp','year not supported by plume model',1) |
|---|
| 96 | ENDIF |
|---|
| 97 | ! |
|---|
| 98 | !--call to sp routine -- 443 nm |
|---|
| 99 | ! |
|---|
| 100 | CALL sp_aop_profile ( & |
|---|
| 101 | klev ,klon ,l443 ,oro ,xlon ,xlat , & |
|---|
| 102 | year_fr ,z ,dz ,dNovrN ,aod_prof ,ssa_prof , & |
|---|
| 103 | asy_prof ) |
|---|
| 104 | ! |
|---|
| 105 | !--AOD calculations for diagnostics |
|---|
| 106 | od443aer(:)= od443aer(:)+SUM(aod_prof(:,:),dim=2) |
|---|
| 107 | ! |
|---|
| 108 | !--call to sp routine -- 550 nm |
|---|
| 109 | ! |
|---|
| 110 | CALL sp_aop_profile ( & |
|---|
| 111 | klev ,klon ,l550 ,oro ,xlon ,xlat , & |
|---|
| 112 | year_fr ,z ,dz ,dNovrN ,aod_prof ,ssa_prof , & |
|---|
| 113 | asy_prof ) |
|---|
| 114 | ! |
|---|
| 115 | !--AOD calculations for diagnostics |
|---|
| 116 | od550aer(:)=od550aer(:)+SUM(aod_prof(:,:),dim=2) |
|---|
| 117 | ! |
|---|
| 118 | !--dry AOD calculation for diagnostics |
|---|
| 119 | dryod550aer(:)=dryod550aer(:)+od550aer(:) |
|---|
| 120 | ! |
|---|
| 121 | !--fine-mode AOD calculation for diagnostics |
|---|
| 122 | od550lt1aer(:)=od550lt1aer(:)+od550aer(:) |
|---|
| 123 | ! |
|---|
| 124 | !--extinction coefficient for diagnostic |
|---|
| 125 | ec550aer(:,:)=ec550aer(:,:)+aod_prof(:,:)/dz(:,:) |
|---|
| 126 | ! |
|---|
| 127 | !--call to sp routine -- 865 nm |
|---|
| 128 | ! |
|---|
| 129 | CALL sp_aop_profile ( & |
|---|
| 130 | klev ,klon ,l865 ,oro ,xlon ,xlat , & |
|---|
| 131 | year_fr ,z ,dz ,dNovrN ,aod_prof ,ssa_prof , & |
|---|
| 132 | asy_prof ) |
|---|
| 133 | ! |
|---|
| 134 | !--AOD calculations for diagnostics |
|---|
| 135 | od865aer(:)=od865aer(:)+SUM(aod_prof(:,:),dim=2) |
|---|
| 136 | ! |
|---|
| 137 | !--re-weighting of piz and cg arrays before adding the anthropogenic aerosols |
|---|
| 138 | !--index 2 = all natural + anthropogenic aerosols |
|---|
| 139 | piz_allaer(:,:,2,:)=piz_allaer(:,:,2,:)*tau_allaer(:,:,2,:) |
|---|
| 140 | cg_allaer(:,:,2,:) =cg_allaer(:,:,2,:)*piz_allaer(:,:,2,:) |
|---|
| 141 | ! |
|---|
| 142 | !--now computing the same at many wavelengths to fill the model bands |
|---|
| 143 | ! |
|---|
| 144 | DO Nwv=0,Nwvmax-1 |
|---|
| 145 | |
|---|
| 146 | IF (Nwv.EQ.0) THEN !--RRTM spectral band 1 |
|---|
| 147 | zlambda=lambda(Nwv) |
|---|
| 148 | zweight=1.0 |
|---|
| 149 | band=1 |
|---|
| 150 | ELSEIF (Nwv.LE.5) THEN !--RRTM spectral band 2 |
|---|
| 151 | zlambda=0.5*(lambda(Nwv)+lambda(Nwv+1)) |
|---|
| 152 | zweight=weight(Nwv)/SUM(weight(1:5)) |
|---|
| 153 | band=2 |
|---|
| 154 | ELSEIF (Nwv.LE.10) THEN !--RRTM spectral band 3 |
|---|
| 155 | zlambda=0.5*(lambda(Nwv)+lambda(Nwv+1)) |
|---|
| 156 | zweight=weight(Nwv)/SUM(weight(6:10)) |
|---|
| 157 | band=3 |
|---|
| 158 | ELSEIF (Nwv.LE.16) THEN !--RRTM spectral band 4 |
|---|
| 159 | zlambda=0.5*(lambda(Nwv)+lambda(Nwv+1)) |
|---|
| 160 | zweight=weight(Nwv)/SUM(weight(11:16)) |
|---|
| 161 | band=4 |
|---|
| 162 | ELSEIF (Nwv.LE.21) THEN !--RRTM spectral band 5 |
|---|
| 163 | zlambda=0.5*(lambda(Nwv)+lambda(Nwv+1)) |
|---|
| 164 | zweight=weight(Nwv)/SUM(weight(17:21)) |
|---|
| 165 | band=5 |
|---|
| 166 | ELSE !--RRTM spectral band 6 |
|---|
| 167 | zlambda=0.5*(lambda(Nwv)+lambda(Nwv+1)) |
|---|
| 168 | zweight=weight(Nwv)/SUM(weight(22:Nwvmax-1)) |
|---|
| 169 | band=6 |
|---|
| 170 | ENDIF |
|---|
| 171 | ! |
|---|
| 172 | CALL sp_aop_profile ( & |
|---|
| 173 | klev ,klon ,zlambda ,oro ,xlon ,xlat , & |
|---|
| 174 | year_fr ,z ,dz ,dNovrN ,aod_prof ,ssa_prof , & |
|---|
| 175 | asy_prof ) |
|---|
| 176 | ! |
|---|
| 177 | !--adding up the quantities tau, piz*tau and cg*piz*tau |
|---|
| 178 | tau_allaer(:,:,2,band)=tau_allaer(:,:,2,band)+zweight*MAX(aod_prof(:,:),1.e-15) |
|---|
| 179 | piz_allaer(:,:,2,band)=piz_allaer(:,:,2,band)+zweight*MAX(aod_prof(:,:),1.e-15)*ssa_prof(:,:) |
|---|
| 180 | cg_allaer(:,:,2,band) =cg_allaer(:,:,2,band) +zweight*MAX(aod_prof(:,:),1.e-15)*ssa_prof(:,:)*asy_prof(:,:) |
|---|
| 181 | ! |
|---|
| 182 | ENDDO |
|---|
| 183 | ! |
|---|
| 184 | !--renpomalizing cg and piz now that MACv2SP increments have been added |
|---|
| 185 | cg_allaer(:,:,2,:) =cg_allaer(:,:,2,:) /piz_allaer(:,:,2,:) |
|---|
| 186 | piz_allaer(:,:,2,:)=piz_allaer(:,:,2,:)/tau_allaer(:,:,2,:) |
|---|
| 187 | ! |
|---|
| 188 | END SUBROUTINE MACv2SP |
|---|
