GCC Code Coverage Report


Directory: ./
File: phys/macv2sp.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 51 0.0%
Branches: 0 90 0.0%

Line Branch Exec Source
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 IMPLICIT NONE
23 !
24 include "YOMCST.h"
25 !
26 REAL,DIMENSION(klon),INTENT(IN) :: pphis ! Geopotentiel de surface
27 REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
28 REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour les interfaces de chaque couche (en Pa)
29 REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point
30 REAL,DIMENSION(klon),INTENT(IN) :: xlon ! longitudes pour chaque point
31 !
32 REAL, DIMENSION(klon,klev,2,nbands_sw_rrtm), INTENT(OUT) :: tau_allaer ! epaisseur optique aerosol
33 REAL, DIMENSION(klon,klev,2,nbands_sw_rrtm), INTENT(OUT) :: piz_allaer ! single scattering albedo aerosol
34 REAL, DIMENSION(klon,klev,2,nbands_sw_rrtm), INTENT(OUT) :: cg_allaer ! asymmetry parameter aerosol
35 !
36 REAL,DIMENSION(klon,klev) :: aod_prof, ssa_prof, asy_prof
37 REAL,DIMENSION(klon,klev) :: z, dz
38 REAL,DIMENSION(klon) :: oro, zrho, zt
39 !
40 INTEGER, PARAMETER :: nmon = 12
41 !
42 REAL, PARAMETER :: l443 = 443.0, l550 = 550.0, l865 = 865.0 !--wavelengths in nm
43 !
44 INTEGER, PARAMETER :: Nwvmax=25
45 REAL, DIMENSION(0:Nwvmax), PARAMETER :: lambda=(/ 240.0, & !--this one is for band 1
46 280.0, 300.0, 330.0, 360.0, 400.0, & !--these are bounds of Streamer bands
47 440.0, 480.0, 520.0, 570.0, 640.0, &
48 690.0, 750.0, 780.0, 870.0, 1000.0, &
49 1100.0, 1190.0, 1280.0, 1530.0, 1640.0, &
50 2130.0, 2380.0, 2910.0, 3420.0, 4000.0 /)
51 !
52 REAL, DIMENSION(1:Nwvmax-1), PARAMETER :: weight =(/ & !--and the weights to be given to the bands
53 0.01, 4.05, 9.51, 15.99, 26.07, 33.10, & !--corresponding to a typical solar spectrum
54 33.07, 39.91, 52.67, 27.89, 43.60, 13.67, &
55 42.22, 40.12, 32.70, 14.44, 19.48, 14.23, &
56 13.43, 16.42, 8.33, 0.95, 0.65, 2.76 /)
57 !
58 REAL :: zlambda, zweight
59 REAL :: year_fr
60 !
61 INTEGER band, i, k, Nwv
62 !
63 ! define the height and dheight arrays
64 !
65 oro(:) = pphis(:)/RG ! surface height in m
66 !
67 DO k = 1, klev
68 zrho(:) = pplay(:,k)/t_seri(:,k)/RD ! air density in kg/m3
69 dz(:,k) = (paprs(:,k)-paprs(:,k+1))/zrho(:)/RG ! layer thickness in m
70 IF (k==1) THEN
71 z(:,1) = oro(:) + (paprs(:,1)-pplay(:,1))/zrho(:)/RG ! altitude middle of first layer in m
72 zt(:) = oro(:) + dz(:,1) ! altitude top of first layer in m
73 ELSE
74 z(:,k) = zt(:) + (paprs(:,k)-pplay(:,k))/zrho(:)/RG ! altitude middle of layer k in m
75 zt(:) = zt(:) + dz(:,k) ! altitude top of layer k in m
76 ENDIF
77 ENDDO
78 !
79 !--fractional year
80 !
81 year_fr = FLOAT(year_cur) + (FLOAT(day_cur)-0.5) / FLOAT(year_len)
82 IF (year_fr.LT.1850.0.OR.year_fr.GE.2017.0) THEN
83 CALL abort_physic ('macv2sp','year not supported by plume model',1)
84 ENDIF
85 !
86 !--call to sp routine -- 443 nm
87 !
88 CALL sp_aop_profile ( &
89 klev ,klon ,l443 ,oro ,xlon ,xlat , &
90 year_fr ,z ,dz ,dNovrN ,aod_prof ,ssa_prof , &
91 asy_prof )
92 !
93 !--AOD calculations for diagnostics
94 od443aer(:)= od443aer(:)+SUM(aod_prof(:,:),dim=2)
95 !
96 !--call to sp routine -- 550 nm
97 !
98 CALL sp_aop_profile ( &
99 klev ,klon ,l550 ,oro ,xlon ,xlat , &
100 year_fr ,z ,dz ,dNovrN ,aod_prof ,ssa_prof , &
101 asy_prof )
102 !
103 !--AOD calculations for diagnostics
104 od550aer(:)=od550aer(:)+SUM(aod_prof(:,:),dim=2)
105 !
106 !--dry AOD calculation for diagnostics
107 dryod550aer(:)=dryod550aer(:)+od550aer(:)
108 !
109 !--fine-mode AOD calculation for diagnostics
110 od550lt1aer(:)=od550lt1aer(:)+od550aer(:)
111 !
112 !--extinction coefficient for diagnostic
113 ec550aer(:,:)=ec550aer(:,:)+aod_prof(:,:)/dz(:,:)
114 !
115 !--call to sp routine -- 865 nm
116 !
117 CALL sp_aop_profile ( &
118 klev ,klon ,l865 ,oro ,xlon ,xlat , &
119 year_fr ,z ,dz ,dNovrN ,aod_prof ,ssa_prof , &
120 asy_prof )
121 !
122 !--AOD calculations for diagnostics
123 od865aer(:)=od865aer(:)+SUM(aod_prof(:,:),dim=2)
124 !
125 !--re-weighting of piz and cg arrays before adding the anthropogenic aerosols
126 !--index 2 = all natural + anthropogenic aerosols
127 piz_allaer(:,:,2,:)=piz_allaer(:,:,2,:)*tau_allaer(:,:,2,:)
128 cg_allaer(:,:,2,:) =cg_allaer(:,:,2,:)*piz_allaer(:,:,2,:)
129 !
130 !--now computing the same at many wavelengths to fill the model bands
131 !
132 DO Nwv=0,Nwvmax-1
133
134 IF (Nwv.EQ.0) THEN !--RRTM spectral band 1
135 zlambda=lambda(Nwv)
136 zweight=1.0
137 band=1
138 ELSEIF (Nwv.LE.5) THEN !--RRTM spectral band 2
139 zlambda=0.5*(lambda(Nwv)+lambda(Nwv+1))
140 zweight=weight(Nwv)/SUM(weight(1:5))
141 band=2
142 ELSEIF (Nwv.LE.10) THEN !--RRTM spectral band 3
143 zlambda=0.5*(lambda(Nwv)+lambda(Nwv+1))
144 zweight=weight(Nwv)/SUM(weight(6:10))
145 band=3
146 ELSEIF (Nwv.LE.16) THEN !--RRTM spectral band 4
147 zlambda=0.5*(lambda(Nwv)+lambda(Nwv+1))
148 zweight=weight(Nwv)/SUM(weight(11:16))
149 band=4
150 ELSEIF (Nwv.LE.21) THEN !--RRTM spectral band 5
151 zlambda=0.5*(lambda(Nwv)+lambda(Nwv+1))
152 zweight=weight(Nwv)/SUM(weight(17:21))
153 band=5
154 ELSE !--RRTM spectral band 6
155 zlambda=0.5*(lambda(Nwv)+lambda(Nwv+1))
156 zweight=weight(Nwv)/SUM(weight(22:Nwvmax-1))
157 band=6
158 ENDIF
159 !
160 CALL sp_aop_profile ( &
161 klev ,klon ,zlambda ,oro ,xlon ,xlat , &
162 year_fr ,z ,dz ,dNovrN ,aod_prof ,ssa_prof , &
163 asy_prof )
164 !
165 !--adding up the quantities tau, piz*tau and cg*piz*tau
166 tau_allaer(:,:,2,band)=tau_allaer(:,:,2,band)+zweight*MAX(aod_prof(:,:),1.e-15)
167 piz_allaer(:,:,2,band)=piz_allaer(:,:,2,band)+zweight*MAX(aod_prof(:,:),1.e-15)*ssa_prof(:,:)
168 cg_allaer(:,:,2,band) =cg_allaer(:,:,2,band) +zweight*MAX(aod_prof(:,:),1.e-15)*ssa_prof(:,:)*asy_prof(:,:)
169 !
170 ENDDO
171 !
172 !--renpomalizing cg and piz now that MACv2SP increments have been added
173 cg_allaer(:,:,2,:) =cg_allaer(:,:,2,:) /piz_allaer(:,:,2,:)
174 piz_allaer(:,:,2,:)=piz_allaer(:,:,2,:)/tau_allaer(:,:,2,:)
175 !
176 END SUBROUTINE MACv2SP
177