Changeset 1150
- Timestamp:
- Apr 17, 2009, 5:34:01 PM (16 years ago)
- Location:
- LMDZ4/branches/LMDZ4-dev/libf/phylmd
- Files:
-
- 2 added
- 7 edited
- 2 copied
- 2 moved
-
aeropt.F (modified) (1 diff)
-
aeropt_2bands.F90 (copied) (copied from LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt.F) (1 diff)
-
aeropt_5wv.F90 (copied) (copied from LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt.F) (1 diff)
-
aerosol_optic.F90 (added)
-
clesphys.h (modified) (1 diff)
-
conf_phys.F90 (modified) (12 diffs)
-
phys_output_write.h (modified) (1 diff)
-
phys_state_var_mod.F90 (modified) (7 diffs)
-
physiq.F (modified) (20 diffs)
-
phytrac.F (modified) (4 diffs)
-
radlwsw_aero.F90 (moved) (moved from LMDZ4/branches/LMDZ4-dev/libf/phylmd/radlwsw.F) (1 diff)
-
readaerosol.F90 (moved) (moved from LMDZ4/branches/LMDZ4-dev/libf/phylmd/readsulfate.F) (1 diff)
-
sw_aero.F90 (added)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt.F
r766 r1150 16 16 c Arguments: 17 17 c 18 REAL paprs(klon,klev+1)19 REAL pplay(klon,klev), t_seri(klon,klev)20 REAL msulfate(klon,klev) ! masse sulfate ug SO4/m3 [ug/m^3]21 REAL RHcl(klon,klev) ! humidite relative ciel clair22 REAL tau_ae(klon,klev,2) ! epaisseur optique aerosol23 REAL piz_ae(klon,klev,2) ! single scattering albedo aerosol24 REAL cg_ae(klon,klev,2) ! asymmetry parameter aerosol25 REAL ai(klon) ! POLDER aerosol index18 REAL, INTENT(in) :: paprs(klon,klev+1) 19 REAL, INTENT(in) :: pplay(klon,klev), t_seri(klon,klev) 20 REAL, INTENT(in) :: msulfate(klon,klev) ! masse sulfate ug SO4/m3 [ug/m^3] 21 REAL, INTENT(in) :: RHcl(klon,klev) ! humidite relative ciel clair 22 REAL, INTENT(out) :: tau_ae(klon,klev,2) ! epaisseur optique aerosol 23 REAL, INTENT(out) :: piz_ae(klon,klev,2) ! single scattering albedo aerosol 24 REAL, INTENT(out) :: cg_ae(klon,klev,2) ! asymmetry parameter aerosol 25 REAL, INTENT(out) :: ai(klon) ! POLDER aerosol index 26 26 c 27 27 c Local -
LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_2bands.F90
r1149 r1150 1 ! 2 ! $Header$ 3 ! 4 SUBROUTINE aeropt(pplay, paprs, t_seri, msulfate, RHcl, 5 . tau_ae, piz_ae, cg_ae, ai ) 6 c 7 USE dimphy 8 IMPLICIT none 9 c 10 c 11 c 12 cym#include "dimensions.h" 13 cym#include "dimphy.h" 14 #include "YOMCST.h" 15 c 16 c Arguments: 17 c 18 REAL paprs(klon,klev+1) 19 REAL pplay(klon,klev), t_seri(klon,klev) 20 REAL msulfate(klon,klev) ! masse sulfate ug SO4/m3 [ug/m^3] 21 REAL RHcl(klon,klev) ! humidite relative ciel clair 22 REAL tau_ae(klon,klev,2) ! epaisseur optique aerosol 23 REAL piz_ae(klon,klev,2) ! single scattering albedo aerosol 24 REAL cg_ae(klon,klev,2) ! asymmetry parameter aerosol 25 REAL ai(klon) ! POLDER aerosol index 26 c 27 c Local 28 c 29 INTEGER i, k, inu 30 INTEGER RH_num, nbre_RH 31 PARAMETER (nbre_RH=12) 32 REAL RH_tab(nbre_RH) 33 REAL RH_MAX, DELTA, rh 34 PARAMETER (RH_MAX=95.) 35 DATA RH_tab/0.,10.,20.,30.,40.,50.,60.,70.,80.,85.,90.,95./ 36 REAL zrho, zdz 37 REAL taue670(klon) ! epaisseur optique aerosol absorption 550 nm 38 REAL taue865(klon) ! epaisseur optique aerosol extinction 865 nm 39 REAL alpha_aer_sulfate(nbre_RH,5) !--unit m2/g SO4 40 REAL alphasulfate 41 c 42 c Proprietes optiques 43 c 44 REAL alpha_aer(nbre_RH,2) !--unit m2/g SO4 45 REAL cg_aer(nbre_RH,2) 46 DATA alpha_aer/.500130E+01, .500130E+01, .500130E+01, 47 . .500130E+01, .500130E+01, .616710E+01, 48 . .826850E+01, .107687E+02, .136976E+02, 49 . .162972E+02, .211690E+02, .354833E+02, 50 . .139460E+01, .139460E+01, .139460E+01, 51 . .139460E+01, .139460E+01, .173910E+01, 52 . .244380E+01, .332320E+01, .440120E+01, 53 . .539570E+01, .734580E+01, .136038E+02 / 54 DATA cg_aer/.619800E+00, .619800E+00, .619800E+00, 55 . .619800E+00, .619800E+00, .662700E+00, 56 . .682100E+00, .698500E+00, .712500E+00, 57 . .721800E+00, .734600E+00, .755800E+00, 58 . .545600E+00, .545600E+00, .545600E+00, 59 . .545600E+00, .545600E+00, .583700E+00, 60 . .607100E+00, .627700E+00, .645800E+00, 61 . .658400E+00, .676500E+00, .708500E+00 / 62 DATA alpha_aer_sulfate/ 63 . 4.910,4.910,4.910,4.910,6.547,7.373, 64 . 8.373,9.788,12.167,14.256,17.924,28.433, 65 . 1.453,1.453,1.453,1.453,2.003,2.321, 66 . 2.711,3.282,4.287,5.210,6.914,12.305, 67 . 4.308,4.308,4.308,4.308,5.753,6.521, 68 . 7.449,8.772,11.014,12.999,16.518,26.772, 69 . 3.265,3.265,3.265,3.265,4.388,5.016, 70 . 5.775,6.868,8.745,10.429,13.457,22.538, 71 . 2.116,2.116,2.116,2.116,2.882,3.330, 72 . 3.876,4.670,6.059,7.327,9.650,16.883/ 73 c 74 DO i=1, klon 75 taue670(i)=0.0 76 taue865(i)=0.0 77 ENDDO 78 c 79 DO k=1, klev 80 DO i=1, klon 81 if (t_seri(i,k).eq.0) write (*,*) 'aeropt T ',i,k,t_seri(i,k) 82 if (pplay(i,k).eq.0) write (*,*) 'aeropt p ',i,k,pplay(i,k) 1 SUBROUTINE AEROPT_2BANDS( & 2 pdel, m_allaer, delt, RHcl, & 3 tau_allaer, piz_allaer, & 4 cg_allaer, fractnat_allaer, & 5 flag_aerosol, pplay, t_seri) 6 7 USE dimphy 8 9 10 ! Yves Balkanski le 12 avril 2006 11 ! Celine Deandreis 12 ! Anne Cozic Avril 2009 13 ! a partir d'une sous-routine de Johannes Quaas pour les sulfates 14 ! 15 IMPLICIT NONE 16 17 INCLUDE "YOMCST.h" 18 INCLUDE "iniprint.h" 19 20 ! 21 ! Input arguments: 22 ! 23 REAL, INTENT(in) :: pdel(KLON,KLEV) 24 REAL, INTENT(in) :: delt 25 REAL, DIMENSION(klon,klev,8), INTENT(in) :: m_allaer 26 REAL, INTENT(in) :: RHcl(KLON,KLEV) ! humidite relative ciel clair 27 REAL, DIMENSION(klon,10), INTENT(in) :: fractnat_allaer 28 INTEGER, INTENT(in) :: flag_aerosol 29 REAL, INTENT(in) :: pplay(klon,klev) 30 REAL, INTENT(in) :: t_seri(klon,klev) 31 32 ! 33 ! Output arguments: 34 ! 35 REAL, DIMENSION(klon,klev,9,2), INTENT(out) :: tau_allaer ! epaisseur optique aerosol 36 REAL, DIMENSION(klon,klev,9,2), INTENT(out) :: piz_allaer ! single scattering albedo aerosol 37 REAL, DIMENSION(klon,klev,9,2), INTENT(out) :: cg_allaer ! asymmetry parameter aerosol 38 39 ! 40 ! Local 41 ! 42 REAL, DIMENSION(klon,klev,10,2) :: tau_ae 43 REAL, DIMENSION(klon,klev,10,2) :: piz_ae 44 REAL, DIMENSION(klon,klev,10,2) :: cg_ae 45 LOGICAL :: soluble 46 INTEGER :: i, k, inu, m, mrfspecies 47 INTEGER :: spsol, spinsol 48 INTEGER :: RH_num, nbre_RH, nbsol_compaer, nbinsol_compaer 49 50 PARAMETER (nbre_RH=12) 51 PARAMETER (nbsol_compaer=5) ! 1- Seasalt AS: 2- Sesalt CS; 3- BC soluble; 4- POM soluble; 5- SO4. 52 PARAMETER (nbinsol_compaer=3) ! 1- Dust; 2- BC insoluble; 3- POM insoluble 53 REAL:: RH_tab(nbre_RH) 54 REAL:: RH_MAX, DELTA, rh 55 REAL:: tau_ae2b_int(KLON,KLEV,2) ! Intermediate computation of epaisseur optique aerosol 56 REAL:: piz_ae2b_int(KLON,KLEV,2) ! Intermediate computation of Single scattering albedo 57 REAL:: cg_ae2b_int(KLON,KLEV,2) ! Intermediate computation of Assymetry parameter 58 PARAMETER (RH_MAX=95.) 59 DATA RH_tab/0.,10.,20.,30.,40.,50.,60.,70.,80.,85.,90.,95./ 60 REAL :: zrho 61 REAL :: fac 62 REAL :: zdp1(klon,klev) 63 REAL, PARAMETER :: gravit = 9.80616 ! m2/s 64 INTEGER, ALLOCATABLE, DIMENSION(:) :: aerosol_name 65 INTEGER, PARAMETER :: id_SSSSM = 1 66 INTEGER, PARAMETER :: id_CSSSM = 2 67 INTEGER, PARAMETER :: id_ASSSM = 3 68 INTEGER, PARAMETER :: id_ASBCM = 4 69 INTEGER, PARAMETER :: id_ASPOMM = 5 70 INTEGER, PARAMETER :: id_ASSO4M = 6 71 INTEGER, PARAMETER :: id_CSSO4M = 7 72 INTEGER, PARAMETER :: id_CIDUSTM = 8 73 INTEGER, PARAMETER :: id_AIBCM = 9 74 INTEGER, PARAMETER :: id_AIPOMM = 10 75 INTEGER :: nb_aer 76 REAL, DIMENSION(klon,klev,8) :: mass_temp 77 78 ! 79 ! Proprietes optiques 80 ! 81 REAL:: alpha_aers_2bands(nbre_RH,2,nbsol_compaer) !--unit m2/g SO4 82 REAL:: alpha_aeri_2bands(2,nbinsol_compaer) 83 REAL:: cg_aers_2bands(nbre_RH,2,nbsol_compaer) !--unit 84 REAL:: cg_aeri_2bands(2,nbinsol_compaer) 85 REAL:: piz_aers_2bands(nbre_RH,2,nbsol_compaer) !-- unit 86 REAL:: piz_aeri_2bands(2,nbinsol_compaer) !-- unit 87 88 89 spsol = 0 90 spinsol = 0 91 92 93 DATA alpha_aers_2bands/ & 94 ! seasalt soluble Coarse Soluble (CS) 95 0.5090,0.6554,0.7129,0.7767,0.8529,1.2728, & 96 1.3820,1.5792,1.9173,2.2002,2.7173,4.1487, & 97 0.5167,0.6613,0.7221,0.7868,0.8622,1.3027, & 98 1.4227,1.6317,1.9887,2.2883,2.8356,4.3453, & 99 ! seasalt soluble Accumulation Soluble (AS) 100 4.125, 4.674, 5.005, 5.434, 5.985, 10.006, & 101 11.175,13.376,17.264,20.540,26.604, 42.349,& 102 4.187, 3.939, 3.919, 3.937, 3.995, 5.078, & 103 5.511, 6.434, 8.317,10.152,14.024, 26.537, & 104 ! bc soluble 105 7.675,7.675,7.675,7.675,7.675,7.675, & 106 7.675,7.675,10.433,11.984,13.767,15.567,& 107 4.720,4.720,4.720,4.720,4.720,4.720, & 108 4.720,4.720,6.081,6.793,7.567,9.344, & 109 ! pom soluble 110 5.503,5.503,5.503,5.503,5.588,5.957, & 111 6.404,7.340,8.545,10.319,13.595,20.398, & 112 1.402,1.402,1.402,1.402,1.431,1.562, & 113 1.715,2.032,2.425,2.991,4.193,7.133, & 114 ! sulfate 115 4.681,5.062,5.460,5.798,6.224,6.733, & 116 7.556,8.613,10.687,12.265,16.32,21.692, & 117 1.107,1.239,1.381,1.490,1.635,1.8030, & 118 2.071,2.407,3.126,3.940,5.539,7.921/ 119 120 121 DATA alpha_aeri_2bands/ & 122 ! dust insoluble 123 0.7661,0.7123,& 124 ! bc insoluble 125 10.360,4.437, & 126 ! pom insoluble 127 3.741,0.606/ 128 129 130 DATA cg_aers_2bands/ & 131 ! seasalt Coarse soluble (CS) 132 0.727, 0.747, 0.755, 0.761, 0.770, 0.788, & 133 0.792, 0.799, 0.805, 0.809, 0.815, 0.826, & 134 0.717, 0.738, 0.745, 0.752, 0.761, 0.779, & 135 0.781, 0.786, 0.793, 0.797, 0.803, 0.813, & 136 ! Sesalt Accumulation Soluble (AS) 137 0.727, 0.741, 0.748, 0.754, 0.761, 0.782, & 138 0.787, 0.792, 0.797, 0.799, 0.801, 0.799, & 139 0.606, 0.645, 0.658, 0.669, 0.681, 0.726, & 140 0.734, 0.746, 0.761, 0.770, 0.782, 0.798, & 141 ! bc soluble 142 .612, .612, .612, .612, .612, .612, & 143 .612, .612, .702, .734, .760, .796, & 144 .433, .433, .433, .433, .433, .433, & 145 .433, .433, .534, .575, .613, .669, & 146 ! pom soluble 147 .663, .663, .663, .663, .666, .674, & 148 .685, .702, .718, .737, .757, .777, & 149 .544, .544, .544, .544, .547, .554, & 150 .565, .583, .604, .631, .661, .698, & 151 ! sulfate 152 .658, .669, .680, .688, .698, .707, & 153 .719, .733, .752, .760, .773, .786, & 154 .544, .555, .565, .573, .583, .593, & 155 .610, .628, .655, .666, .692, .719/ 156 157 158 DATA cg_aeri_2bands/ & 159 ! dust insoluble 160 .701, .670, & 161 ! bc insoluble 162 .471, .297, & 163 ! pom insoluble 164 .568, .365/ 165 166 DATA piz_aers_2bands/& 167 ! seasalt Coarse soluble (CS) 168 1.000,1.000,1.000,1.000,1.000,1.000, & 169 1.000,1.000,1.000,1.000,1.000,1.000, & 170 0.992,0.989,0.987,0.986,0.986,0.980, & 171 0.980,0.978,0.976,0.976,0.974,0.971, & 172 ! seasalt Accumulation Soluble (AS) 173 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, & 174 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, & 175 0.970, 0.975, 0.976, 0.977, 0.978, 0.982, & 176 0.982, 0.983, 0.984, 0.984, 0.985, 0.985, & 177 ! bc soluble 178 .445, .445, .445, .445, .445, .445, & 179 .445, .445, .461, .480, .505, .528, & 180 .362, .362, .362, .362, .362, .362, & 181 .362, .362, .381, .405, .437, .483, & 182 ! pom soluble 183 .972, .972, .972, .972, .972, .974, & 184 .976, .979, .982, .986, .989, .992, & 185 .924, .924, .924, .924, .925, .927, & 186 .932, .938, .945, .952, .961, .970, & 187 ! sulfate 188 1.000,1.000,1.000,1.000,1.000,1.000, & 189 1.000,1.000,1.000,1.000,1.000,1.000, & 190 .992, .988, .988, .987, .986, .985, & 191 .985, .985, .984, .984, .984, .984 / 192 193 194 DATA piz_aeri_2bands/ & 195 ! dust insoluble 196 .963, .987, & 197 ! bc insoluble 198 .395, .264, & 199 ! pom insoluble 200 .966, .859/ 201 202 203 DO k=1, klev 204 DO i=1, klon 205 IF (t_seri(i,k).EQ.0.) THEN 206 WRITE(lunout,*) 't_seri(i,k)=0 for i=',i,'k=',k 207 CALL abort_gcm('aeropt_2bands','t_seri=0',1) 208 END IF 209 IF (pplay(i,k).EQ.0.) THEN 210 WRITE(lunout,*) 'pplay(i,k)=0 for i=',i,'k=',k 211 CALL abort_gcm('aeropt_2bands','pplay=0',1) 212 END IF 213 83 214 zrho=pplay(i,k)/t_seri(i,k)/RD ! kg/m3 84 zdz=(paprs(i,k)-paprs(i,k+1))/zrho/RG ! m 85 rh=MIN(RHcl(i,k)*100.,RH_MAX) 86 RH_num = INT( rh/10. + 1.) 87 IF (rh.LT.0.) STOP 'aeropt: RH < 0 not possible' 88 IF (rh.gt.85.) RH_num=10 89 IF (rh.gt.90.) RH_num=11 90 DELTA=(rh-RH_tab(RH_num))/(RH_tab(RH_num+1)-RH_tab(RH_num)) 91 c 92 inu=1 93 tau_ae(i,k,inu)=alpha_aer(RH_num,inu) + 94 . DELTA*(alpha_aer(RH_num+1,inu)-alpha_aer(RH_num,inu)) 95 tau_ae(i,k,inu)=tau_ae(i,k,inu)*msulfate(i,k)*zdz*1.e-6 96 piz_ae(i,k,inu)=1.0 97 cg_ae(i,k,inu)=cg_aer(RH_num,inu) + 98 . DELTA*(cg_aer(RH_num+1,inu)-cg_aer(RH_num,inu)) 99 c 100 inu=2 101 tau_ae(i,k,inu)=alpha_aer(RH_num,inu) + 102 . DELTA*(alpha_aer(RH_num+1,inu)-alpha_aer(RH_num,inu)) 103 tau_ae(i,k,inu)=tau_ae(i,k,inu)*msulfate(i,k)*zdz*1.e-6 104 piz_ae(i,k,inu)=1.0 105 cg_ae(i,k,inu)=cg_aer(RH_num,inu) + 106 . DELTA*(cg_aer(RH_num+1,inu)-cg_aer(RH_num,inu)) 107 cjq 108 cjq for aerosol index 109 c 110 alphasulfate=alpha_aer_sulfate(RH_num,4) + 111 . DELTA*(alpha_aer_sulfate(RH_num+1,4)- 112 . alpha_aer_sulfate(RH_num,4)) !--m2/g 113 c 114 taue670(i)=taue670(i)+ 115 . alphasulfate*msulfate(i,k)*zdz*1.e-6 116 c 117 alphasulfate=alpha_aer_sulfate(RH_num,5) + 118 . DELTA*(alpha_aer_sulfate(RH_num+1,5)- 119 . alpha_aer_sulfate(RH_num,5)) !--m2/g 120 c 121 taue865(i)=taue865(i)+ 122 . alphasulfate*msulfate(i,k)*zdz*1.e-6 123 124 ENDDO 125 ENDDO 126 c 127 DO i=1, klon 128 ai(i)=(-log(MAX(taue670(i),0.0001)/ 129 . MAX(taue865(i),0.0001))/log(670./865.)) * 130 . taue865(i) 131 ENDDO 132 c 133 RETURN 134 END 215 mass_temp(i,k,:) = m_allaer(i,k,:) / zrho / 1.e+9 216 217 ENDDO 218 ENDDO 219 220 IF (flag_aerosol .EQ. 1) THEN 221 nb_aer = 1 222 ALLOCATE (aerosol_name(nb_aer)) 223 aerosol_name(1) = id_ASSO4M 224 ELSEIF (flag_aerosol .EQ. 2) THEN 225 nb_aer = 2 226 ALLOCATE (aerosol_name(nb_aer)) 227 aerosol_name(1) = id_ASBCM 228 aerosol_name(2) = id_AIBCM 229 ELSEIF (flag_aerosol .EQ. 3) THEN 230 nb_aer = 2 231 ALLOCATE (aerosol_name(nb_aer)) 232 aerosol_name(1) = id_ASPOMM 233 aerosol_name(2) = id_AIPOMM 234 ELSEIF (flag_aerosol .EQ. 4) THEN 235 nb_aer = 5 236 ALLOCATE (aerosol_name(nb_aer)) 237 aerosol_name(1) = id_ASSO4M 238 aerosol_name(2) = id_ASBCM 239 aerosol_name(3) = id_AIBCM 240 aerosol_name(4) = id_ASPOMM 241 aerosol_name(5) = id_AIPOMM 242 ELSEIF (flag_aerosol .EQ. 5) THEN 243 nb_aer = 4 244 ALLOCATE (aerosol_name(nb_aer)) 245 aerosol_name(1) = id_ASBCM 246 aerosol_name(2) = id_AIBCM 247 aerosol_name(3) = id_ASPOMM 248 aerosol_name(4) = id_AIPOMM 249 ELSEIF (flag_aerosol .EQ. 6) THEN 250 nb_aer = 3 251 ALLOCATE (aerosol_name(nb_aer)) 252 aerosol_name(1) = id_ASSO4M 253 aerosol_name(2) = id_ASPOMM 254 aerosol_name(3) = id_AIPOMM 255 ENDIF 256 257 258 ! 259 ! loop over modes, use of precalculated nmd and corresponding sigma 260 ! loop over wavelengths 261 ! for each mass species in mode 262 ! interpolate from Sext to retrieve Sext_at_gridpoint_per_species 263 ! compute optical_thickness_at_gridpoint_per_species 264 265 tau_ae(:,:,:,:)=0. 266 piz_ae(:,:,:,:)=0. 267 cg_ae(:,:,:,:)=0. 268 tau_allaer(:,:,:,:)=0. 269 piz_allaer(:,:,:,:)=0. 270 cg_allaer(:,:,:,:)=0. 271 272 ! 273 ! Calculations that need to be done since we are not in the subroutines INCA 274 ! 275 ! air mass auxiliary variable --> zdp1 [kg/(m^2 *s)] 276 zdp1(:,:)=pdel(:,:)/(gravit*delt) 277 278 279 DO m=1,nb_aer ! tau is only computed for each mass 280 281 fac=1.0 282 IF (aerosol_name(m).EQ.id_SSSSM) THEN ! for now 283 soluble=.TRUE. 284 spsol=1 285 ELSEIF (aerosol_name(m).EQ.id_CSSSM) THEN 286 soluble=.TRUE. 287 spsol=1 288 ELSEIF (aerosol_name(m).EQ.id_ASSSM) THEN 289 soluble=.TRUE. 290 spsol=2 291 ELSEIF (aerosol_name(m).EQ.id_ASBCM) THEN 292 soluble=.TRUE. 293 spsol=3 294 ELSEIF (aerosol_name(m).EQ.id_ASPOMM) THEN 295 soluble=.TRUE. 296 spsol=4 297 ELSEIF ((aerosol_name(m).EQ.id_ASSO4M) .OR. (aerosol_name(m).EQ.id_CSSO4M)) THEN 298 soluble=.TRUE. 299 spsol=5 300 fac=1.375 ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD 301 ELSEIF (aerosol_name(m).EQ.id_CIDUSTM) THEN 302 soluble=.FALSE. 303 spinsol=1 304 ELSEIF (aerosol_name(m).EQ.id_AIBCM) THEN 305 soluble=.FALSE. 306 spinsol=2 307 ELSEIF (aerosol_name(m).EQ.id_AIPOMM) THEN 308 soluble=.FALSE. 309 spinsol=3 310 ELSE 311 CYCLE 312 ENDIF 313 314 315 tau_ae2b_int(:,:,:)=0. 316 piz_ae2b_int(:,:,:)=0. 317 cg_ae2b_int(:,:,:)=0. 318 319 DO k=1, KLEV 320 DO i=1, KLON 321 322 rh=MIN(RHcl(i,k)*100.,RH_MAX) 323 RH_num = INT( rh/10. + 1.) 324 325 IF (rh.GT.85.) RH_num=10 326 IF (rh.GT.90.) RH_num=11 327 DELTA=(rh-RH_tab(RH_num))/(RH_tab(RH_num+1)-RH_tab(RH_num)) 328 329 DO inu=1,2 330 IF (soluble) THEN 331 tau_ae2b_int(i,k,inu)= & 332 alpha_aers_2bands(RH_num,inu,spsol)+ & 333 DELTA* (alpha_aers_2bands(RH_num+1,inu,spsol) - & 334 alpha_aers_2bands(RH_num,inu,spsol)) 335 336 piz_ae2b_int(i,k,inu)= & 337 piz_aers_2bands(RH_num,inu,spsol) + & 338 DELTA* (piz_aers_2bands(RH_num+1,inu,spsol) - & 339 piz_aers_2bands(RH_num,inu,spsol)) 340 341 cg_ae2b_int(i,k,inu)= & 342 cg_aers_2bands(RH_num,inu,spsol) + & 343 DELTA* (cg_aers_2bands(RH_num+1,inu,spsol) - & 344 cg_aers_2bands(RH_num,inu,spsol)) 345 346 tau_ae(i,k,aerosol_name(m),inu) = mass_temp(i,k,spsol)*1000.*zdp1(i,k)*delt*tau_ae2b_int(i,k,inu)*fac 347 348 ELSE 349 tau_ae2b_int(i,k,inu) = alpha_aeri_2bands(inu,spinsol) 350 piz_ae2b_int(i,k,inu) = piz_aeri_2bands(inu,spinsol) 351 cg_ae2b_int(i,k,inu) = cg_aeri_2bands(inu,spinsol) 352 353 tau_ae(i,k,aerosol_name(m),inu) = mass_temp(i,k,5+ spinsol)*1000.*zdp1(i,k)*delt*tau_ae2b_int(i,k,inu)*fac 354 ENDIF 355 356 piz_ae(i,k,aerosol_name(m),inu) = piz_ae2b_int(i,k,inu) 357 358 cg_ae(i,k,aerosol_name(m),inu)= cg_ae2b_int(i,k,inu) 359 360 361 ENDDO ! boucle sur les bandes spectrale 362 ENDDO ! Boucle sur les points géographiques (grille horizontale) 363 ENDDO ! Boucle sur les niveaux verticaux 364 ENDDO ! Boucle sur les masses de traceurs 365 366 367 DO inu=1, 2 368 DO mrfspecies=1,9 369 DO k=1, KLEV 370 DO i=1, KLON 371 IF (mrfspecies .EQ. 2) THEN ! = total aerosol AER 372 tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASSO4M,inu)+tau_ae(i,k,id_CSSO4M,inu)+ & 373 tau_ae(i,k,id_ASBCM,inu)+tau_ae(i,k,id_AIBCM,inu)+ & 374 tau_ae(i,k,id_ASPOMM,inu)+tau_ae(i,k,id_AIPOMM,inu)+ & 375 tau_ae(i,k,id_ASSSM,inu)+tau_ae(i,k,id_CSSSM,inu)+tau_ae(i,k,id_SSSSM,inu)+ & 376 tau_ae(i,k,id_CIDUSTM,inu) 377 tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20) 378 379 piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASSO4M,inu)*piz_ae(i,k,id_ASSO4M,inu)+ & 380 tau_ae(i,k,id_CSSO4M,inu)*piz_ae(i,k,id_CSSO4M,inu)+ & 381 tau_ae(i,k,id_ASBCM,inu)*piz_ae(i,k,id_ASBCM,inu)+ & 382 tau_ae(i,k,id_AIBCM,inu)*piz_ae(i,k,id_AIBCM,inu)+ & 383 tau_ae(i,k,id_ASPOMM,inu)*piz_ae(i,k,id_ASPOMM,inu)+ & 384 tau_ae(i,k,id_AIPOMM,inu)*piz_ae(i,k,id_AIPOMM,inu)+ & 385 tau_ae(i,k,id_ASSSM,inu)*piz_ae(i,k,id_ASSSM,inu)+ & 386 tau_ae(i,k,id_CSSSM,inu)*piz_ae(i,k,id_CSSSM,inu)+ & 387 tau_ae(i,k,id_SSSSM,inu)*piz_ae(i,k,id_SSSSM,inu)+ & 388 tau_ae(i,k,id_CIDUSTM,inu)*piz_ae(i,k,id_CIDUSTM,inu))/tau_allaer(i,k,mrfspecies,inu) 389 piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20) 390 391 cg_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASSO4M,inu)*piz_ae(i,k,id_ASSO4M,inu)*cg_ae(i,k,id_ASSO4M,inu)+ & 392 tau_ae(i,k,id_CSSO4M,inu)*piz_ae(i,k,id_CSSO4M,inu)*cg_ae(i,k,id_CSSO4M,inu)+ & 393 tau_ae(i,k,id_ASBCM,inu)*piz_ae(i,k,id_ASBCM,inu)*cg_ae(i,k,id_ASBCM,inu)+ & 394 tau_ae(i,k,id_AIBCM,inu)*piz_ae(i,k,id_AIBCM,inu)*cg_ae(i,k,id_AIBCM,inu)+ & 395 tau_ae(i,k,id_ASPOMM,inu)*piz_ae(i,k,id_ASPOMM,inu)*cg_ae(i,k,id_ASPOMM,inu)+ & 396 tau_ae(i,k,id_AIPOMM,inu)*piz_ae(i,k,id_AIPOMM,inu)*cg_ae(i,k,id_AIPOMM,inu)+ & 397 tau_ae(i,k,id_ASSSM,inu)*piz_ae(i,k,id_ASSSM,inu)*cg_ae(i,k,id_ASSSM,inu)+ & 398 tau_ae(i,k,id_CSSSM,inu)*piz_ae(i,k,id_CSSSM,inu)*cg_ae(i,k,id_CSSSM,inu)+ & 399 tau_ae(i,k,id_SSSSM,inu)*piz_ae(i,k,id_SSSSM,inu)*cg_ae(i,k,id_SSSSM,inu)+ & 400 tau_ae(i,k,id_CIDUSTM,inu)*piz_ae(i,k,id_CIDUSTM,inu)*cg_ae(i,k,id_CIDUSTM,inu))/ & 401 (tau_allaer(i,k,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu)) 402 403 404 ELSEIF (mrfspecies .EQ. 3) THEN ! = natural aerosol NAT 405 tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASSO4M,inu)*fractnat_allaer(i,id_ASSO4M)+ & 406 tau_ae(i,k,id_CSSO4M,inu)*fractnat_allaer(i,id_CSSO4M)+ & 407 tau_ae(i,k,id_ASBCM,inu)*fractnat_allaer(i,id_ASBCM)+ & 408 tau_ae(i,k,id_AIBCM,inu)*fractnat_allaer(i,id_AIBCM)+ & 409 tau_ae(i,k,id_ASPOMM,inu)*fractnat_allaer(i,id_ASPOMM)+ & 410 tau_ae(i,k,id_AIPOMM,inu)*fractnat_allaer(i,id_AIPOMM)+ & 411 tau_ae(i,k,id_ASSSM,inu)*fractnat_allaer(i,id_ASSSM)+ & 412 tau_ae(i,k,id_CSSSM,inu)*fractnat_allaer(i,id_CSSSM)+ & 413 tau_ae(i,k,id_SSSSM,inu)*fractnat_allaer(i,id_SSSSM)+ & 414 tau_ae(i,k,id_CIDUSTM,inu)*fractnat_allaer(i,id_CIDUSTM) 415 tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20) 416 417 piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASSO4M,inu)*fractnat_allaer(i,id_ASSO4M)*piz_ae(i,k,id_ASSO4M,inu)+ & 418 tau_ae(i,k,id_CSSO4M,inu)*fractnat_allaer(i,id_CSSO4M)*piz_ae(i,k,id_CSSO4M,inu)+ & 419 tau_ae(i,k,id_ASBCM,inu)*fractnat_allaer(i,id_ASBCM)*piz_ae(i,k,id_ASBCM,inu)+ & 420 tau_ae(i,k,id_AIBCM,inu)*fractnat_allaer(i,id_AIBCM)*piz_ae(i,k,id_AIBCM,inu)+ & 421 tau_ae(i,k,id_ASPOMM,inu)*fractnat_allaer(i,id_ASPOMM)*piz_ae(i,k,id_ASPOMM,inu)+ & 422 tau_ae(i,k,id_AIPOMM,inu)*fractnat_allaer(i,id_AIPOMM)*piz_ae(i,k,id_AIPOMM,inu)+ & 423 tau_ae(i,k,id_ASSSM,inu)*fractnat_allaer(i,id_ASSSM)*piz_ae(i,k,id_ASSSM,inu)+ & 424 tau_ae(i,k,id_CSSSM,inu)*fractnat_allaer(i,id_CSSSM)*piz_ae(i,k,id_CSSSM,inu)+ & 425 tau_ae(i,k,id_SSSSM,inu)*fractnat_allaer(i,id_SSSSM)*piz_ae(i,k,id_SSSSM,inu)+ & 426 tau_ae(i,k,id_CIDUSTM,inu)*fractnat_allaer(i,id_CIDUSTM)*piz_ae(i,k,id_CIDUSTM,inu)) & 427 /tau_allaer(i,k,mrfspecies,inu) 428 piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20) 429 430 cg_allaer(i,k,mrfspecies,inu)=(& 431 tau_ae(i,k,id_ASSO4M,inu)*fractnat_allaer(i,id_ASSO4M)*piz_ae(i,k,id_ASSO4M,inu)*cg_ae(i,k,id_ASSO4M,inu)+ & 432 tau_ae(i,k,id_CSSO4M,inu)*fractnat_allaer(i,id_CSSO4M)*piz_ae(i,k,id_CSSO4M,inu)*cg_ae(i,k,id_CSSO4M,inu)+ & 433 tau_ae(i,k,id_ASBCM,inu)*fractnat_allaer(i,id_ASBCM)*piz_ae(i,k,id_ASBCM,inu)*cg_ae(i,k,id_ASBCM,inu)+ & 434 tau_ae(i,k,id_AIBCM,inu)*fractnat_allaer(i,id_AIBCM)*piz_ae(i,k,id_AIBCM,inu)*cg_ae(i,k,id_AIBCM,inu)+ & 435 tau_ae(i,k,id_ASPOMM,inu)*fractnat_allaer(i,id_ASPOMM)*piz_ae(i,k,id_ASPOMM,inu)*cg_ae(i,k,id_ASPOMM,inu)+ & 436 tau_ae(i,k,id_AIPOMM,inu)*fractnat_allaer(i,id_AIPOMM)*piz_ae(i,k,id_AIPOMM,inu)*cg_ae(i,k,id_AIPOMM,inu)+ & 437 tau_ae(i,k,id_ASSSM,inu)*fractnat_allaer(i,id_ASSSM)*piz_ae(i,k,id_ASSSM,inu)*cg_ae(i,k,id_ASSSM,inu)+ & 438 tau_ae(i,k,id_CSSSM,inu)*fractnat_allaer(i,id_CSSSM)*piz_ae(i,k,id_CSSSM,inu)*cg_ae(i,k,id_CSSSM,inu)+ & 439 tau_ae(i,k,id_SSSSM,inu)*fractnat_allaer(i,id_SSSSM)*piz_ae(i,k,id_SSSSM,inu)*cg_ae(i,k,id_SSSSM,inu)+ & 440 tau_ae(i,k,id_CIDUSTM,inu)*fractnat_allaer(i,id_CIDUSTM)*piz_ae(i,k,id_CIDUSTM,inu)*& 441 cg_ae(i,k,id_CIDUSTM,inu))/ & 442 (tau_allaer(i,k,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu)) 443 444 445 ELSEIF (mrfspecies .EQ. 4) THEN ! = BC 446 tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASBCM,inu)+tau_ae(i,k,id_AIBCM,inu) 447 tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20) 448 piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASBCM,inu)*piz_ae(i,k,id_ASBCM,inu) & 449 +tau_ae(i,k,id_AIBCM,inu)*piz_ae(i,k,id_AIBCM,inu))/ & 450 tau_allaer(i,k,mrfspecies,inu) 451 piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20) 452 cg_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASBCM,inu)*piz_ae(i,k,id_ASBCM,inu) *cg_ae(i,k,id_ASBCM,inu)& 453 +tau_ae(i,k,id_AIBCM,inu)*piz_ae(i,k,id_AIBCM,inu)*cg_ae(i,k,id_AIBCM,inu))/ & 454 (tau_allaer(i,k,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu)) 455 ELSEIF (mrfspecies .EQ. 5) THEN ! = SO4 456 tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASSO4M,inu)+tau_ae(i,k,id_CSSO4M,inu) 457 tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20) 458 piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_CSSO4M,inu)*piz_ae(i,k,id_CSSO4M,inu) & 459 +tau_ae(i,k,id_ASSO4M,inu)*piz_ae(i,k,id_ASSO4M,inu))/ & 460 tau_allaer(i,k,mrfspecies,inu) 461 piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20) 462 cg_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_CSSO4M,inu)*piz_ae(i,k,id_CSSO4M,inu) *cg_ae(i,k,id_CSSO4M,inu)& 463 +tau_ae(i,k,id_ASSO4M,inu)*piz_ae(i,k,id_ASSO4M,inu)*cg_ae(i,k,id_ASSO4M,inu))/ & 464 (tau_allaer(i,k,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu)) 465 466 ELSEIF (mrfspecies .EQ. 6) THEN ! = POM 467 tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASPOMM,inu)+tau_ae(i,k,id_AIPOMM,inu) 468 tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20) 469 piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASPOMM,inu)*piz_ae(i,k,id_ASPOMM,inu) & 470 +tau_ae(i,k,id_AIPOMM,inu)*piz_ae(i,k,id_AIPOMM,inu))/ & 471 tau_allaer(i,k,mrfspecies,inu) 472 piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20) 473 cg_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASPOMM,inu)*piz_ae(i,k,id_ASPOMM,inu) *cg_ae(i,k,id_ASPOMM,inu)& 474 +tau_ae(i,k,id_AIPOMM,inu)*piz_ae(i,k,id_AIPOMM,inu)*cg_ae(i,k,id_AIPOMM,inu))/ & 475 (tau_allaer(i,k,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu)) 476 ELSEIF (mrfspecies .EQ. 7) THEN ! = DUST 477 tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_CIDUSTM,inu) 478 tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20) 479 piz_allaer(i,k,mrfspecies,inu)=piz_ae(i,k,id_CIDUSTM,inu) 480 cg_allaer(i,k,mrfspecies,inu)=cg_ae(i,k,id_CIDUSTM,inu) 481 482 ELSEIF (mrfspecies .EQ. 8) THEN ! = SS 483 tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASSSM,inu)+tau_ae(i,k,id_CSSSM,inu)+tau_ae(i,k,id_SSSSM,inu) 484 tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20) 485 piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASSSM,inu)*piz_ae(i,k,id_ASSSM,inu) & 486 +tau_ae(i,k,id_CSSSM,inu)*piz_ae(i,k,id_CSSSM,inu) & 487 +tau_ae(i,k,id_SSSSM,inu)*piz_ae(i,k,id_SSSSM,inu))/ & 488 tau_allaer(i,k,mrfspecies,inu) 489 piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20) 490 cg_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASSSM,inu)*piz_ae(i,k,id_ASSSM,inu) *cg_ae(i,k,id_ASSSM,inu)& 491 +tau_ae(i,k,id_CSSSM,inu)*piz_ae(i,k,id_CSSSM,inu)*cg_ae(i,k,id_CSSSM,inu) & 492 +tau_ae(i,k,id_SSSSM,inu)*piz_ae(i,k,id_SSSSM,inu)*cg_ae(i,k,id_SSSSM,inu))/ & 493 (tau_allaer(i,k,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu)) 494 495 ELSEIF (mrfspecies .EQ. 9) THEN ! = NO3 496 tau_allaer(i,k,mrfspecies,inu)=0. ! preliminary 497 piz_allaer(i,k,mrfspecies,inu)=0. 498 cg_allaer(i,k,mrfspecies,inu)=0. 499 ENDIF 500 ENDDO 501 ENDDO 502 ENDDO 503 ENDDO 504 505 DEALLOCATE(aerosol_name) 506 507 END SUBROUTINE AEROPT_2BANDS -
LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_5wv.F90
r1149 r1150 1 ! 2 ! $Header$ 3 ! 4 SUBROUTINE aeropt(pplay, paprs, t_seri, msulfate, RHcl, 5 . tau_ae, piz_ae, cg_ae, ai ) 6 c 7 USE dimphy 8 IMPLICIT none 9 c 10 c 11 c 12 cym#include "dimensions.h" 13 cym#include "dimphy.h" 14 #include "YOMCST.h" 15 c 16 c Arguments: 17 c 18 REAL paprs(klon,klev+1) 19 REAL pplay(klon,klev), t_seri(klon,klev) 20 REAL msulfate(klon,klev) ! masse sulfate ug SO4/m3 [ug/m^3] 21 REAL RHcl(klon,klev) ! humidite relative ciel clair 22 REAL tau_ae(klon,klev,2) ! epaisseur optique aerosol 23 REAL piz_ae(klon,klev,2) ! single scattering albedo aerosol 24 REAL cg_ae(klon,klev,2) ! asymmetry parameter aerosol 25 REAL ai(klon) ! POLDER aerosol index 26 c 27 c Local 28 c 29 INTEGER i, k, inu 30 INTEGER RH_num, nbre_RH 31 PARAMETER (nbre_RH=12) 32 REAL RH_tab(nbre_RH) 33 REAL RH_MAX, DELTA, rh 34 PARAMETER (RH_MAX=95.) 35 DATA RH_tab/0.,10.,20.,30.,40.,50.,60.,70.,80.,85.,90.,95./ 36 REAL zrho, zdz 37 REAL taue670(klon) ! epaisseur optique aerosol absorption 550 nm 38 REAL taue865(klon) ! epaisseur optique aerosol extinction 865 nm 39 REAL alpha_aer_sulfate(nbre_RH,5) !--unit m2/g SO4 40 REAL alphasulfate 41 c 42 c Proprietes optiques 43 c 44 REAL alpha_aer(nbre_RH,2) !--unit m2/g SO4 45 REAL cg_aer(nbre_RH,2) 46 DATA alpha_aer/.500130E+01, .500130E+01, .500130E+01, 47 . .500130E+01, .500130E+01, .616710E+01, 48 . .826850E+01, .107687E+02, .136976E+02, 49 . .162972E+02, .211690E+02, .354833E+02, 50 . .139460E+01, .139460E+01, .139460E+01, 51 . .139460E+01, .139460E+01, .173910E+01, 52 . .244380E+01, .332320E+01, .440120E+01, 53 . .539570E+01, .734580E+01, .136038E+02 / 54 DATA cg_aer/.619800E+00, .619800E+00, .619800E+00, 55 . .619800E+00, .619800E+00, .662700E+00, 56 . .682100E+00, .698500E+00, .712500E+00, 57 . .721800E+00, .734600E+00, .755800E+00, 58 . .545600E+00, .545600E+00, .545600E+00, 59 . .545600E+00, .545600E+00, .583700E+00, 60 . .607100E+00, .627700E+00, .645800E+00, 61 . .658400E+00, .676500E+00, .708500E+00 / 62 DATA alpha_aer_sulfate/ 63 . 4.910,4.910,4.910,4.910,6.547,7.373, 64 . 8.373,9.788,12.167,14.256,17.924,28.433, 65 . 1.453,1.453,1.453,1.453,2.003,2.321, 66 . 2.711,3.282,4.287,5.210,6.914,12.305, 67 . 4.308,4.308,4.308,4.308,5.753,6.521, 68 . 7.449,8.772,11.014,12.999,16.518,26.772, 69 . 3.265,3.265,3.265,3.265,4.388,5.016, 70 . 5.775,6.868,8.745,10.429,13.457,22.538, 71 . 2.116,2.116,2.116,2.116,2.882,3.330, 72 . 3.876,4.670,6.059,7.327,9.650,16.883/ 73 c 74 DO i=1, klon 75 taue670(i)=0.0 76 taue865(i)=0.0 1 2 3 SUBROUTINE AEROPT_5WV(& 4 pdel, m_allaer, delt, & 5 RHcl, ai, flag_aerosol, & 6 pplay, t_seri) 7 8 USE DIMPHY 9 10 ! 11 ! Yves Balkanski le 12 avril 2006 12 ! Celine Deandreis 13 ! Anne Cozic Avril 2009 14 ! a partir d'une sous-routine de Johannes Quaas pour les sulfates 15 ! 16 ! 17 ! Refractive indices for seasalt come from Shettle and Fenn (1979) 18 ! 19 ! Refractive indices from water come from Hale and Querry (1973) 20 ! 21 ! Refractive indices from Ammonium Sulfate Toon and Pollack (1976) 22 ! 23 ! Refractive indices for Dust, internal mixture of minerals coated with 1.5% hematite 24 ! by Volume (Balkanski et al., 2006) 25 ! 26 ! Refractive indices for POM: Kinne (pers. Communication 27 ! 28 ! Refractive index for BC from Shettle and Fenn (1979) 29 ! 30 ! Shettle, E. P., & Fenn, R. W. (1979), Models for the aerosols of the lower atmosphere and 31 ! the effects of humidity variations on their optical properties, U.S. Air Force Geophysics 32 ! Laboratory Rept. AFGL-TR-79-0214, Hanscomb Air Force Base, MA. 33 ! 34 ! Hale, G. M. and M. R. Querry, Optical constants of water in the 200-nm to 200-m 35 ! wavelength region, Appl. Opt., 12, 555-563, 1973. 36 ! 37 ! Toon, O. B. and J. B. Pollack, The optical constants of several atmospheric aerosol species: 38 ! Ammonium sulfate, aluminum oxide, and sodium chloride, J. Geohys. Res., 81, 5733-5748, 39 ! 1976. 40 ! 41 ! Balkanski, Y., M. Schulz, T. Claquin And O. Boucher, Reevaluation of mineral aerosol 42 ! radiative forcings suggests a better agreement with satellite and AERONET data, Atmospheric 43 ! Chemistry and Physics Discussions., 6, pp 8383-8419, 2006. 44 ! 45 IMPLICIT NONE 46 INCLUDE "YOMCST.h" 47 ! 48 ! Input arguments: 49 ! 50 REAL, DIMENSION(klon,klev), INTENT(in) :: pdel 51 REAL, INTENT(in) :: delt 52 REAL, DIMENSION(klon,klev,8), INTENT(in) :: m_allaer 53 REAL, DIMENSION(klon,klev), INTENT(in) :: RHcl ! humidite relative ciel clair 54 INTEGER,INTENT(in) :: flag_aerosol 55 REAL, DIMENSION(klon,klev), INTENT(in) :: pplay 56 REAL, DIMENSION(klon,klev), INTENT(in) :: t_seri 57 58 ! 59 ! Output arguments: 60 ! 61 REAL, DIMENSION(klon), INTENT(out) :: ai ! POLDER aerosol index 62 63 ! 64 ! Local 65 ! 66 INTEGER, PARAMETER :: las = 5 67 LOGICAL :: soluble 68 69 INTEGER :: i, k, m 70 INTEGER :: spsol, spinsol, la 71 INTEGER :: RH_num 72 INTEGER, PARAMETER :: la443 = 1 73 INTEGER, PARAMETER :: la550 = 2 74 INTEGER, PARAMETER :: la670 = 3 75 INTEGER, PARAMETER :: la765 = 4 76 INTEGER, PARAMETER :: la865 = 5 77 INTEGER, PARAMETER :: nbre_RH=12 78 INTEGER, PARAMETER :: nbsol_compaer=5 ! 1- Seasalt AS: 2- Sesalt CS; 3- BC soluble; 4- POM soluble; 5- SO4. 79 INTEGER, PARAMETER :: nbinsol_compaer=3 ! 1- Dust; 2- BC insoluble; 3- POM insoluble 80 REAL :: zrho 81 REAL :: RH_tab(nbre_RH) 82 REAL :: DELTA, rh 83 REAL :: tau_ae5wv_int(KLON,KLEV,las) ! Intermediate computation of epaisseur optique aerosol 84 REAL :: piz_ae5wv_int(KLON,KLEV,las) ! Intermediate single scattering albedo aerosol 85 REAL :: cg_ae5wv_int(KLON,KLEV,las) ! Intermediate asymmetry parameter aerosol 86 REAL, PARAMETER :: RH_MAX=95. 87 DATA RH_tab/0.,10.,20.,30.,40.,50.,60.,70.,80.,85.,90.,95./ 88 REAL :: taue670(KLON) ! epaisseur optique aerosol absorption 550 nm 89 REAL :: taue865(KLON) ! epaisseur optique aerosol extinction 865 nm 90 REAL :: fac 91 REAL :: zdp1(klon,klev) 92 REAL, PARAMETER :: gravit = 9.80616 ! m2/s 93 INTEGER, ALLOCATABLE, DIMENSION(:) :: aerosol_name 94 INTEGER, PARAMETER :: id_SSSSM = 1 95 INTEGER, PARAMETER :: id_CSSSM = 2 96 INTEGER, PARAMETER :: id_ASSSM = 3 97 INTEGER, PARAMETER :: id_ASBCM = 4 98 INTEGER, PARAMETER :: id_ASPOMM = 5 99 INTEGER, PARAMETER :: id_ASSO4M = 6 100 INTEGER, PARAMETER :: id_CSSO4M = 7 101 INTEGER, PARAMETER :: id_CIDUSTM = 8 102 INTEGER, PARAMETER :: id_AIBCM = 9 103 INTEGER, PARAMETER :: id_AIPOMM = 10 104 INTEGER :: nb_aer 105 106 REAL :: tau3d(KLON,KLEV), piz3d(KLON,KLEV), cg3d(KLON,KLEV) 107 REAL :: abs3d(KLON,KLEV) ! epaisseur optique d'absorption 108 109 110 REAL :: alpha_aers_5wv(nbre_RH,las,nbsol_compaer) ! ext. coeff. Soluble comp. units *** m2/g 111 ! 1- Seasalt AS: 2- Sesalt CS; 3- BC; 4- POM; 5- SO4. 112 REAL :: alpha_aeri_5wv(las,nbinsol_compaer) ! ext. coeff. Insoluble comp. 1- Dust: 2- BC; 3- POM 113 REAL :: cg_aers_5wv(nbre_RH,las,nbsol_compaer) ! Asym. param. soluble comp. 114 ! 1- Seasalt AS: 2- Sesalt CS; 3- BC; 4- POM; 5- SO4. 115 REAL :: cg_aeri_5wv(las,nbinsol_compaer) ! Asym. param. insoluble comp. 1- Dust: 2- BC; 3- POM 116 REAL :: piz_aers_5wv(nbre_RH,las,nbsol_compaer) 117 ! 1- Seasalt AS: 2- Sesalt CS; 3- BC; 4- POM; 5- SO4. 118 REAL :: piz_aeri_5wv(las,nbinsol_compaer) ! Insoluble comp. 1- Dust: 2- BC; 3- POM 119 120 REAL, ALLOCATABLE :: TAUSUM(:,:,:) ! Aerosol optical thickness per species 121 REAL, ALLOCATABLE :: TAU(:,:,:,:) ! Aerosol optical thickness per species 122 REAL, DIMENSION(klon,klev,8) :: mass_temp 123 124 ! 125 ! Proprietes optiques 126 ! 127 REAL :: radry = 287.054 ! dry air mass constant 128 129 ! 130 ! 131 ! 132 ! From here on we look at the optical parameters at 5 wavelengths: 133 ! 443nm, 550, 670, 765 and 865 nm 134 ! le 12 AVRIL 2006 135 ! 136 DATA alpha_aers_5wv/ & 137 ! seasalt soluble CS 138 0.50,0.90,1.05,1.21,1.40,2.41, & 139 2.66,3.11,3.88,4.52,5.69,8.84, & 140 0.51,0.92,1.07,1.23,1.42,2.45, & 141 2.70,3.16,3.94,4.58,5.76,8.94, & 142 0.52,0.93,1.08,1.24,1.43,2.47, & 143 2.73,3.20,3.99,4.64,5.84,9.04, & 144 0.52,0.93,1.09,1.25,1.44,2.50, & 145 2.76,3.23,4.03,4.68,5.89,9.14, & 146 0.52,0.94,1.09,1.26,1.45,2.51, & 147 2.78,3.25,4.06,4.72,5.94,9.22, & 148 ! seasalt soluble AS 149 4.28, 7.17, 8.44, 9.85,11.60,22.44, & 150 25.34,30.54,39.38,46.52,59.33,91.77, & 151 3.40, 5.67, 6.69, 7.85, 9.32,19.03, & 152 21.78,26.88,35.87,43.40,57.33,93.43, & 153 2.48, 4.22, 5.02, 5.94, 7.11,15.29, & 154 17.70,22.31,30.73,38.06,52.15,90.59, & 155 1.90, 3.29, 3.94, 4.69, 5.65, 12.58, & 156 14.68,18.77,26.41,33.25,46.77,85.50, & 157 1.47, 2.59, 3.12, 3.74, 4.54, 10.42, & 158 12.24,15.82,22.66,28.91,41.54,79.33, & 159 ! bc soluble 160 7.930,7.930,7.930,7.930,7.930,7.930, & 161 7.930,7.930,10.893,12.618,14.550,16.613, & 162 7.658,7.658,7.658,7.658,7.658,7.658, & 163 7.658,7.658,10.351,11.879,13.642,15.510, & 164 7.195,7.195,7.195,7.195,7.195,7.195, & 165 7.195,7.195,9.551,10.847,12.381,13.994, & 166 6.736,6.736,6.736,6.736,6.736,6.736, & 167 6.736,6.736,8.818,9.938,11.283,12.687, & 168 6.277,6.277,6.277,6.277,6.277,6.277, & 169 6.277,6.277,8.123,9.094,10.275,11.501, & 170 ! pom soluble 171 6.676,6.676,6.676,6.676,6.710,6.934, & 172 7.141,7.569,8.034,8.529,9.456,10.511, & 173 5.109,5.109,5.109,5.109,5.189,5.535, & 174 5.960,6.852,8.008,9.712,12.897,19.676, & 175 3.718,3.718,3.718,3.718,3.779,4.042, & 176 4.364,5.052,5.956,7.314,9.896,15.688, & 177 2.849,2.849,2.849,2.849,2.897,3.107, & 178 3.365,3.916,4.649,5.760,7.900,12.863, & 179 2.229,2.229,2.229,2.229,2.268,2.437, & 180 2.645,3.095,3.692,4.608,6.391,10.633, & 181 ! Sulfate 182 5.751,6.215,6.690,7.024,7.599,8.195, & 183 9.156,10.355,12.660,14.823,18.908,24.508, & 184 4.320,4.675,5.052,5.375,5.787,6.274, & 185 7.066,8.083,10.088,12.003,15.697,21.133, & 186 3.079,3.351,3.639,3.886,4.205,4.584, & 187 5.206,6.019,7.648,9.234,12.391,17.220, & 188 2.336,2.552,2.781,2.979,3.236,3.540, & 189 4.046,4.711,6.056,7.388,10.093,14.313, & 190 1.777,1.949,2.134,2.292,2.503,2.751, & 191 3.166,3.712,4.828,5.949,8.264,11.922/ 192 193 DATA alpha_aeri_5wv/ & 194 ! dust insoluble 195 0.759, 0.770, 0.775, 0.775, 0.772, & 196 !!jb bc insoluble 197 11.536,10.033, 8.422, 7.234, 6.270, & 198 ! pom insoluble 199 5.042, 3.101, 1.890, 1.294, 0.934/ 200 201 DATA cg_aers_5wv/ & 202 ! seasalt soluble (CS) 203 0.730,0.753,0.760,0.766,0.772,0.793, & 204 0.797,0.802,0.809,0.813,0.820,0.830, & 205 0.719,0.744,0.751,0.757,0.764,0.786, & 206 0.791,0.796,0.803,0.808,0.815,0.826, & 207 0.721,0.744,0.750,0.756,0.762,0.784, & 208 0.787,0.793,0.800,0.804,0.811,0.822, & 209 0.717,0.741,0.747,0.753,0.759,0.780, & 210 0.784,0.789,0.795,0.800,0.806,0.817, & 211 0.715,0.739,0.745,0.751,0.757,0.777, & 212 0.781,0.786,0.793,0.797,0.803,0.814, & 213 ! seasalt soluble (AS) 214 0.698,0.722,0.729,0.736,0.743,0.765, & 215 0.768,0.773,0.777,0.779,0.781,0.779, & 216 0.682,0.710,0.719,0.727,0.735,0.764, & 217 0.769,0.776,0.783,0.787,0.791,0.792, & 218 0.658,0.691,0.701,0.710,0.720,0.756, & 219 0.763,0.771,0.782,0.788,0.795,0.801, & 220 0.632,0.668,0.679,0.690,0.701,0.743, & 221 0.750,0.762,0.775,0.783,0.792,0.804, & 222 0.605,0.644,0.656,0.669,0.681,0.729, & 223 0.737,0.750,0.765,0.775,0.787,0.803, & 224 ! bc soluble 225 .651, .651, .651, .651, .651, .651, & 226 .651, .651, .738, .764, .785, .800, & 227 .597, .597, .597, .597, .597, .597, & 228 .597, .597, .695, .725, .751, .770, & 229 .543, .543, .543, .543, .543, .543, & 230 .543, .543, .650, .684, .714, .736, & 231 .504, .504, .504, .504, .504, .504, & 232 .504, .504, .614, .651, .683, .708, & 233 .469, .469, .469, .469, .469, .469, & 234 .469, .469, .582, .620, .655, .681, & 235 ! pom soluble 236 .679, .679, .679, .679, .683, .691, & 237 .703, .720, .736, .751, .766, .784, & 238 .656, .656, .656, .656, .659, .669, & 239 .681, .699, .717, .735, .750, .779, & 240 .623, .623, .623, .623, .627, .637, & 241 .649, .668, .688, .709, .734, .762, & 242 .592, .592, .592, .592, .595, .605, & 243 .618, .639, .660, .682, .711, .743, & 244 .561, .561, .561, .561, .565, .575, & 245 .588, .609, .632, .656, .688, .724, & 246 ! sulfate 247 .671, .684, .697, .704, .714, .723, & 248 .734, .746, .762, .771, .781, .789, & 249 .653, .666, .678, .687, .697, .707, & 250 .719, .732, .751, .762, .775, .789, & 251 .622, .635, .648, .657, .667, .678, & 252 .691, .705, .728, .741, .758, .777, & 253 .591, .604, .617, .627, .638, .650, & 254 .664, .679, .704, .719, .739, .761, & 255 .560, .574, .587, .597, .609, .621, & 256 .637, .653, .680, .697, .719, .745/ 257 ! 258 259 DATA cg_aeri_5wv/& 260 ! dust insoluble 261 0.714, 0.697, 0.688, 0.683, 0.679, & 262 ! bc insoluble 263 0.511, 0.445, 0.384, 0.342, 0.307, & 264 !c pom insoluble 265 0.596, 0.536, 0.466, 0.409, 0.359/ 266 ! 267 DATA piz_aers_5wv/& 268 ! seasalt soluble (CS) 269 1.000,1.000,1.000,1.000,1.000,1.000, & 270 1.000,1.000,1.000,1.000,1.000,1.000, & 271 1.000,1.000,1.000,1.000,1.000,1.000, & 272 1.000,1.000,1.000,1.000,1.000,1.000, & 273 1.000,1.000,1.000,1.000,1.000,1.000, & 274 1.000,1.000,1.000,1.000,1.000,1.000, & 275 1.000,1.000,1.000,1.000,1.000,1.000, & 276 1.000,1.000,1.000,1.000,1.000,1.000, & 277 1.000,1.000,1.000,1.000,1.000,1.000, & 278 1.000,1.000,1.000,1.000,1.000,1.000, & 279 ! seasalt soluble (AS) 280 1.000,1.000,1.000,1.000,1.000,1.000, & 281 1.000,1.000,1.000,1.000,1.000,1.000, & 282 1.000,1.000,1.000,1.000,1.000,1.000, & 283 1.000,1.000,1.000,1.000,1.000,1.000, & 284 1.000,1.000,1.000,1.000,1.000,1.000, & 285 1.000,1.000,1.000,1.000,1.000,1.000, & 286 1.000,1.000,1.000,1.000,1.000,1.000, & 287 1.000,1.000,1.000,1.000,1.000,1.000, & 288 1.000,1.000,1.000,1.000,1.000,1.000, & 289 1.000,1.000,1.000,1.000,1.000,1.000, & 290 ! bc soluble 291 .445, .445, .445, .445, .445, .445, & 292 .445, .445, .470, .487, .508, .531, & 293 .442, .442, .442, .442, .442, .442, & 294 .442, .442, .462, .481, .506, .533, & 295 .427, .427, .427, .427, .427, .427, & 296 .427, .427, .449, .470, .497, .526, & 297 .413, .413, .413, .413, .413, .413, & 298 .413, .413, .437, .458, .486, .516, & 299 .399, .399, .399, .399, .399, .399, & 300 .399, .399, .423, .445, .473, .506, & 301 ! pom soluble 302 .975, .975, .975, .975, .975, .977, & 303 .979, .982, .984, .987, .990, .994, & 304 .972, .972, .972, .972, .973, .974, & 305 .977, .980, .983, .986, .989, .993, & 306 .963, .963, .963, .963, .964, .966, & 307 .969, .974, .977, .982, .986, .991, & 308 .955, .955, .955, .955, .955, .958, & 309 .962, .967, .972, .977, .983, .989, & 310 .944, .944, .944, .944, .944, .948, & 311 .952, .959, .962, .972, .979, .987, & 312 ! sulfate 313 1.000,1.000,1.000,1.000,1.000,1.000, & 314 1.000,1.000,1.000,1.000,1.000,1.000, & 315 1.000,1.000,1.000,1.000,1.000,1.000, & 316 1.000,1.000,1.000,1.000,1.000,1.000, & 317 1.000,1.000,1.000,1.000,1.000,1.000, & 318 1.000,1.000,1.000,1.000,1.000,1.000, & 319 1.000,1.000,1.000,1.000,1.000,1.000, & 320 1.000,1.000,1.000,1.000,1.000,1.000, & 321 1.000,1.000,1.000,1.000,1.000,1.000, & 322 1.000,1.000,1.000,1.000,1.000,1.000/ 323 ! 324 DATA piz_aeri_5wv/& 325 ! dust insoluble 326 0.944, 0.970, 0.977, 0.982, 0.987, & 327 ! bc insoluble 328 0.415, 0.387, 0.355, 0.328, 0.301, & 329 ! pom insoluble 330 0.972, 0.963, 0.943, 0.923, 0.897/ 331 332 333 334 DO k=1, klev 335 DO i=1, klon 336 IF (t_seri(i,k).EQ.0) stop 'stop aeropt_5wv T ' 337 IF (pplay(i,k).EQ.0) stop 'stop aeropt_5wv p ' 338 339 zrho=pplay(i,k)/t_seri(i,k)/RD ! kg/m3 340 mass_temp(i,k,:) = m_allaer(i,k,:) / zrho / 1.e+9 341 342 ENDDO 343 ENDDO 344 345 346 347 IF (flag_aerosol .EQ. 1) THEN 348 nb_aer = 1 349 ALLOCATE (aerosol_name(nb_aer)) 350 aerosol_name(1) = id_ASSO4M 351 ELSEIF (flag_aerosol .EQ. 2) THEN 352 nb_aer = 2 353 ALLOCATE (aerosol_name(nb_aer)) 354 aerosol_name(1) = id_ASBCM 355 aerosol_name(2) = id_AIBCM 356 ELSEIF (flag_aerosol .EQ. 3) THEN 357 nb_aer = 2 358 ALLOCATE (aerosol_name(nb_aer)) 359 aerosol_name(1) = id_ASPOMM 360 aerosol_name(2) = id_AIPOMM 361 ELSEIF (flag_aerosol .EQ. 4) THEN 362 nb_aer = 5 363 ALLOCATE (aerosol_name(nb_aer)) 364 aerosol_name(1) = id_ASSO4M 365 aerosol_name(2) = id_ASBCM 366 aerosol_name(3) = id_AIBCM 367 aerosol_name(4) = id_ASPOMM 368 aerosol_name(5) = id_AIPOMM 369 ELSEIF (flag_aerosol .EQ. 5) THEN 370 nb_aer = 4 371 ALLOCATE (aerosol_name(nb_aer)) 372 aerosol_name(1) = id_ASBCM 373 aerosol_name(2) = id_AIBCM 374 aerosol_name(3) = id_ASPOMM 375 aerosol_name(4) = id_AIPOMM 376 ELSEIF (flag_aerosol .EQ. 6) THEN 377 nb_aer = 3 378 ALLOCATE (aerosol_name(nb_aer)) 379 aerosol_name(1) = id_ASSO4M 380 aerosol_name(2) = id_ASPOMM 381 aerosol_name(3) = id_AIPOMM 382 ENDIF 383 384 ALLOCATE (tausum(klon,las,nb_aer)) 385 ALLOCATE (tau(klon,klev,las,nb_aer)) 386 387 388 389 390 ! 391 ! loop over modes, use of precalculated nmd and corresponding sigma 392 ! loop over wavelengths 393 ! for each mass species in mode 394 ! interpolate from Sext to retrieve Sext_at_gridpoint_per_species 395 ! compute optical_thickness_at_gridpoint_per_species 396 397 ai(:)=0. 398 399 tau_ae5wv_int(:,:,:)=0. 400 piz_ae5wv_int(:,:,:)=0. 401 cg_ae5wv_int(:,:,:)=0. 402 tausum(:,:,:)=0. 403 404 tau(:,:,:,:)=0. 405 ! 406 ! Calculations that need to be done since we are not in the subroutines INCA 407 ! 408 ! air mass auxiliary variable --> zdp1 [kg/(m^2 *s)] 409 zdp1=pdel/(gravit*delt) 410 411 DO m=1,nb_aer ! tau is only computed for each mass 412 413 fac=1.0 414 IF (aerosol_name(m).EQ.id_SSSSM) THEN ! for now 415 soluble=.TRUE. 416 spsol=1 417 ELSEIF (aerosol_name(m).EQ.id_CSSSM) THEN 418 soluble=.TRUE. 419 spsol=1 420 ELSEIF (aerosol_name(m).EQ.id_ASSSM) THEN 421 soluble=.TRUE. 422 spsol=2 423 ELSEIF (aerosol_name(m).EQ.id_ASBCM) THEN 424 soluble=.TRUE. 425 spsol=3 426 ELSEIF (aerosol_name(m).EQ.id_ASPOMM) THEN 427 soluble=.TRUE. 428 spsol=4 429 ELSEIF ((aerosol_name(m).EQ.id_ASSO4M) .OR. (aerosol_name(m).EQ.id_CSSO4M)) THEN 430 soluble=.TRUE. 431 spsol=5 432 fac=1.375 ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD 433 ELSEIF (aerosol_name(m).EQ.id_CIDUSTM) THEN 434 soluble=.FALSE. 435 spinsol=1 436 ELSEIF (aerosol_name(m).EQ.id_AIBCM) THEN 437 soluble=.FALSE. 438 spinsol=2 439 ELSEIF (aerosol_name(m).EQ.id_AIPOMM) THEN 440 soluble=.FALSE. 441 spinsol=3 442 ELSE 443 CYCLE 444 ENDIF 445 446 DO la=1,las 447 tau3d(:,:)=0. 448 piz3d(:,:)=0. 449 cg3d(:,:)=0. 450 abs3d(:,:)=0. 451 452 DO k=1, KLEV 453 DO i=1, KLON 454 455 rh=MIN(RHcl(i,k)*100.,RH_MAX) 456 RH_num = INT( rh/10. + 1.) 457 458 IF (rh.GT.85.) RH_num=10 459 IF (rh.GT.90.) RH_num=11 460 DELTA=(rh-RH_tab(RH_num))/(RH_tab(RH_num+1)-RH_tab(RH_num)) 461 462 IF (soluble) THEN 463 tau_ae5wv_int(i,k,la) = & 464 alpha_aers_5wv(RH_num,la,spsol)+DELTA* & 465 (alpha_aers_5wv(RH_num+1,la,spsol) - & 466 alpha_aers_5wv(RH_num,la,spsol)) 467 468 piz_ae5wv_int(i,k,la) = & 469 piz_aers_5wv(RH_num,la,spsol)+DELTA* & 470 (piz_aers_5wv(RH_num+1,la,spsol) - & 471 piz_aers_5wv(RH_num,la,spsol)) 472 473 cg_ae5wv_int(i,k,la) = & 474 cg_aers_5wv(RH_num,la,spsol)+DELTA* & 475 (cg_aers_5wv(RH_num+1,la,spsol) - & 476 cg_aers_5wv(RH_num,la,spsol)) 477 478 tau3d(i,k) = & 479 mass_temp(i,k,spsol)*1000.*zdp1(i,k)*tau_ae5wv_int(i,k,la)*delt*fac 480 481 ELSE 482 tau_ae5wv_int(i,k,la) = alpha_aeri_5wv(la,spinsol) 483 piz_ae5wv_int(i,k,la) = piz_aeri_5wv(la,spinsol) 484 cg_ae5wv_int(i,k,la) = cg_aeri_5wv(la,spinsol) 485 486 tau3d(i,k) = & 487 mass_temp(i,k,5+spinsol)*1000.*zdp1(i,k)*tau_ae5wv_int(i,k,la)*delt*fac 488 ENDIF 489 490 491 ENDDO ! Boucle sur les points géographiques (grille horizontale) 492 ENDDO ! Boucle sur les niveaux verticaux 493 494 tau(:,:,la,m)=tau3d(:,:) 495 496 DO k=1, KLEV 497 DO i=1,KLON 498 tausum(i,la,m)=tausum(i,la,m)+tau3d(i,k) 499 ENDDO 77 500 ENDDO 78 c 79 DO k=1, klev 80 DO i=1, klon 81 if (t_seri(i,k).eq.0) write (*,*) 'aeropt T ',i,k,t_seri(i,k) 82 if (pplay(i,k).eq.0) write (*,*) 'aeropt p ',i,k,pplay(i,k) 83 zrho=pplay(i,k)/t_seri(i,k)/RD ! kg/m3 84 zdz=(paprs(i,k)-paprs(i,k+1))/zrho/RG ! m 85 rh=MIN(RHcl(i,k)*100.,RH_MAX) 86 RH_num = INT( rh/10. + 1.) 87 IF (rh.LT.0.) STOP 'aeropt: RH < 0 not possible' 88 IF (rh.gt.85.) RH_num=10 89 IF (rh.gt.90.) RH_num=11 90 DELTA=(rh-RH_tab(RH_num))/(RH_tab(RH_num+1)-RH_tab(RH_num)) 91 c 92 inu=1 93 tau_ae(i,k,inu)=alpha_aer(RH_num,inu) + 94 . DELTA*(alpha_aer(RH_num+1,inu)-alpha_aer(RH_num,inu)) 95 tau_ae(i,k,inu)=tau_ae(i,k,inu)*msulfate(i,k)*zdz*1.e-6 96 piz_ae(i,k,inu)=1.0 97 cg_ae(i,k,inu)=cg_aer(RH_num,inu) + 98 . DELTA*(cg_aer(RH_num+1,inu)-cg_aer(RH_num,inu)) 99 c 100 inu=2 101 tau_ae(i,k,inu)=alpha_aer(RH_num,inu) + 102 . DELTA*(alpha_aer(RH_num+1,inu)-alpha_aer(RH_num,inu)) 103 tau_ae(i,k,inu)=tau_ae(i,k,inu)*msulfate(i,k)*zdz*1.e-6 104 piz_ae(i,k,inu)=1.0 105 cg_ae(i,k,inu)=cg_aer(RH_num,inu) + 106 . DELTA*(cg_aer(RH_num+1,inu)-cg_aer(RH_num,inu)) 107 cjq 108 cjq for aerosol index 109 c 110 alphasulfate=alpha_aer_sulfate(RH_num,4) + 111 . DELTA*(alpha_aer_sulfate(RH_num+1,4)- 112 . alpha_aer_sulfate(RH_num,4)) !--m2/g 113 c 114 taue670(i)=taue670(i)+ 115 . alphasulfate*msulfate(i,k)*zdz*1.e-6 116 c 117 alphasulfate=alpha_aer_sulfate(RH_num,5) + 118 . DELTA*(alpha_aer_sulfate(RH_num+1,5)- 119 . alpha_aer_sulfate(RH_num,5)) !--m2/g 120 c 121 taue865(i)=taue865(i)+ 122 . alphasulfate*msulfate(i,k)*zdz*1.e-6 123 124 ENDDO 125 ENDDO 126 c 127 DO i=1, klon 128 ai(i)=(-log(MAX(taue670(i),0.0001)/ 129 . MAX(taue865(i),0.0001))/log(670./865.)) * 130 . taue865(i) 131 ENDDO 132 c 133 RETURN 134 END 501 502 ENDDO ! boucle sur les longueurs d'onde 503 ENDDO ! Boucle sur les masses de traceurs 504 505 506 taue670(:) = SUM(tausum(:,la670,:),dim=2) 507 taue865(:) = SUM(tausum(:,la865,:),dim=2) 508 509 DO i=1, klon 510 ai(i)=-LOG(MAX(taue670(i),0.0001)/ & 511 MAX(taue865(i),0.0001))/LOG(670./865.) 512 ENDDO 513 514 515 END SUBROUTINE AEROPT_5WV -
LMDZ4/branches/LMDZ4-dev/libf/phylmd/clesphys.h
r1143 r1150 40 40 INTEGER lev_histhf, lev_histday, lev_histmth 41 41 CHARACTER*4 type_run 42 ! aer_type: pour utiliser un fichier constant dans read sulfate42 ! aer_type: pour utiliser un fichier constant dans readaerosol 43 43 CHARACTER*8 :: aer_type 44 44 LOGICAL ok_isccp, ok_regdyn -
LMDZ4/branches/LMDZ4-dev/libf/phylmd/conf_phys.F90
r1143 r1150 12 12 & iflag_cldcon, & 13 13 & iflag_ratqs,ratqsbas,ratqshaut, & 14 & ok_ade, ok_aie, aerosol_couple, & 14 & ok_ade, ok_aie, aerosol_couple, & 15 & flag_aerosol, new_aod, & 15 16 & bl95_b0, bl95_b1,& 16 17 & iflag_thermals,nsplit_thermals,tau_thermals, & … … 60 61 logical :: ok_LES 61 62 LOGICAL :: ok_ade, ok_aie, aerosol_couple 63 INTEGER :: flag_aerosol 64 LOGICAL :: new_aod 62 65 REAL :: bl95_b0, bl95_b1 63 66 real :: fact_cldcon, facttemps,ratqsbas,ratqshaut … … 71 74 logical,SAVE :: ok_LES_omp 72 75 LOGICAL,SAVE :: ok_ade_omp, ok_aie_omp, aerosol_couple_omp 76 INTEGER, SAVE :: flag_aerosol_omp 77 LOGICAL, SAVE :: new_aod_omp 73 78 REAL,SAVE :: bl95_b0_omp, bl95_b1_omp 74 79 REAL,SAVE :: freq_ISCCP_omp, ecrit_ISCCP_omp … … 239 244 CALL getin('aerosol_couple',aerosol_couple_omp) 240 245 246 ! 247 !Config Key = flag_aerosol 248 !Config Desc = which aerosol is use for coupled model 249 !Config Def = 1 250 !Config Help = Used in physiq.F 251 ! 252 ! - flag_aerosol=1 => so4 only (defaut) 253 ! - flag_aerosol=2 => bc only 254 ! - flag_aerosol=3 => pom only 255 ! - flag_aerosol=4 => all aerosol 256 ! - flag_aerosol=5 => bcpom 257 ! - flag_aerosol=6 => pomsulf 258 259 flag_aerosol_omp = 1 260 CALL getin('flag_aerosol',flag_aerosol_omp) 261 262 ! Temporary variable for testing purpose!! 263 !Config Key = new_aod 264 !Config Desc = which calcul of aeropt 265 !Config Def = false 266 !Config Help = Used in physiq.F 267 ! 268 new_aod_omp = .true. 269 CALL getin('new_aod',new_aod_omp) 270 241 271 ! 242 272 !Config Key = aer_type 243 273 !Config Desc = Use a constant field for the aerosols 244 274 !Config Def = scenario 245 !Config Help = Used in read sulfate.F275 !Config Help = Used in readaerosol.F90 246 276 ! 247 277 aer_type_omp = 'scenario' … … 918 948 !Config Help = 919 949 ! 920 iflag_coupl = 0950 iflag_coupl_omp = 0 921 951 call getin('iflag_coupl',iflag_coupl_omp) 922 952 … … 927 957 !Config Help = 928 958 ! 929 iflag_clos = 1959 iflag_clos_omp = 1 930 960 call getin('iflag_clos',iflag_clos_omp) 931 961 ! … … 935 965 !Config Help = 936 966 ! 937 iflag_cvl_sigd = 0967 iflag_cvl_sigd_omp = 0 938 968 call getin('iflag_cvl_sigd',iflag_cvl_sigd_omp) 939 969 … … 943 973 !Config Help = 944 974 ! 945 iflag_wake = 0975 iflag_wake_omp = 0 946 976 call getin('iflag_wake',iflag_wake_omp) 947 977 … … 1266 1296 ok_aie = ok_aie_omp 1267 1297 aerosol_couple = aerosol_couple_omp 1298 flag_aerosol=flag_aerosol_omp 1299 new_aod=new_aod_omp 1268 1300 aer_type = aer_type_omp 1269 1301 bl95_b0 = bl95_b0_omp … … 1326 1358 WRITE(numout,*)' ERROR version_ocean=',version_ocean,' not valid with slab ocean' 1327 1359 CALL abort_gcm('conf_phys','version_ocean not valid',1) 1360 END IF 1361 1362 ! Test sur new_aod. Ce flag permet de retrouver les resultats de l'AR4 1363 ! il n'est utilisable que lors du couplage avec le SO4 seul 1364 IF (ok_ade .OR. ok_aie) THEN 1365 IF ( .NOT. new_aod .AND. flag_aerosol .NE. 1) THEN 1366 CALL abort_gcm('conf_phys','new_aod=.FALSE. not compatible avec flag_aerosol=1',1) 1367 END IF 1328 1368 END IF 1329 1369 … … 1396 1436 write(numout,*)' ok_aie = ',ok_aie 1397 1437 write(numout,*)' aerosol_couple = ', aerosol_couple 1438 write(numout,*)' flag_aerosol = ', flag_aerosol 1439 write(numout,*)' new_aod = ', new_aod 1398 1440 write(numout,*)' aer_type = ',aer_type 1399 1441 write(numout,*)' bl95_b0 = ',bl95_b0 … … 1409 1451 write(numout,*)' type_run = ',type_run 1410 1452 write(numout,*)' ok_isccp = ',ok_isccp 1411 WRITE(numout,*)' solarlong0 = ', solarlong01453 write(numout,*)' solarlong0 = ', solarlong0 1412 1454 write(numout,*)' qsol0 = ', qsol0 1413 1455 write(numout,*)' inertie_sol = ', inertie_sol -
LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_output_write.h
r1123 r1150 832 832 833 833 IF (ok_ade) THEN 834 IF (o_topswad%flag(iff)<=lev_files(iff)) THEN 835 CALL histwrite_phy(nid_files(iff),o_topswad%name,itau_w,topswad) 836 ENDIF 837 IF (o_solswad%flag(iff)<=lev_files(iff)) THEN 838 CALL histwrite_phy(nid_files(iff),o_solswad%name,itau_w,solswad) 839 ENDIF 834 IF (o_topswad%flag(iff)<=lev_files(iff)) THEN 835 CALL histwrite_phy(nid_files(iff),o_topswad%name,itau_w, 836 $ topswad_aero) 837 ENDIF 838 IF (o_solswad%flag(iff)<=lev_files(iff)) THEN 839 CALL histwrite_phy(nid_files(iff),o_solswad%name,itau_w, 840 $ solswad_aero) 841 ENDIF 840 842 ENDIF 841 843 842 844 IF (ok_aie) THEN 843 IF (o_topswai%flag(iff)<=lev_files(iff)) THEN 844 CALL histwrite_phy(nid_files(iff),o_topswai%name,itau_w,topswai) 845 ENDIF 846 IF (o_solswai%flag(iff)<=lev_files(iff)) THEN 847 CALL histwrite_phy(nid_files(iff),o_solswai%name,itau_w,solswai) 848 ENDIF 845 IF (o_topswai%flag(iff)<=lev_files(iff)) THEN 846 CALL histwrite_phy(nid_files(iff),o_topswai%name,itau_w, 847 $ topswai_aero) 848 ENDIF 849 IF (o_solswai%flag(iff)<=lev_files(iff)) THEN 850 CALL histwrite_phy(nid_files(iff),o_solswai%name,itau_w, 851 $ solswai_aero) 852 ENDIF 849 853 ENDIF 850 854 -
LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_state_var_mod.F90
r1054 r1150 255 255 !$OMP THREADPRIVATE(snow_con) 256 256 ! 257 ! sulfate_pi : SO4 aerosol concentration [ug/m3] (pre-industrial value)258 259 REAL,SAVE,ALLOCATABLE :: sulfate_pi(:, :)260 !$OMP THREADPRIVATE(sulfate_pi)261 257 REAL,SAVE,ALLOCATABLE :: rlonPOS(:) 262 258 !$OMP THREADPRIVATE(rlonPOS) … … 269 265 ! ok_aie=T -> 270 266 ! ok_ade=T -AIE=topswai-topswad 271 ! 267 ! ok_ade=F -AIE=topswai-topsw 272 268 ! 273 269 !topswad, solswad : Aerosol direct effect … … 277 273 REAL,SAVE,ALLOCATABLE :: topswai(:), solswai(:) 278 274 !$OMP THREADPRIVATE(topswai,solswai) 279 REAL,SAVE,ALLOCATABLE :: tau_ae(:,:,:), piz_ae(:,:,:) 280 !$OMP THREADPRIVATE(tau_ae,piz_ae) 281 REAL,SAVE,ALLOCATABLE :: cg_ae(:,:,:) 282 !$OMP THREADPRIVATE(cg_ae) 283 284 ! Les variables suivants uniquement pour un configuration avec INCA 285 ! topswad_inca, solswad_inca : Aerosol direct effect 286 REAL,SAVE,ALLOCATABLE :: topswad_inca(:), solswad_inca(:) 287 !$OMP THREADPRIVATE(topswad_inca,solswad_inca) 288 ! topswad0_inca, solswad0_inca : Aerosol direct effect 289 REAL,SAVE,ALLOCATABLE :: topswad0_inca(:), solswad0_inca(:) 290 !$OMP THREADPRIVATE(topswad0_inca,solswad0_inca) 291 ! topswai_inca, solswai_inca : Aerosol indirect effect 292 REAL,SAVE,ALLOCATABLE :: topswai_inca(:), solswai_inca(:) 293 !$OMP THREADPRIVATE(topswai_inca,solswai_inca) 294 REAL,SAVE,ALLOCATABLE :: topsw_inca(:,:), solsw_inca(:,:) 295 !$OMP THREADPRIVATE(topsw_inca,solsw_inca) 296 REAL,SAVE,ALLOCATABLE :: topsw0_inca(:,:), solsw0_inca(:,:) 297 !$OMP THREADPRIVATE(topsw0_inca,solsw0_inca) 298 REAL,SAVE,ALLOCATABLE :: tau_inca(:,:,:,:) 299 !$OMP THREADPRIVATE(tau_inca) 300 REAL,SAVE,ALLOCATABLE :: piz_inca(:,:,:,:) 301 !$OMP THREADPRIVATE(piz_inca) 302 REAL,SAVE,ALLOCATABLE :: cg_inca(:,:,:,:) 303 !$OMP THREADPRIVATE(cg_inca) 275 304 276 REAL,SAVE,ALLOCATABLE :: ccm(:,:,:) 305 277 !$OMP THREADPRIVATE(ccm) … … 402 374 ALLOCATE(ibas_con(klon), itop_con(klon)) 403 375 ALLOCATE(rain_con(klon), snow_con(klon)) 404 !405 ALLOCATE(sulfate_pi(klon, klev))406 376 ALLOCATE(rlonPOS(klon)) 407 377 ALLOCATE(newsst(klon)) … … 409 379 ALLOCATE(topswad(klon), solswad(klon)) 410 380 ALLOCATE(topswai(klon), solswai(klon)) 411 ALLOCATE(tau_ae(klon,klev,2), piz_ae(klon,klev,2)) 412 ALLOCATE(cg_ae(klon,klev,2)) 413 414 IF (config_inca /= 'none') THEN 415 ALLOCATE(topswad_inca(klon), solswad_inca(klon)) 416 ALLOCATE(topswad0_inca(klon), solswad0_inca(klon)) 417 ALLOCATE(topswai_inca(klon), solswai_inca(klon)) 418 ALLOCATE(topsw_inca(klon,9), solsw_inca(klon,9)) 419 ALLOCATE(topsw0_inca(klon,9), solsw0_inca(klon,9)) 420 END IF 421 ! Following 4 variables are needed only by INCA but must be 422 ! allocated as they exist in the phytrac argument list 423 ALLOCATE(tau_inca(klon,klev,9,2)) 424 ALLOCATE(piz_inca(klon,klev,9,2)) 425 ALLOCATE(cg_inca(klon,klev,9,2)) 381 426 382 ALLOCATE(ccm(klon,klev,2)) 427 383 … … 505 461 deallocate(ibas_con, itop_con) 506 462 deallocate(rain_con, snow_con) 507 !508 deallocate(sulfate_pi)509 463 deallocate(rlonPOS) 510 464 deallocate(newsst) … … 513 467 deallocate(topswai, solswai) 514 468 515 deallocate(tau_ae, piz_ae)516 deallocate(cg_ae)517 518 IF (config_inca /= 'none') THEN519 deallocate(topswad_inca, solswad_inca)520 deallocate(topswad0_inca, solswad0_inca)521 deallocate(topswai_inca, solswai_inca)522 deallocate(topsw_inca, solsw_inca)523 deallocate(topsw0_inca, solsw0_inca)524 END IF525 deallocate(tau_inca)526 deallocate(piz_inca)527 deallocate(cg_inca)528 469 deallocate(ccm) 529 470 -
LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F
r1134 r1150 50 50 c 51 51 c nlon----input-I-nombre de points horizontaux 52 c nlev----input-I-nombre de couches verticales 52 c nlev----input-I-nombre de couches verticales, doit etre egale a klev 53 53 c debut---input-L-variable logique indiquant le premier passage 54 54 c lafin---input-L-variable logique indiquant le dernier passage … … 736 736 EXTERNAL phyetat0 ! lire l'etat initial de la physique 737 737 EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique 738 EXTERNAL radlwsw ! rayonnements solaire et infrarouge739 738 EXTERNAL suphel ! initialiser certaines constantes 740 739 EXTERNAL transp ! transport total de l'eau et de l'energie … … 1056 1055 CHARACTER*40 t2mincels, t2maxcels !t2m min., t2m max 1057 1056 CHARACTER*40 tinst, tave, typeval 1058 cjq Aerosol effects (Johannes Quaas, 27/11/2003)1059 REAL sulfate(klon, klev) ! SO4 aerosol concentration [ug/m3]1060 1061 1057 REAL cldtaupi(klon,klev) ! Cloud optical thickness for pre-industrial (pi) aerosols 1062 1058 … … 1067 1063 1068 1064 ! Aerosol optical properties 1069 1070 ! Aerosol optical properties by INCA model 1071 CHARACTER*4 :: rfname(9) 1072 REAL aerindex(klon) ! POLDER aerosol index 1073 1065 CHARACTER*4, DIMENSION(9) :: rfname 1066 REAL, DIMENSION(klon) :: aerindex ! POLDER aerosol index 1067 REAL, DIMENSION(klon,klev) :: maerosol ! aerosol concentration [ug/m3] 1068 REAL, DIMENSION(klon,klev) :: maerosol_pi ! aerosol concentration [ug/m3] (pre-industrial value) 1069 REAL, DIMENSION(klon,klev,9,2) :: tau_aero, piz_aero, cg_aero 1070 REAL, DIMENSION(klon) :: topswad_aero, solswad_aero ! diag 1071 REAL, DIMENSION(klon) :: topswai_aero, solswai_aero ! diag 1072 REAL, DIMENSION(klon) :: topswad0_aero, solswad0_aero ! pas utilise, eventuellment pour diag 1073 REAL, DIMENSION(klon,9) :: topsw_aero, solsw_aero ! pas utilise 1074 REAL, DIMENSION(klon,9) :: topsw0_aero, solsw0_aero ! pas utilise 1075 1076 1074 1077 ! Parameters 1075 1078 LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not … … 1080 1083 ! false : lecture des aerosol dans un fichier 1081 1084 c$OMP THREADPRIVATE(aerosol_couple) 1082 1085 INTEGER, SAVE :: flag_aerosol 1086 c$OMP THREADPRIVATE(flag_aerosol) 1087 LOGICAL, SAVE :: new_aod 1088 c$OMP THREADPRIVATE(new_aod) 1089 1083 1090 c 1084 1091 c Declaration des constantes et des fonctions thermodynamiques … … 1106 1113 write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!' 1107 1114 write(lunout,*) 1108 s 'nlon, nlev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys'1115 s 'nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys' 1109 1116 write(lunout,*) 1110 s nlon, nlev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys1117 s nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys 1111 1118 1112 1119 write(lunout,*) 'papers, play, phi, u, v, t, omega' 1113 do k=1, nlev1120 do k=1,klev 1114 1121 write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k), 1115 1122 s u(igout,k),v(igout,k),t(igout,k),omega(igout,k) 1116 1123 enddo 1117 1124 write(lunout,*) 'ovap (g/kg), oliq (g/kg)' 1118 do k=1, nlev1125 do k=1,klev 1119 1126 write(lunout,*) qx(igout,k,1)*1000,qx(igout,k,2)*1000. 1120 1127 enddo … … 1190 1197 u10m(:,:)=0. 1191 1198 v10m(:,:)=0. 1192 piz_ae(:,:,:)=0.1193 tau_ae(:,:,:)=0.1194 cg_ae(:,:,:)=0.1195 1199 rain_con(:)=0. 1196 1200 snow_con(:)=0. … … 1205 1209 wmax_th(:)=0. 1206 1210 tau_overturning_th(:)=0. 1207 IF (config_inca /= 'none') THEN 1208 tau_inca(:,:,:,:) = 0. 1209 piz_inca(:,:,:,:) = 0. 1210 cg_inca(:,:,:,:) = 0. 1211 ccm(:,:,:) = 0. 1212 topswai_inca(:) = 0. 1213 topswad_inca(:) = 0. 1214 topswad0_inca(:) = 0. 1215 topsw_inca(:,:) = 0. 1216 topsw0_inca(:,:) = 0. 1217 solswai_inca(:) = 0. 1218 solswad_inca(:) = 0. 1219 solswad0_inca(:) = 0. 1220 solsw_inca(:,:) = 0. 1221 solsw0_inca(:,:) = 0. 1222 END IF 1211 1212 IF (config_inca /= 'none') ccm(:,:,:) = 0. 1223 1213 1224 1214 rnebcon0(:,:) = 0.0 … … 1239 1229 . iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut, 1240 1230 . ok_ade, ok_aie, aerosol_couple, 1231 . flag_aerosol, new_aod, 1241 1232 . bl95_b0, bl95_b1, 1242 1233 . iflag_thermals,nsplit_thermals,tau_thermals, … … 2658 2649 cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr) 2659 2650 IF (ok_ade.OR.ok_aie) THEN 2660 IF ( .NOT. aerosol_couple ) THEN 2661 ! Get sulfate aerosol distribution 2662 CALL readsulfate(rjourvrai, debut, sulfate) 2663 CALL readsulfate_preind(rjourvrai, debut, sulfate_pi) 2664 2665 ! Calculate aerosol optical properties (Olivier Boucher) 2666 CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl, 2667 . tau_ae, piz_ae, cg_ae, aerindex) 2668 ENDIF 2651 IF (.NOT. aerosol_couple) 2652 & CALL aerosol_optic( 2653 & debut, new_aod, flag_aerosol, rjourvrai, pdtphys, 2654 & pplay, paprs, t_seri, rhcl, 2655 & maerosol, maerosol_pi, 2656 & tau_aero, piz_aero, cg_aero ) 2669 2657 ELSE 2670 tau_ae(:,:,:)=0.02671 piz_ae(:,:,:)=0.0 2672 cg_ae(:,:,:)=0.0 2658 tau_aero(:,:,:,:) = 0. 2659 piz_aero(:,:,:,:) = 0. 2660 cg_aero(:,:,:,:) = 0. 2673 2661 ENDIF 2674 2662 … … 2847 2835 2848 2836 IF (aerosol_couple) THEN 2849 sulfate(:,:)= ccm(:,:,1)2850 sulfate_pi(:,:) = ccm(:,:,2)2851 END IF2837 maerosol(:,:) = ccm(:,:,1) 2838 maerosol_pi(:,:) = ccm(:,:,2) 2839 END IF 2852 2840 2853 2841 if (ok_newmicro) then … … 2857 2845 . flwp, fiwp, flwc, fiwc, 2858 2846 e ok_aie, 2859 e sulfate, sulfate_pi,2847 e maerosol, maerosol_pi, 2860 2848 e bl95_b0, bl95_b1, 2861 2849 s cldtaupi, re, fl) … … 2865 2853 . cldh, cldl, cldm, cldt, cldq, 2866 2854 e ok_aie, 2867 e sulfate, sulfate_pi,2855 e maerosol, maerosol_pi, 2868 2856 e bl95_b0, bl95_b1, 2869 2857 s cldtaupi, re, fl) … … 2895 2883 IF (aerosol_couple) THEN 2896 2884 #ifdef INCA 2897 CALL radlwsw_inca 2898 e (kdlon,kflev,dist, rmu0, fract, solaire, 2899 e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, 2900 e wo, 2901 e cldfra, cldemi, cldtau, 2902 s heat,heat0,cool,cool0,radsol,albpla, 2903 s topsw,toplw,solsw,sollw, 2904 s sollwdown, 2905 s topsw0,toplw0,solsw0,sollw0, 2906 s lwdn0, lwdn, lwup0, lwup, 2907 s swdn0, swdn, swup0, swup, 2908 e ok_ade, ok_aie, 2909 e tau_inca, piz_inca, cg_inca, 2910 s topswad_inca, solswad_inca, 2911 s topswad0_inca, solswad0_inca, 2912 s topsw_inca, topsw0_inca, 2913 s solsw_inca, solsw0_inca, 2914 e cldtaupi, 2915 s topswai_inca, solswai_inca) 2885 CALL radlwsw_inca 2886 e (kdlon,kflev,dist, rmu0, fract, solaire, 2887 e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, 2888 e wo, 2889 e cldfra, cldemi, cldtau, 2890 s heat,heat0,cool,cool0,radsol,albpla, 2891 s topsw,toplw,solsw,sollw, 2892 s sollwdown, 2893 s topsw0,toplw0,solsw0,sollw0, 2894 s lwdn0, lwdn, lwup0, lwup, 2895 s swdn0, swdn, swup0, swup, 2896 e ok_ade, ok_aie, 2897 e tau_aero, piz_aero, cg_aero, 2898 s topswad_aero, solswad_aero, 2899 s topswad0_aero, solswad0_aero, 2900 s topsw_aero, topsw0_aero, 2901 s solsw_aero, solsw0_aero, 2902 e cldtaupi, 2903 s topswai_aero, solswai_aero) 2904 2916 2905 #endif 2917 2906 ELSE 2918 CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS) 2919 e (dist, rmu0, fract, 2920 e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, 2921 e wo, 2922 e cldfra, cldemi, cldtau, 2923 s heat,heat0,cool,cool0,radsol,albpla, 2924 s topsw,toplw,solsw,sollw, 2925 s sollwdown, 2926 s topsw0,toplw0,solsw0,sollw0, 2927 s lwdn0, lwdn, lwup0, lwup, 2928 s swdn0, swdn, swup0, swup, 2929 e ok_ade, ok_aie, ! new for aerosol radiative effects 2930 e tau_ae, piz_ae, cg_ae, ! ="= 2931 s topswad, solswad, ! ="= 2932 e cldtaupi, ! ="= 2933 s topswai, solswai,zqsat,flwc,fiwc) ! ="= 2907 2908 CALL radlwsw_aero 2909 e (dist, rmu0, fract, solaire, 2910 e paprs, pplay,zxtsol,albsol1, albsol2, 2911 e t_seri,q_seri,wo, 2912 e cldfra, cldemi, cldtau, 2913 e ok_ade, ok_aie, 2914 e tau_aero, piz_aero, cg_aero, 2915 e cldtaupi,new_aod, 2916 s heat,heat0,cool,cool0,radsol,albpla, 2917 s topsw,toplw,solsw,sollw, 2918 s sollwdown, 2919 s topsw0,toplw0,solsw0,sollw0, 2920 s lwdn0, lwdn, lwup0, lwup, 2921 s swdn0, swdn, swup0, swup, 2922 s topswad_aero, solswad_aero, 2923 s topswai_aero, solswai_aero, 2924 o topswad0_aero, solswad0_aero, 2925 o topsw_aero, topsw0_aero, 2926 o solsw_aero, solsw0_aero) 2927 2928 2934 2929 ENDIF 2935 2930 itaprad = 0 … … 3159 3154 I lafin, 3160 3155 I nlon, 3161 I nlev,3156 I klev, 3162 3157 I dtime, 3163 3158 I u, … … 3207 3202 I aerosol_couple, 3208 3203 I flxmass_w, 3209 I tau_ inca,3210 I piz_ inca,3211 I cg_ inca,3204 I tau_aero, 3205 I piz_aero, 3206 I cg_aero, 3212 3207 I ccm, 3213 3208 I rfname, … … 3218 3213 print*,'Attention on met a 0 les thermiques pour phystoke' 3219 3214 call phystokenc ( 3220 I nlon, nlev,pdtphys,rlon,rlat,3215 I nlon,klev,pdtphys,rlon,rlat, 3221 3216 I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, 3222 3217 I fm_therm,entr_therm, … … 3415 3410 write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!' 3416 3411 write(lunout,*) 3417 s 'nlon, nlev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys pct tlos'3412 s 'nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys pct tlos' 3418 3413 write(lunout,*) 3419 s nlon, nlev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys,3414 s nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys, 3420 3415 s pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), 3421 3416 s pctsrf(igout,is_sic) 3422 3417 write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva' 3423 do k=1, nlev3418 do k=1,klev 3424 3419 write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), 3425 3420 s d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), … … 3427 3422 enddo 3428 3423 write(lunout,*) 'cool,heat' 3429 do k=1, nlev3424 do k=1,klev 3430 3425 write(lunout,*) cool(igout,k),heat(igout,k) 3431 3426 enddo 3432 3427 3433 3428 write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec' 3434 do k=1, nlev3429 do k=1,klev 3435 3430 write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), 3436 3431 s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k) … … 3439 3434 write(lunout,*) 'd_ps ',d_ps(igout) 3440 3435 write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 ' 3441 do k=1, nlev3436 do k=1,klev 3442 3437 write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), 3443 3438 s d_qx(igout,k,1),d_qx(igout,k,2) -
LMDZ4/branches/LMDZ4-dev/libf/phylmd/phytrac.F
r1141 r1150 57 57 I aerosol_couple, 58 58 I flxmass_w, 59 I tau_ inca,60 I piz_ inca,61 I cg_ inca,59 I tau_aero, 60 I piz_aero, 61 I cg_aero, 62 62 I ccm, 63 63 I rfname, … … 138 138 CHARACTER(len=8) :: solsym(nbtr) 139 139 integer la 140 REAL :: tau_ inca(klon,klev,9,2)141 REAL :: piz_ inca(klon,klev,9,2)142 REAL :: cg_ inca(klon,klev,9,2)140 REAL :: tau_aero(klon,klev,9,2) 141 REAL :: piz_aero(klon,klev,9,2) 142 REAL :: cg_aero(klon,klev,9,2) 143 143 character*4 :: rfname(9) 144 144 REAL :: ccm(klon,klev,2) … … 458 458 $ t_seri, ! for chimiaq 459 459 $ rh, 460 $ tau_ inca,461 $ piz_ inca,462 $ cg_ inca,460 $ tau_aero, 461 $ piz_aero, 462 $ cg_aero, 463 463 $ rfname, 464 464 $ ccm, … … 497 497 $ sh, !sh 498 498 $ rh, !rh 499 $ .false., !wrhstts500 $ hbuf, !hbuf501 $ obuf, !obuf502 499 $ iip1, !nx 503 500 $ jjp1, !ny -
LMDZ4/branches/LMDZ4-dev/libf/phylmd/radlwsw_aero.F90
r1149 r1150 1 ! 2 ! $Header$ 3 ! 4 SUBROUTINE radlwsw(dist, rmu0, fract, 5 . paprs, pplay,tsol,alb1, alb2, t,q,wo, 6 . cldfra, cldemi, cldtaupd, 7 . heat,heat0,cool,cool0,radsol,albpla, 8 . topsw,toplw,solsw,sollw, 9 . sollwdown, 10 . topsw0,toplw0,solsw0,sollw0, 11 . lwdn0, lwdn, lwup0, lwup, 12 . swdn0, swdn, swup0, swup, 13 . ok_ade, ok_aie, 14 . tau_ae, piz_ae, cg_ae, 15 . topswad, solswad, 16 . cldtaupi, topswai, solswai,qsat,flwc,fiwc) 17 c 18 USE dimphy 19 IMPLICIT none 20 c====================================================================== 21 c Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719 22 c Objet: interface entre le modele et les rayonnements 23 c Arguments: 24 c dist-----input-R- distance astronomique terre-soleil 25 c rmu0-----input-R- cosinus de l'angle zenithal 26 c fract----input-R- duree d'ensoleillement normalisee 27 c co2_ppm--input-R- concentration du gaz carbonique (en ppm) 28 c solaire--input-R- constante solaire (W/m**2) 29 c paprs----input-R- pression a inter-couche (Pa) 30 c pplay----input-R- pression au milieu de couche (Pa) 31 c tsol-----input-R- temperature du sol (en K) 32 c alb1-----input-R- albedo du sol(entre 0 et 1) dans l'interval visible 33 c alb2-----input-R- albedo du sol(entre 0 et 1) dans l'interval proche infra-rouge 34 c t--------input-R- temperature (K) 35 c q--------input-R- vapeur d'eau (en kg/kg) 36 c wo-------input-R- contenu en ozone (en kg/kg) correction MPL 100505 37 c cldfra---input-R- fraction nuageuse (entre 0 et 1) 38 c cldtaupd---input-R- epaisseur optique des nuages dans le visible (present-day value) 39 c cldemi---input-R- emissivite des nuages dans l'IR (entre 0 et 1) 40 c ok_ade---input-L- apply the Aerosol Direct Effect or not? 41 c ok_aie---input-L- apply the Aerosol Indirect Effect or not? 42 c tau_ae, piz_ae, cg_ae-input-R- aerosol optical properties (calculated in aeropt.F) 43 c cldtaupi-input-R- epaisseur optique des nuages dans le visible 44 c calculated for pre-industrial (pi) aerosol concentrations, i.e. with smaller 45 c droplet concentration, thus larger droplets, thus generally cdltaupi cldtaupd 46 c it is needed for the diagnostics of the aerosol indirect radiative forcing 47 c 48 c heat-----output-R- echauffement atmospherique (visible) (K/jour) 49 c cool-----output-R- refroidissement dans l'IR (K/jour) 50 c radsol---output-R- bilan radiatif net au sol (W/m**2) (+ vers le bas) 51 c albpla---output-R- albedo planetaire (entre 0 et 1) 52 c topsw----output-R- flux solaire net au sommet de l'atm. 53 c toplw----output-R- ray. IR montant au sommet de l'atmosphere 54 c solsw----output-R- flux solaire net a la surface 55 c sollw----output-R- ray. IR montant a la surface 56 c solswad---output-R- ray. solaire net absorbe a la surface (aerosol dir) 57 c topswad---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol dir) 58 c solswai---output-R- ray. solaire net absorbe a la surface (aerosol ind) 59 c topswai---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol ind) 60 c 61 c ATTENTION: swai and swad have to be interpreted in the following manner: 62 c --------- 63 c ok_ade=F & ok_aie=F -both are zero 64 c ok_ade=T & ok_aie=F -aerosol direct forcing is F_{AD} = topsw-topswad 65 c indirect is zero 66 c ok_ade=F & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai 67 c direct is zero 68 c ok_ade=T & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai 69 c aerosol direct forcing is F_{AD} = topswai-topswad 70 c 71 72 c====================================================================== 73 cym#include "dimensions.h" 74 cym#include "dimphy.h" 75 cym#include "raddim.h" 76 #include "YOETHF.h" 77 c 78 real rmu0(klon), fract(klon), dist 79 cIM real co2_ppm 80 cIM real solaire 81 #include "clesphys.h" 82 c 83 real paprs(klon,klev+1), pplay(klon,klev) 84 real alb1(klon), alb2(klon), tsol(klon) 85 real t(klon,klev), q(klon,klev), wo(klon,klev) 86 real cldfra(klon,klev), cldemi(klon,klev), cldtaupd(klon,klev) 87 real heat(klon,klev), cool(klon,klev) 88 real heat0(klon,klev), cool0(klon,klev) 89 real radsol(klon), topsw(klon), toplw(klon) 90 real solsw(klon), sollw(klon), albpla(klon) 91 real topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon) 92 real sollwdown(klon) 93 cIM output 3D 94 REAL*8 ZFSUP(KDLON,KFLEV+1) 95 REAL*8 ZFSDN(KDLON,KFLEV+1) 96 REAL*8 ZFSUP0(KDLON,KFLEV+1) 97 REAL*8 ZFSDN0(KDLON,KFLEV+1) 98 c 99 REAL*8 ZFLUP(KDLON,KFLEV+1) 100 REAL*8 ZFLDN(KDLON,KFLEV+1) 101 REAL*8 ZFLUP0(KDLON,KFLEV+1) 102 REAL*8 ZFLDN0(KDLON,KFLEV+1) 103 c 104 REAL*8 zx_alpha1, zx_alpha2 105 c 106 #include "YOMCST.h" 107 c 108 INTEGER k, kk, i, j, iof, nb_gr 109 EXTERNAL LW_LMDAR4,SW_LMDAR4 110 c 111 cIM ctes ds clesphys.h REAL*8 RCO2, RCH4, RN2O, RCFC11, RCFC12 112 REAL*8 PSCT 113 c 114 REAL*8 PALBD(kdlon,2), PALBP(kdlon,2) 115 REAL*8 PEMIS(kdlon), PDT0(kdlon), PVIEW(kdlon) 116 REAL*8 PPSOL(kdlon), PDP(kdlon,klev) 117 REAL*8 PTL(kdlon,kflev+1), PPMB(kdlon,kflev+1) 118 REAL*8 PTAVE(kdlon,kflev) 119 REAL*8 PWV(kdlon,kflev), PQS(kdlon,kflev), POZON(kdlon,kflev) 120 REAL*8 PAER(kdlon,kflev,5) 121 REAL*8 PCLDLD(kdlon,kflev) 122 REAL*8 PCLDLU(kdlon,kflev) 123 REAL*8 PCLDSW(kdlon,kflev) 124 REAL*8 PTAU(kdlon,2,kflev) 125 REAL*8 POMEGA(kdlon,2,kflev) 126 REAL*8 PCG(kdlon,2,kflev) 127 c 128 REAL*8 zfract(kdlon), zrmu0(kdlon), zdist 129 c 130 REAL*8 zheat(kdlon,kflev), zcool(kdlon,kflev) 131 REAL*8 zheat0(kdlon,kflev), zcool0(kdlon,kflev) 132 REAL*8 ztopsw(kdlon), ztoplw(kdlon) 133 REAL*8 zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon) 134 cIM 135 REAL*8 zsollwdown(kdlon) 136 c 137 REAL*8 ztopsw0(kdlon), ztoplw0(kdlon) 138 REAL*8 zsolsw0(kdlon), zsollw0(kdlon) 139 REAL*8 zznormcp 140 cIM output 3D : SWup, SWdn, LWup, LWdn 141 REAL swdn(klon,kflev+1),swdn0(klon,kflev+1) 142 REAL swup(klon,kflev+1),swup0(klon,kflev+1) 143 REAL lwdn(klon,kflev+1),lwdn0(klon,kflev+1) 144 REAL lwup(klon,kflev+1),lwup0(klon,kflev+1) 145 REAL qsat(klon,klev),flwc(klon,klev),fiwc(klon,klev) 146 c-OB 147 cjq the following quantities are needed for the aerosol radiative forcings 148 149 real topswad(klon), solswad(klon) ! output: aerosol direct forcing at TOA and surface 150 real topswai(klon), solswai(klon) ! output: aerosol indirect forcing atTOA and surface 151 real tau_ae(klon,klev,2), piz_ae(klon,klev,2), cg_ae(klon,klev,2) ! aerosol optical properties (see aeropt.F) 152 real cldtaupi(klon,klev) ! cloud optical thickness for pre-industrial aerosol concentrations 153 ! (i.e., with a smaller droplet concentrationand thus larger droplet radii) 154 logical ok_ade, ok_aie ! switches whether to use aerosol direct (indirect) effects or not 155 real*8 tauae(kdlon,kflev,2) ! aer opt properties 156 real*8 pizae(kdlon,kflev,2) 157 real*8 cgae(kdlon,kflev,2) 158 REAL*8 PTAUA(kdlon,2,kflev) ! present-day value of cloud opt thickness (PTAU is pre-industrial value), local use 159 REAL*8 POMEGAA(kdlon,2,kflev) ! dito for single scatt albedo 160 REAL*8 ztopswad(kdlon), zsolswad(kdlon) ! Aerosol direct forcing at TOAand surface 161 REAL*8 ztopswai(kdlon), zsolswai(kdlon) ! dito, indirect 162 cjq-end 163 !rv 164 tauae(:,:,:)=0. 165 pizae(:,:,:)=0. 166 cgae(:,:,:)=0. 167 !rv 168 169 c 170 c------------------------------------------- 171 nb_gr = klon / kdlon 172 IF (nb_gr*kdlon .NE. klon) THEN 173 PRINT*, "kdlon mauvais:", klon, kdlon, nb_gr 174 CALL abort 175 ENDIF 176 IF (kflev .NE. klev) THEN 177 PRINT*, "kflev differe de klev, kflev, klev" 178 CALL abort 179 ENDIF 180 c------------------------------------------- 181 DO k = 1, klev 182 DO i = 1, klon 183 heat(i,k)=0. 184 cool(i,k)=0. 185 heat0(i,k)=0. 186 cool0(i,k)=0. 187 ENDDO 188 ENDDO 189 c 190 zdist = dist 191 c 192 cIM anciennes valeurs 193 c RCO2 = co2_ppm * 1.0e-06 * 44.011/28.97 194 c 195 cIM : on met RCO2, RCH4, RN2O, RCFC11 et RCFC12 dans clesphys.h /lecture ds conf_phys.F90 196 c RCH4 = 1.65E-06* 16.043/28.97 197 c RN2O = 306.E-09* 44.013/28.97 198 c RCFC11 = 280.E-12* 137.3686/28.97 199 c RCFC12 = 484.E-12* 120.9140/28.97 200 cIM anciennes valeurs 201 c RCH4 = 1.72E-06* 16.043/28.97 202 c RN2O = 310.E-09* 44.013/28.97 203 c 204 c PRINT*,'IMradlwsw : solaire, co2= ', solaire, co2_ppm 205 PSCT = solaire/zdist/zdist 206 c 207 DO 99999 j = 1, nb_gr 208 iof = kdlon*(j-1) 209 c 1 SUBROUTINE radlwsw_aero( & 2 dist, rmu0, fract, solaire, & 3 paprs, pplay,tsol,albedo, alblw, & 4 t,q,wo,& 5 cldfra, cldemi, cldtaupd,& 6 ok_ade, ok_aie,& 7 tau_aero, piz_aero, cg_aero,& 8 cldtaupi, new_aod, & 9 heat,heat0,cool,cool0,radsol,albpla,& 10 topsw,toplw,solsw,sollw,& 11 sollwdown,& 12 topsw0,toplw0,solsw0,sollw0,& 13 lwdn0, lwdn, lwup0, lwup,& 14 swdn0, swdn, swup0, swup,& 15 topswad_aero, solswad_aero,& 16 topswai_aero, solswai_aero, & 17 topswad0_aero, solswad0_aero,& 18 topsw_aero, topsw0_aero,& 19 solsw_aero, solsw0_aero) 20 21 22 23 USE DIMPHY 24 25 IMPLICIT NONE 26 !====================================================================== 27 ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719 28 ! Objet: interface entre le modele et les rayonnements 29 ! Arguments: 30 ! dist-----input-R- distance astronomique terre-soleil 31 ! rmu0-----input-R- cosinus de l'angle zenithal 32 ! fract----input-R- duree d'ensoleillement normalisee 33 ! co2_ppm--input-R- concentration du gaz carbonique (en ppm) 34 ! solaire--input-R- constante solaire (W/m**2) 35 ! paprs----input-R- pression a inter-couche (Pa) 36 ! pplay----input-R- pression au milieu de couche (Pa) 37 ! tsol-----input-R- temperature du sol (en K) 38 ! albedo---input-R- albedo du sol (entre 0 et 1) 39 ! t--------input-R- temperature (K) 40 ! q--------input-R- vapeur d'eau (en kg/kg) 41 ! wo-------input-R- contenu en ozone (en kg/kg) correction MPL 100505 42 ! cldfra---input-R- fraction nuageuse (entre 0 et 1) 43 ! cldtaupd---input-R- epaisseur optique des nuages dans le visible (present-day value) 44 ! cldemi---input-R- emissivite des nuages dans l'IR (entre 0 et 1) 45 ! ok_ade---input-L- apply the Aerosol Direct Effect or not? 46 ! ok_aie---input-L- apply the Aerosol Indirect Effect or not? 47 ! tau_ae, piz_ae, cg_ae-input-R- aerosol optical properties (calculated in aeropt.F) 48 ! cldtaupi-input-R- epaisseur optique des nuages dans le visible 49 ! calculated for pre-industrial (pi) aerosol concentrations, i.e. with smaller 50 ! droplet concentration, thus larger droplets, thus generally cdltaupi cldtaupd 51 ! it is needed for the diagnostics of the aerosol indirect radiative forcing 52 ! 53 ! heat-----output-R- echauffement atmospherique (visible) (K/jour) 54 ! cool-----output-R- refroidissement dans l'IR (K/jour) 55 ! radsol---output-R- bilan radiatif net au sol (W/m**2) (+ vers le bas) 56 ! albpla---output-R- albedo planetaire (entre 0 et 1) 57 ! topsw----output-R- flux solaire net au sommet de l'atm. 58 ! toplw----output-R- ray. IR montant au sommet de l'atmosphere 59 ! solsw----output-R- flux solaire net a la surface 60 ! sollw----output-R- ray. IR montant a la surface 61 ! solswad---output-R- ray. solaire net absorbe a la surface (aerosol dir) 62 ! topswad---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol dir) 63 ! solswai---output-R- ray. solaire net absorbe a la surface (aerosol ind) 64 ! topswai---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol ind) 65 ! 66 ! ATTENTION: swai and swad have to be interpreted in the following manner: 67 ! --------- 68 ! ok_ade=F & ok_aie=F -both are zero 69 ! ok_ade=T & ok_aie=F -aerosol direct forcing is F_{AD} = topsw-topswad 70 ! indirect is zero 71 ! ok_ade=F & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai 72 ! direct is zero 73 ! ok_ade=T & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai 74 ! aerosol direct forcing is F_{AD} = topswai-topswad 75 ! 76 77 !====================================================================== 78 79 ! ==================================================================== 80 ! Adapte au modele de chimie INCA par Celine Deandreis & Anne Cozic -- 2009 81 ! 1 = ZERO 82 ! 2 = AER total 83 ! 3 = NAT 84 ! 4 = BC 85 ! 5 = SO4 86 ! 6 = POM 87 ! 7 = DUST 88 ! 8 = SS 89 ! 9 = NO3 90 ! 91 ! ==================================================================== 92 include "YOETHF.h" 93 include "YOMCST.h" 94 95 ! Input arguments 96 REAL, INTENT(in) :: solaire 97 REAL, INTENT(in) :: dist 98 REAL, INTENT(in) :: rmu0(KLON), fract(KLON) 99 REAL, INTENT(in) :: paprs(KLON,KLEV+1), pplay(KLON,KLEV) 100 REAL, INTENT(in) :: albedo(KLON), alblw(KLON), tsol(KLON) 101 REAL, INTENT(in) :: t(KLON,KLEV), q(KLON,KLEV), wo(KLON,KLEV) 102 LOGICAL, INTENT(in) :: ok_ade, ok_aie ! switches whether to use aerosol direct (indirect) effects or not 103 REAL, INTENT(in) :: cldfra(KLON,KLEV), cldemi(KLON,KLEV), cldtaupd(KLON,KLEV) 104 REAL, INTENT(in) :: tau_aero(KLON,KLEV,9,2) ! aerosol optical properties (see aeropt.F) 105 REAL, INTENT(in) :: piz_aero(KLON,KLEV,9,2) ! aerosol optical properties (see aeropt.F) 106 REAL, INTENT(in) :: cg_aero(KLON,KLEV,9,2) ! aerosol optical properties (see aeropt.F) 107 REAL, INTENT(in) :: cldtaupi(KLON,KLEV) ! cloud optical thickness for pre-industrial aerosol concentrations 108 LOGICAL, INTENT(in) :: new_aod ! flag pour retrouver les resultats exacts de l'AR4 dans le cas ou l'on ne travaille qu'avec les sulfates 109 110 ! Output arguments 111 REAL, INTENT(out) :: heat(KLON,KLEV), cool(KLON,KLEV) 112 REAL, INTENT(out) :: heat0(KLON,KLEV), cool0(KLON,KLEV) 113 REAL, INTENT(out) :: radsol(KLON), topsw(KLON), toplw(KLON) 114 REAL, INTENT(out) :: solsw(KLON), sollw(KLON), albpla(KLON) 115 REAL, INTENT(out) :: topsw0(KLON), toplw0(KLON), solsw0(KLON), sollw0(KLON) 116 REAL, INTENT(out) :: sollwdown(KLON) 117 REAL, INTENT(out) :: swdn(KLON,kflev+1),swdn0(KLON,kflev+1) 118 REAL, INTENT(out) :: swup(KLON,kflev+1),swup0(KLON,kflev+1) 119 REAL, INTENT(out) :: lwdn(KLON,kflev+1),lwdn0(KLON,kflev+1) 120 REAL, INTENT(out) :: lwup(KLON,kflev+1),lwup0(KLON,kflev+1) 121 REAL, INTENT(out) :: topswad_aero(KLON), solswad_aero(KLON) ! output: aerosol direct forcing at TOA and surface 122 REAL, INTENT(out) :: topswai_aero(KLON), solswai_aero(KLON) ! output: aerosol indirect forcing atTOA and surface 123 REAL, DIMENSION(klon), INTENT(out) :: topswad0_aero 124 REAL, DIMENSION(klon), INTENT(out) :: solswad0_aero 125 REAL, DIMENSION(kdlon,9), INTENT(out) :: topsw_aero 126 REAL, DIMENSION(kdlon,9), INTENT(out) :: topsw0_aero 127 REAL, DIMENSION(kdlon,9), INTENT(out) :: solsw_aero 128 REAL, DIMENSION(kdlon,9), INTENT(out) :: solsw0_aero 129 130 ! Local variables 131 REAL*8 ZFSUP(KDLON,KFLEV+1) 132 REAL*8 ZFSDN(KDLON,KFLEV+1) 133 REAL*8 ZFSUP0(KDLON,KFLEV+1) 134 REAL*8 ZFSDN0(KDLON,KFLEV+1) 135 REAL*8 ZFLUP(KDLON,KFLEV+1) 136 REAL*8 ZFLDN(KDLON,KFLEV+1) 137 REAL*8 ZFLUP0(KDLON,KFLEV+1) 138 REAL*8 ZFLDN0(KDLON,KFLEV+1) 139 REAL*8 zx_alpha1, zx_alpha2 140 INTEGER k, kk, i, j, iof, nb_gr 141 REAL*8 PSCT 142 REAL*8 PALBD(kdlon,2), PALBP(kdlon,2) 143 REAL*8 PEMIS(kdlon), PDT0(kdlon), PVIEW(kdlon) 144 REAL*8 PPSOL(kdlon), PDP(kdlon,KLEV) 145 REAL*8 PTL(kdlon,kflev+1), PPMB(kdlon,kflev+1) 146 REAL*8 PTAVE(kdlon,kflev) 147 REAL*8 PWV(kdlon,kflev), PQS(kdlon,kflev), POZON(kdlon,kflev) 148 REAL*8 PAER(kdlon,kflev,5) 149 REAL*8 PCLDLD(kdlon,kflev) 150 REAL*8 PCLDLU(kdlon,kflev) 151 REAL*8 PCLDSW(kdlon,kflev) 152 REAL*8 PTAU(kdlon,2,kflev) 153 REAL*8 POMEGA(kdlon,2,kflev) 154 REAL*8 PCG(kdlon,2,kflev) 155 REAL*8 zfract(kdlon), zrmu0(kdlon), zdist 156 REAL*8 zheat(kdlon,kflev), zcool(kdlon,kflev) 157 REAL*8 zheat0(kdlon,kflev), zcool0(kdlon,kflev) 158 REAL*8 ztopsw(kdlon), ztoplw(kdlon) 159 REAL*8 zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon) 160 REAL*8 zsollwdown(kdlon) 161 REAL*8 ztopsw0(kdlon), ztoplw0(kdlon) 162 REAL*8 zsolsw0(kdlon), zsollw0(kdlon) 163 REAL*8 zznormcp 164 REAL*8 tauaero(kdlon,kflev,9,2) ! aer opt properties 165 REAL*8 pizaero(kdlon,kflev,9,2) 166 REAL*8 cgaero(kdlon,kflev,9,2) 167 REAL*8 PTAUA(kdlon,2,kflev) ! present-day value of cloud opt thickness (PTAU is pre-industrial value), local use 168 REAL*8 POMEGAA(kdlon,2,kflev) ! dito for single scatt albedo 169 REAL*8 ztopswadaero(kdlon), zsolswadaero(kdlon) ! Aerosol direct forcing at TOAand surface 170 REAL*8 ztopswad0aero(kdlon), zsolswad0aero(kdlon) ! Aerosol direct forcing at TOAand surface 171 REAL*8 ztopswaiaero(kdlon), zsolswaiaero(kdlon) ! dito, indirect 172 REAL*8 ztopsw_aero(kdlon,9), ztopsw0_aero(kdlon,9) 173 REAL*8 zsolsw_aero(kdlon,9), zsolsw0_aero(kdlon,9) 174 175 ! initialisation 176 tauaero(:,:,:,:)=0. 177 pizaero(:,:,:,:)=0. 178 cgaero(:,:,:,:)=0. 179 180 ! 181 !------------------------------------------- 182 nb_gr = KLON / kdlon 183 IF (nb_gr*kdlon .NE. KLON) THEN 184 PRINT*, "kdlon mauvais:", KLON, kdlon, nb_gr 185 CALL abort 186 ENDIF 187 IF (kflev .NE. KLEV) THEN 188 PRINT*, "kflev differe de KLEV, kflev, KLEV" 189 CALL abort 190 ENDIF 191 !------------------------------------------- 192 DO k = 1, KLEV 193 DO i = 1, KLON 194 heat(i,k)=0. 195 cool(i,k)=0. 196 heat0(i,k)=0. 197 cool0(i,k)=0. 198 ENDDO 199 ENDDO 200 ! 201 zdist = dist 202 ! 203 PSCT = solaire/zdist/zdist 204 DO j = 1, nb_gr 205 iof = kdlon*(j-1) 206 DO i = 1, kdlon 207 zfract(i) = fract(iof+i) 208 zrmu0(i) = rmu0(iof+i) 209 PALBD(i,1) = albedo(iof+i) 210 PALBD(i,2) = alblw(iof+i) 211 PALBP(i,1) = albedo(iof+i) 212 PALBP(i,2) = alblw(iof+i) 213 PEMIS(i) = 1.0 214 PVIEW(i) = 1.66 215 PPSOL(i) = paprs(iof+i,1) 216 zx_alpha1 = (paprs(iof+i,1)-pplay(iof+i,2))/(pplay(iof+i,1)-pplay(iof+i,2)) 217 zx_alpha2 = 1.0 - zx_alpha1 218 PTL(i,1) = t(iof+i,1) * zx_alpha1 + t(iof+i,2) * zx_alpha2 219 PTL(i,KLEV+1) = t(iof+i,KLEV) 220 PDT0(i) = tsol(iof+i) - PTL(i,1) 221 ENDDO 222 DO k = 2, kflev 210 223 DO i = 1, kdlon 211 zfract(i) = fract(iof+i) 212 zrmu0(i) = rmu0(iof+i) 213 PALBD(i,1) = alb1(iof+i) 214 ! PALBD(i,2) = alb1(iof+i) 215 PALBD(i,2) = alb2(iof+i) 216 PALBP(i,1) = alb1(iof+i) 217 ! PALBP(i,2) = alb1(iof+i) 218 PALBP(i,2) = alb2(iof+i) 219 cIM cf. JLD pour etre en accord avec ORCHIDEE il faut mettre PEMIS(i) = 0.96 220 PEMIS(i) = 1.0 221 PVIEW(i) = 1.66 222 PPSOL(i) = paprs(iof+i,1) 223 zx_alpha1 = (paprs(iof+i,1)-pplay(iof+i,2)) 224 . / (pplay(iof+i,1)-pplay(iof+i,2)) 225 zx_alpha2 = 1.0 - zx_alpha1 226 PTL(i,1) = t(iof+i,1) * zx_alpha1 + t(iof+i,2) * zx_alpha2 227 PTL(i,klev+1) = t(iof+i,klev) 228 PDT0(i) = tsol(iof+i) - PTL(i,1) 229 ENDDO 230 DO k = 2, kflev 224 PTL(i,k) = (t(iof+i,k)+t(iof+i,k-1))*0.5 225 ENDDO 226 ENDDO 227 DO k = 1, kflev 231 228 DO i = 1, kdlon 232 PTL(i,k) = (t(iof+i,k)+t(iof+i,k-1))*0.5 233 ENDDO 234 ENDDO 229 PDP(i,k) = paprs(iof+i,k)-paprs(iof+i,k+1) 230 PTAVE(i,k) = t(iof+i,k) 231 PWV(i,k) = MAX (q(iof+i,k), 1.0e-12) 232 PQS(i,k) = PWV(i,k) 233 ! wo: cm.atm (epaisseur en cm dans la situation standard) 234 ! POZON: kg/kg 235 POZON(i,k) = MAX(wo(iof+i,k),1.0e-12)*RG/46.6968 & 236 /(paprs(iof+i,k)-paprs(iof+i,k+1))& 237 *(paprs(iof+i,1)/101325.0) 238 PCLDLD(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k) 239 PCLDLU(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k) 240 PCLDSW(i,k) = cldfra(iof+i,k) 241 PTAU(i,1,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! 1e-12 serait instable 242 PTAU(i,2,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! pour 32-bit machines 243 POMEGA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAU(i,1,k)) 244 POMEGA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i,2,k)) 245 PCG(i,1,k) = 0.865 246 PCG(i,2,k) = 0.910 247 !- 248 ! Introduced for aerosol indirect forcings. 249 ! The following values use the cloud optical thickness calculated from 250 ! present-day aerosol concentrations whereas the quantities without the 251 ! "A" at the end are for pre-industial (natural-only) aerosol concentrations 252 ! 253 PTAUA(i,1,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! 1e-12 serait instable 254 PTAUA(i,2,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! pour 32-bit machines 255 POMEGAA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAUA(i,1,k)) 256 POMEGAA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i,2,k)) 257 ENDDO 258 ENDDO 259 ! 260 DO k = 1, kflev+1 261 DO i = 1, kdlon 262 PPMB(i,k) = paprs(iof+i,k)/100.0 263 ENDDO 264 ENDDO 265 ! 266 DO kk = 1, 5 235 267 DO k = 1, kflev 268 DO i = 1, kdlon 269 PAER(i,k,kk) = 1.0E-15 270 ENDDO 271 ENDDO 272 ENDDO 273 DO k = 1, kflev 236 274 DO i = 1, kdlon 237 PDP(i,k) = paprs(iof+i,k)-paprs(iof+i,k+1) 238 PTAVE(i,k) = t(iof+i,k) 239 PWV(i,k) = MAX (q(iof+i,k), 1.0e-12) 240 PQS(i,k) = PWV(i,k) 241 c wo: cm.atm (epaisseur en cm dans la situation standard) 242 c POZON: kg/kg 243 POZON(i,k) = MAX(wo(iof+i,k),1.0e-12)*RG/46.6968 244 . /(paprs(iof+i,k)-paprs(iof+i,k+1)) 245 . *(paprs(iof+i,1)/101325.0) 246 PCLDLD(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k) 247 PCLDLU(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k) 248 PCLDSW(i,k) = cldfra(iof+i,k) 249 PTAU(i,1,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! 1e-12 serait instable 250 PTAU(i,2,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! pour 32-bit machines 251 POMEGA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAU(i,1,k)) 252 POMEGA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i,2,k)) 253 PCG(i,1,k) = 0.865 254 PCG(i,2,k) = 0.910 255 c-OB 256 cjq Introduced for aerosol indirect forcings. 257 cjq The following values use the cloud optical thickness calculated from 258 cjq present-day aerosol concentrations whereas the quantities without the 259 cjq "A" at the end are for pre-industial (natural-only) aerosol concentrations 260 cjq 261 PTAUA(i,1,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! 1e-12 serait instable 262 PTAUA(i,2,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! pour 32-bit machines 263 POMEGAA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAUA(i,1,k)) 264 POMEGAA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i,2,k)) 265 cjq-end 266 ENDDO 267 ENDDO 268 c 275 tauaero(i,k,:,1)=tau_aero(iof+i,k,:,1) 276 pizaero(i,k,:,1)=piz_aero(iof+i,k,:,1) 277 cgaero(i,k,:,1) =cg_aero(iof+i,k,:,1) 278 tauaero(i,k,:,2)=tau_aero(iof+i,k,:,2) 279 pizaero(i,k,:,2)=piz_aero(iof+i,k,:,2) 280 cgaero(i,k,:,2) =cg_aero(iof+i,k,:,2) 281 ENDDO 282 ENDDO 283 ! 284 !====================================================================== 285 CALL LW_LMDAR4(& 286 PPMB, PDP,& 287 PPSOL,PDT0,PEMIS,& 288 PTL, PTAVE, PWV, POZON, PAER,& 289 PCLDLD,PCLDLU,& 290 PVIEW,& 291 zcool, zcool0,& 292 ztoplw,zsollw,ztoplw0,zsollw0,& 293 zsollwdown,& 294 ZFLUP, ZFLDN, ZFLUP0,ZFLDN0) 295 296 297 IF (.NOT. new_aod) THEN 298 ! use old version 299 CALL SW_LMDAR4(PSCT, zrmu0, zfract,& 300 PPMB, PDP, & 301 PPSOL, PALBD, PALBP,& 302 PTAVE, PWV, PQS, POZON, PAER,& 303 PCLDSW, PTAU, POMEGA, PCG,& 304 zheat, zheat0,& 305 zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,& 306 ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,& 307 tau_aero(:,:,5,:), piz_aero(:,:,5,:), cg_aero(:,:,5,:),& 308 PTAUA, POMEGAA,& 309 ztopswadaero,zsolswadaero,& 310 ztopswaiaero,zsolswaiaero,& 311 ok_ade, ok_aie) 312 ELSE 313 314 CALL SW_AERO(PSCT, zrmu0, zfract,& 315 PPMB, PDP,& 316 PPSOL, PALBD, PALBP,& 317 PTAVE, PWV, PQS, POZON, PAER,& 318 PCLDSW, PTAU, POMEGA, PCG,& 319 zheat, zheat0,& 320 zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,& 321 ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,& 322 tauaero, pizaero, cgaero, & 323 PTAUA, POMEGAA,& 324 ztopswadaero,zsolswadaero,& 325 ztopswad0aero,zsolswad0aero,& 326 ztopswaiaero,zsolswaiaero, & 327 ztopsw_aero,ztopsw0_aero,& 328 zsolsw_aero,zsolsw0_aero,& 329 ok_ade, ok_aie) 330 331 ENDIF 332 !====================================================================== 333 DO i = 1, kdlon 334 radsol(iof+i) = zsolsw(i) + zsollw(i) 335 topsw(iof+i) = ztopsw(i) 336 toplw(iof+i) = ztoplw(i) 337 solsw(iof+i) = zsolsw(i) 338 sollw(iof+i) = zsollw(i) 339 sollwdown(iof+i) = zsollwdown(i) 269 340 DO k = 1, kflev+1 341 lwdn0 ( iof+i,k) = ZFLDN0 ( i,k) 342 lwdn ( iof+i,k) = ZFLDN ( i,k) 343 lwup0 ( iof+i,k) = ZFLUP0 ( i,k) 344 lwup ( iof+i,k) = ZFLUP ( i,k) 345 ENDDO 346 topsw0(iof+i) = ztopsw0(i) 347 toplw0(iof+i) = ztoplw0(i) 348 solsw0(iof+i) = zsolsw0(i) 349 sollw0(iof+i) = zsollw0(i) 350 albpla(iof+i) = zalbpla(i) 351 352 DO k = 1, kflev+1 353 swdn0 ( iof+i,k) = ZFSDN0 ( i,k) 354 swdn ( iof+i,k) = ZFSDN ( i,k) 355 swup0 ( iof+i,k) = ZFSUP0 ( i,k) 356 swup ( iof+i,k) = ZFSUP ( i,k) 357 ENDDO 358 ENDDO 359 !-transform the aerosol forcings, if they have 360 ! to be calculated 361 IF (ok_ade) THEN 362 DO i = 1, kdlon 363 topswad_aero(iof+i) = ztopswadaero(i) 364 topswad0_aero(iof+i) = ztopswad0aero(i) 365 solswad_aero(iof+i) = zsolswadaero(i) 366 solswad0_aero(iof+i) = zsolswad0aero(i) 367 topsw_aero(iof+i,:) = ztopsw_aero(iof+i,:) 368 topsw0_aero(iof+i,:) = ztopsw0_aero(iof+i,:) 369 solsw_aero(iof+i,:) = zsolsw_aero(iof+i,:) 370 solsw0_aero(iof+i,:) = zsolsw0_aero(iof+i,:) 371 372 ENDDO 373 ELSE 374 DO i = 1, kdlon 375 topswad_aero(iof+i) = 0.0 376 solswad_aero(iof+i) = 0.0 377 topswad0_aero(iof+i) = 0.0 378 solswad0_aero(iof+i) = 0.0 379 topsw_aero(iof+i,:) = 0. 380 topsw0_aero(iof+i,:) =0. 381 solsw_aero(iof+i,:) = 0. 382 solsw0_aero(iof+i,:) = 0. 383 ENDDO 384 ENDIF 385 IF (ok_aie) THEN 386 DO i = 1, kdlon 387 topswai_aero(iof+i) = ztopswaiaero(i) 388 solswai_aero(iof+i) = zsolswaiaero(i) 389 ENDDO 390 ELSE 391 DO i = 1, kdlon 392 topswai_aero(iof+i) = 0.0 393 solswai_aero(iof+i) = 0.0 394 ENDDO 395 ENDIF 396 DO k = 1, kflev 270 397 DO i = 1, kdlon 271 PPMB(i,k) = paprs(iof+i,k)/100.0 272 ENDDO 273 ENDDO 274 c 275 DO kk = 1, 5 276 DO k = 1, kflev 277 DO i = 1, kdlon 278 PAER(i,k,kk) = 1.0E-15 279 ENDDO 280 ENDDO 281 ENDDO 282 c-OB 283 DO k = 1, kflev 284 DO i = 1, kdlon 285 tauae(i,k,1)=tau_ae(iof+i,k,1) 286 pizae(i,k,1)=piz_ae(iof+i,k,1) 287 cgae(i,k,1) =cg_ae(iof+i,k,1) 288 tauae(i,k,2)=tau_ae(iof+i,k,2) 289 pizae(i,k,2)=piz_ae(iof+i,k,2) 290 cgae(i,k,2) =cg_ae(iof+i,k,2) 291 ENDDO 292 ENDDO 293 c 294 c===== si iflag_rrtm=0 ================================================ 295 cIM ctes ds clesphys.h CALL LW(RCO2,RCH4,RN2O,RCFC11,RCFC12, 296 cIM ctes ds clesphys.h CALL SW(PSCT, RCO2, zrmu0, zfract, 297 c 298 if (iflag_rrtm.eq.0) then 299 CALL LW_LMDAR4( 300 . PPMB, PDP, 301 . PPSOL,PDT0,PEMIS, 302 . PTL, PTAVE, PWV, POZON, PAER, 303 . PCLDLD,PCLDLU, 304 . PVIEW, 305 . zcool, zcool0, 306 . ztoplw,zsollw,ztoplw0,zsollw0, 307 . zsollwdown, 308 . ZFLUP, ZFLDN, ZFLUP0,ZFLDN0) 309 CALL SW_LMDAR4(PSCT, zrmu0, zfract, 310 S PPMB, PDP, 311 S PPSOL, PALBD, PALBP, 312 S PTAVE, PWV, PQS, POZON, PAER, 313 S PCLDSW, PTAU, POMEGA, PCG, 314 S zheat, zheat0, 315 S zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0, 316 S ZFSUP,ZFSDN,ZFSUP0,ZFSDN0, 317 S tauae, pizae, cgae, ! aerosol optical properties 318 s PTAUA, POMEGAA, 319 s ztopswad,zsolswad,ztopswai,zsolswai, ! diagnosed aerosol forcing 320 J ok_ade, ok_aie) ! apply aerosol effects or not? 321 else 322 c===== si iflag_rrtm=1, on passe dans SW via RECMWFL =============== 323 PRINT*, "Cette option ne fonctionne pas encore !!!" 324 CALL abort 325 endif ! if(iflag_rrtm=0) 326 327 c====================================================================== 328 DO i = 1, kdlon 329 radsol(iof+i) = zsolsw(i) + zsollw(i) 330 topsw(iof+i) = ztopsw(i) 331 toplw(iof+i) = ztoplw(i) 332 solsw(iof+i) = zsolsw(i) 333 sollw(iof+i) = zsollw(i) 334 sollwdown(iof+i) = zsollwdown(i) 335 cIM 336 DO k = 1, kflev+1 337 lwdn0 ( iof+i,k) = ZFLDN0 ( i,k) 338 lwdn ( iof+i,k) = ZFLDN ( i,k) 339 lwup0 ( iof+i,k) = ZFLUP0 ( i,k) 340 lwup ( iof+i,k) = ZFLUP ( i,k) 341 ENDDO 342 c 343 topsw0(iof+i) = ztopsw0(i) 344 toplw0(iof+i) = ztoplw0(i) 345 solsw0(iof+i) = zsolsw0(i) 346 sollw0(iof+i) = zsollw0(i) 347 albpla(iof+i) = zalbpla(i) 348 cIM 349 DO k = 1, kflev+1 350 swdn0 ( iof+i,k) = ZFSDN0 ( i,k) 351 swdn ( iof+i,k) = ZFSDN ( i,k) 352 swup0 ( iof+i,k) = ZFSUP0 ( i,k) 353 swup ( iof+i,k) = ZFSUP ( i,k) 354 ENDDO !k=1, kflev+1 355 ENDDO 356 cjq-transform the aerosol forcings, if they have 357 cjq to be calculated 358 IF (ok_ade) THEN 359 DO i = 1, kdlon 360 topswad(iof+i) = ztopswad(i) 361 solswad(iof+i) = zsolswad(i) 362 ENDDO 363 ELSE 364 DO i = 1, kdlon 365 topswad(iof+i) = 0.0 366 solswad(iof+i) = 0.0 367 ENDDO 368 ENDIF 369 IF (ok_aie) THEN 370 DO i = 1, kdlon 371 topswai(iof+i) = ztopswai(i) 372 solswai(iof+i) = zsolswai(i) 373 ENDDO 374 ELSE 375 DO i = 1, kdlon 376 topswai(iof+i) = 0.0 377 solswai(iof+i) = 0.0 378 ENDDO 379 ENDIF 380 cjq-end 381 DO k = 1, kflev 382 c DO i = 1, kdlon 383 c heat(iof+i,k) = zheat(i,k) 384 c cool(iof+i,k) = zcool(i,k) 385 c heat0(iof+i,k) = zheat0(i,k) 386 c cool0(iof+i,k) = zcool0(i,k) 387 c ENDDO 388 DO i = 1, kdlon 389 C scale factor to take into account the difference between 390 C dry air and watter vapour scpecific heat capacity 391 zznormcp=1.0+RVTMP2*PWV(i,k) 392 heat(iof+i,k) = zheat(i,k)/zznormcp 393 cool(iof+i,k) = zcool(i,k)/zznormcp 394 heat0(iof+i,k) = zheat0(i,k)/zznormcp 395 cool0(iof+i,k) = zcool0(i,k)/zznormcp 396 ENDDO 397 ENDDO 398 c 399 99999 CONTINUE 400 RETURN 401 END 398 ! scale factor to take into account the difference between 399 ! dry air and watter vapour scpecifi! heat capacity 400 zznormcp=1.0+RVTMP2*PWV(i,k) 401 heat(iof+i,k) = zheat(i,k)/zznormcp 402 cool(iof+i,k) = zcool(i,k)/zznormcp 403 heat0(iof+i,k) = zheat0(i,k)/zznormcp 404 cool0(iof+i,k) = zcool0(i,k)/zznormcp 405 ENDDO 406 ENDDO 407 ! 408 ENDDO 409 410 411 ENDSUBROUTINE radlwsw_aero 412 413 -
LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol.F90
r1149 r1150 2 2 ! $Header$ 3 3 ! 4 SUBROUTINE readsulfate (r_day, first, sulfate_p) 5 USE dimphy, ONLY : klev 6 USE mod_grid_phy_lmdz, klon=>klon_glo 7 USE mod_phys_lmdz_para 8 IMPLICIT none 9 10 c Content: 11 c -------- 12 c This routine reads in monthly mean values of sulfate aerosols and 13 c returns a linearly interpolated dayly-mean field. 14 c 15 c 16 c Author: 17 c ------- 18 c Johannes Quaas (quaas@lmd.jussieu.fr) 19 c 26/04/01 20 c 21 c Modifications: 22 c -------------- 23 c 21/06/01: Make integrations of more than one year possible ;-) 24 c ATTENTION!! runs are supposed to start with Jan, 1. 1930 25 c (rday=1) 26 c 27 c 27/06/01: Correction: The model always has 360 days per year! 28 c 27/06/01: SO4 concentration rather than mixing ratio 29 c 27/06/01: 10yr-mean-values to interpolate 30 c 20/08/01: Correct the error through integer-values in interpolations 31 c 21/08/01: Introduce flag to read in just one decade 32 c 33 c Include-files: 34 c -------------- 35 #include "YOMCST.h" 36 #include "chem.h" 37 #include "dimensions.h" 38 #include "temps.h" 39 #include "clesphys.h" 40 #include "iniprint.h" 41 c 42 c Input: 43 c ------ 44 REAL r_day ! Day of integration 45 LOGICAL first ! First timestep 46 ! (and therefore initialization necessary) 47 c 48 c Output: 49 c ------- 50 REAL sulfate_p(klon_omp,klev) 51 REAL sulfate (klon, klev) ! Mass of sulfate (monthly mean data, 52 ! from file) [ug SO4/m3] 53 c 54 c Local Variables: 55 c ---------------- 56 INTEGER i, ig, k, it 57 INTEGER j, iday, ny, iyr, iyr1, iyr2 58 parameter (ny=jjm+1) 59 60 CJLD INTEGER idec1, idec2 ! The two decadal data read ini 61 CHARACTER*4 cyear 62 63 INTEGER im, day1, day2, im2 64 REAL so4_1(iim, jjm+1, klev, 12) 65 REAL so4_2(iim, jjm+1, klev, 12) ! The sulfate distributions 66 67 REAL, allocatable,save :: so4(:, :, :) ! SO4 in right dimension 68 REAL, allocatable,save :: so4_out(:, :) 69 c$OMP THREADPRIVATE(so4,so4_out) 70 71 LOGICAL lnewday 72 LOGICAL lonlyone 73 PARAMETER (lonlyone=.FALSE.) 74 logical,save :: first2=.true. 75 c$OMP THREADPRIVATE(first2) 76 77 c$OMP MASTER 78 if (first2) then 79 80 allocate( so4(klon, klev, 12) ) 81 allocate( so4_out(klon, klev)) 82 first2=.false. 83 84 endif 85 86 if (is_mpi_root) then 87 88 IF (aer_type /= 'actuel ' .AND. aer_type /= 'preind ' .AND. & 89 & aer_type /= 'scenario') THEN 90 WRITE(lunout,*)' *** Warning ***' 91 WRITE(lunout,*)'Option aer_type pour les aerosols = ', & 92 & aer_type 93 WRITE(lunout,*)'Cas non prevu, force a preind' 94 aer_type = 'preind ' 95 ENDIF 96 4 SUBROUTINE readaerosol (id_aero, r_day, first, massvar_p) 5 6 USE dimphy, ONLY : klev 7 USE mod_grid_phy_lmdz, klon=>klon_glo 8 USE mod_phys_lmdz_para 9 10 IMPLICIT NONE 11 12 ! Content: 13 ! -------- 14 ! This routine reads in monthly mean values of massvar aerosols and 15 ! returns a linearly interpolated dayly-mean field. 16 ! 17 ! 18 ! Author: 19 ! ------- 20 ! Johannes Quaas (quaas@lmd.jussieu.fr) 21 ! Celine Deandreis & Anne Cozic 22 ! 19/02/09 23 ! 24 25 ! 26 ! Include-files: 27 ! -------------- 28 INCLUDE "YOMCST.h" 29 INCLUDE "chem.h" 30 INCLUDE "dimensions.h" 31 INCLUDE "temps.h" 32 INCLUDE "clesphys.h" 33 INCLUDE "iniprint.h" 34 ! 35 ! Input: 36 ! ------ 37 INTEGER, INTENT(in) :: id_aero 38 REAL, INTENT(in) :: r_day ! Day of integration 39 LOGICAL, INTENT(in) :: first ! First timestep 40 ! (and therefore initialization necessary) 41 ! 42 ! Output: 43 ! ------- 44 REAL, INTENT(out) :: massvar_p(klon_omp,klev) ! Mass of massvar (monthly mean data,from file) [ug AIBCM/m3] 45 46 ! 47 ! Local Variables: 48 ! ---------------- 49 INTEGER :: i, ig, k, it 50 INTEGER :: j, iday, iyr, iyr1, iyr2 51 INTEGER :: im, day1, day2, im2 52 INTEGER, PARAMETER :: ny=jjm+1 53 CHARACTER(len=4) :: cyear 54 CHARACTER(len=7),DIMENSION(8) :: name_aero 55 REAL :: var_1(iim, jjm+1, klev, 12) 56 REAL, DIMENSION(klon,klev) :: massvar 57 REAL, DIMENSION(iim, jjm+1, klev, 12) :: var_2 ! The massvar distributions 58 REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: var ! VAR in right dimension 59 REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE :: var_out 60 !$OMP THREADPRIVATE(var,var_out) 61 62 LOGICAL :: lnewday 63 LOGICAL, PARAMETER :: lonlyone=.FALSE. 64 LOGICAL,SAVE :: first2=.TRUE. 65 !$OMP THREADPRIVATE(first2) 66 67 ! Variable pour definir a partir d'un numero d'aerosol, son nom 68 name_aero(1) = "SSSSM " 69 name_aero(2) = "ASSSM " 70 name_aero(3) = "ASBCM " 71 name_aero(4) = "ASPOMM " 72 name_aero(5) = "SO4 " 73 name_aero(6) = "CIDUSTM" 74 name_aero(7) = "AIBCM " 75 name_aero(8) = "AIPOMM " 76 77 !$OMP MASTER 78 IF (first2) THEN 79 ALLOCATE( var(klon, klev, 12,8) ) 80 ALLOCATE( var_out(klon, klev,8)) 81 first2=.FALSE. 82 83 IF (aer_type /= 'actuel ' .AND. aer_type /= 'preind ' .AND. & 84 aer_type /= 'scenario') THEN 85 WRITE(lunout,*)' *** Warning ***' 86 WRITE(lunout,*)'Option aer_type non prevu pour les aerosols = ', & 87 aer_type 88 CALL abort_gcm('readaerosol','Error : aer_type parameter not accepted',1) 89 ENDIF 90 ENDIF 91 92 IF (is_mpi_root) THEN 93 97 94 iday = INT(r_day) 98 99 95 ! Get the year of the run 100 96 iyr = iday/360 101 102 97 ! Get the day of the actual year: 103 98 iday = iday-iyr*360 104 105 99 ! 0.02 is about 0.5/24, namly less than half an hour 106 100 lnewday = (r_day-FLOAT(iday).LT.0.02) 107 101 108 !---------------------------------------------109 !All has to be done only, if a new day begins!110 !---------------------------------------------111 102 ! --------------------------------------------- 103 ! All has to be done only, if a new day begins! 104 ! --------------------------------------------- 105 112 106 IF (lnewday.OR.first) THEN 113 114 im = iday/30 +1 ! the actual month 115 ! annee_ref is the initial year (defined in temps.h) 116 iyr = iyr + annee_ref 117 118 ! Do I have to read new data? (Is this the first day of a year?) 119 IF (first.OR.iday.EQ.1.) THEN 120 ! Initialize values 121 DO it=1,12 122 DO k=1,klev 123 DO i=1,klon 124 so4(i,k,it)=0. 125 ENDDO 126 ENDDO 127 ENDDO 128 129 130 131 IF (aer_type == 'actuel ') then 132 cyear='1980' 133 CALL getso4fromfile(cyear, so4_1) 134 ELSE IF (aer_type == 'preind ') THEN 135 cyear='.nat' 136 CALL getso4fromfile(cyear, so4_1) 137 ELSE 138 IF (iyr .lt. 1850) THEN 139 cyear='.nat' 140 WRITE(*,*) 'getso4 iyr=', iyr,' ',cyear 141 CALL getso4fromfile(cyear, so4_1) 142 ELSE IF (iyr .ge. 2100) THEN 143 cyear='2100' 144 WRITE(*,*) 'getso4 iyr=', iyr,' ',cyear 145 CALL getso4fromfile(cyear, so4_1) 146 ELSE 147 148 ! Read in data: 149 ! a) from actual 10-yr-period 150 151 IF (iyr.LT.1900) THEN 152 iyr1 = 1850 153 iyr2 = 1900 154 ELSE IF (iyr.ge.1900.and.iyr.lt.1920) THEN 155 iyr1 = 1900 156 iyr2 = 1920 107 108 im = iday/30 +1 ! the actual month 109 ! annee_ref is the initial year (defined in temps.h) 110 iyr = iyr + annee_ref 111 112 ! Do I have to read new data? (Is this the first day of a year?) 113 IF (first.OR.iday.EQ.1.) THEN 114 ! Initialize values 115 DO it=1,12 116 DO k=1,klev 117 DO i=1,klon 118 var(i,k,it,id_aero)=0. 119 ENDDO 120 ENDDO 121 ENDDO 122 123 124 IF (aer_type == 'actuel ') then 125 126 cyear='1980' 127 CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero)) 128 129 ELSE IF (aer_type == 'preind ') THEN 130 131 cyear='.nat' 132 CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero)) 133 134 ELSE ! aer_type=scenario 135 136 IF (iyr .LT. 1850) THEN 137 cyear='.nat' 138 WRITE(lunout,*) 'get_aero iyr=', iyr,' ',cyear 139 CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero)) 140 141 ELSE IF (iyr .GE. 2100) THEN 142 cyear='2100' 143 WRITE(lunout,*) 'get_aero iyr=', iyr,' ',cyear 144 CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero)) 145 ELSE 146 ! Read in data: 147 ! a) from actual 10-yr-period 148 IF (iyr.LT.1900) THEN 149 iyr1 = 1850 150 iyr2 = 1900 151 ELSE IF (iyr.GE.1900.AND.iyr.LT.1920) THEN 152 iyr1 = 1900 153 iyr2 = 1920 154 ELSE 155 iyr1 = INT(iyr/10)*10 156 iyr2 = INT(1+iyr/10)*10 157 ENDIF 158 159 WRITE(cyear,'(I4)') iyr1 160 WRITE(lunout,*) 'get_aero iyr=', iyr,' ',cyear 161 162 CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero)) 163 164 ENDIF 165 ! If to read two decades: 166 IF (.NOT.lonlyone) THEN 167 168 ! b) from the next following one 169 WRITE(cyear,'(I4)') iyr2 170 WRITE(lunout,*) 'get_aero iyr=', iyr,' ',cyear 171 172 CALL get_aero_fromfile(cyear, var_2, name_aero(id_aero)) 173 174 ! Interpolate linarily to the actual year: 175 DO it=1,12 176 DO k=1,klev 177 DO j=1,jjm 178 DO i=1,iim 179 var_1(i,j,k,it) = & 180 var_1(i,j,k,it) - FLOAT(iyr-iyr1)/FLOAT(iyr2-iyr1) * & 181 (var_1(i,j,k,it) - var_2(i,j,k,it)) 182 ENDDO 183 ENDDO 184 ENDDO 185 ENDDO 186 187 ENDIF !lonlyone 188 ENDIF ! aer_type 189 190 ! Transform the horizontal 2D-field into the physics-field 191 ! (Also the levels and the latitudes have to be inversed) 192 193 DO it=1,12 194 DO k=1, klev 195 ! a) at the poles, use the zonal mean: 196 DO i=1,iim 197 ! North pole 198 var(1,k,it,id_aero) = & 199 var(1,k,it,id_aero)+ var_1(i,jjm+1,klev+1-k,it) 200 ! South pole 201 var(klon,k,it,id_aero)= & 202 var(klon,k,it,id_aero)+ var_1(i,1,klev+1-k,it) 203 ENDDO 204 var(1,k,it,id_aero) = var(1,k,it,id_aero)/FLOAT(iim) 205 var(klon,k,it,id_aero)= var(klon,k,it,id_aero)/FLOAT(iim) 206 207 ! b) the values between the poles: 208 ig=1 209 DO j=2,jjm 210 DO i=1,iim 211 ig=ig+1 212 213 IF (ig.GT.klon) THEN 214 WRITE(lunout,*) 'Attention dans readaerosol', & 215 name_aero(id_aero), 'ig > klon', ig, klon 216 ENDIF 217 218 var(ig,k,it,id_aero) = var_1(i,jjm+1+1-j,klev+1-k,it) 219 220 ENDDO 221 ENDDO 222 IF (ig.NE.klon-1) THEN 223 PRINT *,'Error in readaerosol', name_aero(id_aero), 'conversion' 224 CALL abort_gcm('readaerosol','Error in readaerosol 1',1) 225 ENDIF 226 ENDDO ! Loop over k (vertical) 227 ENDDO ! Loop over it (months) 228 229 ENDIF ! Had to read new data? 230 231 232 ! Interpolate to actual day: 233 IF (iday.LT.im*30-15) THEN 234 ! in the first half of the month use month before and actual month 235 im2=im-1 236 day2 = im2*30-15 237 day1 = im2*30+15 238 IF (im2.LE.0) THEN 239 ! the month is january, thus the month before december 240 im2=12 241 ENDIF 242 DO k=1,klev 243 DO i=1,klon 244 massvar(i,k) = & 245 var(i,k,im2,id_aero)- FLOAT(iday-day2)/FLOAT(day1-day2) * & 246 (var(i,k,im2,id_aero) - var(i,k,im,id_aero)) 247 248 IF (massvar(i,k).LT.0.) THEN 249 IF (iday-day2.LT.0.) WRITE(lunout,*) 'iday-day2',iday-day2 250 IF (var(i,k,im2,id_aero) - var(i,k,im,id_aero).LT.0.) THEN 251 PRINT*, name_aero(id_aero),'(i,k,im2)- ', & 252 name_aero(id_aero),'(i,k,im)', & 253 var(i,k,im2,id_aero) - var(i,k,im,id_aero) 254 ENDIF 255 256 IF (day1-day2.LT.0.) WRITE(lunout,*) 'day1-day2',day1-day2 257 WRITE(lunout,*) 'stop',name_aero(id_aero) 258 CALL abort_gcm('readaerosol','Error in interpolation 2',1) 259 ENDIF 260 ENDDO 261 ENDDO 157 262 ELSE 158 iyr1 = INT(iyr/10)*10 159 iyr2 = INT(1+iyr/10)*10 263 ! the second half of the month 264 im2=im+1 265 IF (im2.GT.12) THEN 266 ! the month is december, the following thus january 267 im2=1 268 ENDIF 269 day2 = im*30-15 270 day1 = im*30+15 271 DO k=1,klev 272 DO i=1,klon 273 massvar(i,k) = & 274 var(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * & 275 (var(i,k,im2,id_aero) - var(i,k,im,id_aero)) 276 277 IF (massvar(i,k).LT.0.) THEN 278 IF (iday-day2.LT.0.) & 279 WRITE(lunout,*) 'iday-day2',iday-day2 280 IF (var(i,k,im2,id_aero) - var(i,k,im,id_aero).LT.0.) THEN 281 WRITE(lunout,*) name_aero(id_aero),'(i,k,im2)-', & 282 name_aero(id_aero),'(i,k,im)', & 283 var(i,k,im2,id_aero) - var(i,k,im,id_aero) 284 ENDIF 285 IF (day1-day2.LT.0.) & 286 WRITE(lunout,*) 'day1-day2',day1-day2 287 WRITE(lunout,*) 'stop',name_aero(id_aero) 288 CALL abort_gcm('readaerosol','Error in interpolation 3',1) 289 ENDIF 290 ENDDO 291 ENDDO 160 292 ENDIF 161 WRITE(cyear,'(I4)') iyr1 162 ENDIF 163 WRITE(*,*) 'getso4 iyr=', iyr,' ',cyear 164 CALL getso4fromfile(cyear, so4_1) 165 166 167 ! If to read two decades: 168 IF (.NOT.lonlyone) THEN 169 170 ! b) from the next following one 171 WRITE(cyear,'(I4)') iyr2 172 WRITE(*,*) 'getso4 iyr=', iyr,' ',cyear 173 CALL getso4fromfile(cyear, so4_2) 174 293 294 295 ! The massvar concentration [molec cm-3] is read in. 296 ! Convert it into mass [ug VAR/m3] 297 ! masse_var in [g/mol], n_avogadro in [molec/mol] 298 ! The sulfate mass [ug VAR/m3] is read in. 299 DO k=1,klev 300 DO i=1,klon 301 var_out(i,k,id_aero) = massvar(i,k) 302 IF (var_out(i,k,id_aero).LT.0) THEN 303 PRINT*, 'Attention il y a un probleme dans readaerosol', & 304 name_aero(id_aero), 'la masse est negative' 305 CALL abort_gcm('readaerosol','Error in readaerosol 4',1) 306 ENDIF 307 ENDDO 308 ENDDO 309 ELSE ! if no new day, use old data: 310 DO k=1,klev 311 DO i=1,klon 312 massvar(i,k) = var_out(i,k,id_aero) 313 IF (var_out(i,k,id_aero).LT.0) THEN 314 PRINT*, 'Attention il y a un probleme dans readaerosol', & 315 name_aero(id_aero), 'la masse est negative' 316 CALL abort_gcm('readaerosol','Error in readaerosol 5',1) 317 ENDIF 318 ENDDO 319 ENDDO 320 321 322 ENDIF ! Did I have to do anything (was it a new day?) 323 324 ENDIF ! phy_rank==0 325 326 !$OMP END MASTER 327 CALL Scatter(massvar,massvar_p) 328 329 END SUBROUTINE readaerosol 330 331 332 333 334 335 !----------------------------------------------------------------------------- 336 ! Read in /calculate pre-industrial values of sulfate 337 !----------------------------------------------------------------------------- 338 339 SUBROUTINE readaerosol_preind (id_aero, r_day, first, pi_massvar_p) 340 341 USE dimphy, ONLY : klev 342 USE mod_grid_phy_lmdz, klon=>klon_glo 343 USE mod_phys_lmdz_para 344 IMPLICIT NONE 345 346 ! Content: 347 ! -------- 348 ! This routine reads in monthly mean values of massvar aerosols and 349 ! returns a linearly interpolated dayly-mean field. 350 ! 351 ! It does so for the preindustriel values of the massvar, to a large part 352 ! analogous to the routine readmassvar above. 353 ! 354 ! Only Pb: Variables must be saved and don t have to be overwritten! 355 ! 356 ! Author: 357 ! ------- 358 ! Celine Deandreis & Anne Cozic (LSCE) 2009 359 ! Johannes Quaas (quaas@lmd.jussieu.fr) 360 ! 26/06/01 361 ! 362 ! Modifications: 363 ! -------------- 364 ! see above 365 366 367 ! Include-files: 368 ! -------------- 369 370 INCLUDE "YOMCST.h" 371 INCLUDE "chem.h" 372 INCLUDE "dimensions.h" 373 INCLUDE "temps.h" 374 INCLUDE "iniprint.h" 175 375 176 ! Interpolate linarily to the actual year: 177 DO it=1,12 178 DO k=1,klev 179 DO j=1,jjm 180 DO i=1,iim 181 so4_1(i,j,k,it)=so4_1(i,j,k,it) 182 . - FLOAT(iyr-iyr1)/FLOAT(iyr2-iyr1) 183 . * (so4_1(i,j,k,it) - so4_2(i,j,k,it)) 184 ENDDO 185 ENDDO 186 ENDDO 187 ENDDO 188 189 190 ENDIF !lonlyone 191 ENDIF !aer_type 192 193 ! Transform the horizontal 2D-field into the physics-field 194 ! (Also the levels and the latitudes have to be inversed) 195 196 DO it=1,12 197 DO k=1, klev 198 ! a) at the poles, use the zonal mean: 199 DO i=1,iim 200 ! North pole 201 so4(1,k,it)=so4(1,k,it)+so4_1(i,jjm+1,klev+1-k,it) 202 ! South pole 203 so4(klon,k,it)=so4(klon,k,it)+so4_1(i,1,klev+1-k,it) 204 ENDDO 205 so4(1,k,it)=so4(1,k,it)/FLOAT(iim) 206 so4(klon,k,it)=so4(klon,k,it)/FLOAT(iim) 207 208 ! b) the values between the poles: 209 ig=1 210 DO j=2,jjm 211 DO i=1,iim 212 ig=ig+1 213 if (ig.gt.klon) write (*,*) 'shit' 214 so4(ig,k,it) = so4_1(i,jjm+1+1-j,klev+1-k,it) 215 ENDDO 216 ENDDO 217 IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)' 218 ENDDO ! Loop over k (vertical) 219 ENDDO ! Loop over it (months) 220 221 222 ENDIF ! Had to read new data? 223 224 225 ! Interpolate to actual day: 226 IF (iday.LT.im*30-15) THEN 227 ! in the first half of the month use month before and actual month 228 im2=im-1 229 day2 = im2*30-15 230 day1 = im2*30+15 231 IF (im2.LE.0) THEN 232 ! the month is january, thus the month before december 233 im2=12 234 ENDIF 235 DO k=1,klev 236 DO i=1,klon 237 sulfate(i,k) = so4(i,k,im2) 238 . - FLOAT(iday-day2)/FLOAT(day1-day2) 239 . * (so4(i,k,im2) - so4(i,k,im)) 240 IF (sulfate(i,k).LT.0.) THEN 241 IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2 242 IF (so4(i,k,im2) - so4(i,k,im).LT.0.) 243 . write(*,*) 'so4(i,k,im2) - so4(i,k,im)', 244 . so4(i,k,im2) - so4(i,k,im) 245 IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2 246 stop 'sulfate' 247 endif 248 ENDDO 249 ENDDO 250 ELSE 251 ! the second half of the month 252 im2=im+1 253 IF (im2.GT.12) THEN 254 ! the month is december, the following thus january 255 im2=1 256 ENDIF 257 day2 = im*30-15 258 day1 = im*30+15 259 DO k=1,klev 260 DO i=1,klon 261 sulfate(i,k) = so4(i,k,im2) 262 . - FLOAT(iday-day2)/FLOAT(day1-day2) 263 . * (so4(i,k,im2) - so4(i,k,im)) 264 IF (sulfate(i,k).LT.0.) THEN 265 IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2 266 IF (so4(i,k,im2) - so4(i,k,im).LT.0.) 267 . write(*,*) 'so4(i,k,im2) - so4(i,k,im)', 268 . so4(i,k,im2) - so4(i,k,im) 269 IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2 270 stop 'sulfate' 271 endif 272 ENDDO 273 ENDDO 274 ENDIF 275 276 277 CJLD ! The sulfate concentration [molec cm-3] is read in. 278 CJLD ! Convert it into mass [ug SO4/m3] 279 CJLD ! masse_so4 in [g/mol], n_avogadro in [molec/mol] 280 ! The sulfate mass [ug SO4/m3] is read in. 281 DO k=1,klev 282 DO i=1,klon 283 CJLD sulfate(i,k) = sulfate(i,k)*masse_so4 284 CJLD . /n_avogadro*1.e+12 285 so4_out(i,k) = sulfate(i,k) 286 IF (so4_out(i,k).LT.0) 287 . stop 'WAS SOLL DER SCHEISS ? ' 288 ENDDO 289 ENDDO 290 ELSE ! if no new day, use old data: 291 DO k=1,klev 292 DO i=1,klon 293 sulfate(i,k) = so4_out(i,k) 294 IF (so4_out(i,k).LT.0) 295 . stop 'WAS SOLL DER SCHEISS ? ' 296 ENDDO 297 ENDDO 298 299 300 ENDIF ! Did I have to do anything (was it a new day?) 301 302 endif ! phy_rank==0 303 304 c$OMP END MASTER 305 call Scatter(sulfate,sulfate_p) 306 307 RETURN 308 END 309 310 311 312 313 314 c----------------------------------------------------------------------------- 315 c Read in /calculate pre-industrial values of sulfate 316 c----------------------------------------------------------------------------- 317 318 SUBROUTINE readsulfate_preind (r_day, first, pi_sulfate_p) 319 USE dimphy, ONLY : klev 320 USE mod_grid_phy_lmdz, klon=>klon_glo 321 USE mod_phys_lmdz_para 322 IMPLICIT none 323 324 c Content: 325 c -------- 326 c This routine reads in monthly mean values of sulfate aerosols and 327 c returns a linearly interpolated dayly-mean field. 328 c 329 c It does so for the preindustriel values of the sulfate, to a large part 330 c analogous to the routine readsulfate above. 331 c 332 c Only Pb: Variables must be saved and don t have to be overwritten! 333 c 334 c Author: 335 c ------- 336 c Johannes Quaas (quaas@lmd.jussieu.fr) 337 c 26/06/01 338 c 339 c Modifications: 340 c -------------- 341 c see above 342 c 343 c Include-files: 344 c -------------- 345 #include "YOMCST.h" 346 #include "chem.h" 347 #include "dimensions.h" 348 #include "temps.h" 349 c 350 c Input: 351 c ------ 352 REAL r_day ! Day of integration 353 LOGICAL first ! First timestep 354 ! (and therefore initialization necessary) 355 c 356 c Output: 357 c ------- 358 REAL pi_sulfate_p (klon_omp, klev) 359 360 REAL pi_sulfate (klon, klev) ! Number conc. sulfate (monthly mean data, 361 ! from fil 362 c 363 c Local Variables: 364 c ---------------- 365 INTEGER i, ig, k, it 366 INTEGER j, iday, ny, iyr 367 parameter (ny=jjm+1) 368 369 INTEGER im, day1, day2, im2 370 REAL pi_so4_1(iim, jjm+1, klev, 12) 371 372 REAL, allocatable,save :: pi_so4(:, :, :) ! SO4 in right dimension 373 REAL, allocatable,save :: pi_so4_out(:, :) 374 c$OMP THREADPRIVATE(pi_so4,pi_so4_out) 375 376 CHARACTER*4 cyear 377 LOGICAL lnewday 378 logical,save :: first2=.true. 379 c$OMP THREADPRIVATE(first2) 380 381 c$OMP MASTER 382 if (first2) then 383 384 allocate( pi_so4(klon, klev, 12) ) 385 allocate( pi_so4_out(klon, klev)) 386 first2=.false. 387 388 endif 389 390 if (is_mpi_root) then 391 392 393 376 ! Input: 377 ! ------ 378 INTEGER, INTENT(in) :: id_aero 379 REAL, INTENT(in) :: r_day ! Day of integration 380 LOGICAL, INTENT(in) :: first ! First timestep (and therefore initialization necessary) 381 382 383 ! 384 ! Output: 385 ! ------- 386 REAL, DIMENSION(klon_omp,klev), INTENT(out) :: pi_massvar_p ! Number conc. massvar (monthly mean data, from file) 387 388 389 ! 390 ! Local Variables: 391 ! ---------------- 392 INTEGER :: i, ig, k, it 393 INTEGER :: j, iday, iyr 394 INTEGER, PARAMETER :: ny=jjm+1 395 INTEGER :: im, day1, day2, im2 396 397 REAL, DIMENSION(iim, jjm+1, klev, 12) :: pi_var_1 398 REAL, DIMENSION(klon,klev) :: pi_massvar 399 REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: pi_var ! VAR in right dimension 400 REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: pi_var_out 401 !$OMP THREADPRIVATE(pi_var,pi_var_out) 402 403 CHARACTER(len=4) :: cyear 404 CHARACTER(len=7),DIMENSION(8) :: name_aero 405 LOGICAL :: lnewday 406 LOGICAL,SAVE :: first2=.TRUE. 407 !$OMP THREADPRIVATE(first2) 408 409 410 ! Variable pour definir a partir d'un numero d'aerosol, son nom 411 name_aero(1) = "SSSSM " 412 name_aero(2) = "ASSSM " 413 name_aero(3) = "ASBCM " 414 name_aero(4) = "ASPOMM " 415 name_aero(5) = "SO4 " 416 name_aero(6) = "CIDUSTM" 417 name_aero(7) = "AIBCM " 418 name_aero(8) = "AIPOMM " 419 420 421 !$OMP MASTER 422 IF (first2) THEN 423 ALLOCATE( pi_var(klon, klev, 12,8) ) 424 ALLOCATE( pi_var_out(klon, klev,8)) 425 first2=.FALSE. 426 ENDIF 427 428 IF (is_mpi_root) THEN 394 429 iday = INT(r_day) 395 396 430 ! Get the year of the run 397 431 iyr = iday/360 398 399 432 ! Get the day of the actual year: 400 433 iday = iday-iyr*360 401 402 434 ! 0.02 is about 0.5/24, namly less than half an hour 403 435 lnewday = (r_day-FLOAT(iday).LT.0.02) 404 436 405 !---------------------------------------------406 !All has to be done only, if a new day begins!407 !---------------------------------------------408 437 ! --------------------------------------------- 438 ! All has to be done only, if a new day begins! 439 ! --------------------------------------------- 440 409 441 IF (lnewday.OR.first) THEN 442 443 im = iday/30 +1 ! the actual month 444 ! annee_ref is the initial year (defined in temps.h) 445 iyr = iyr + annee_ref 446 447 IF (first) THEN 448 cyear='.nat' 449 CALL get_aero_fromfile(cyear,pi_var_1, name_aero(id_aero)) 450 451 ! Transform the horizontal 2D-field into the physics-field 452 ! (Also the levels and the latitudes have to be inversed) 453 454 ! Initialize field 455 DO it=1,12 456 DO k=1,klev 457 DO i=1,klon 458 pi_var(i,k,it,id_aero)=0. 459 ENDDO 460 ENDDO 461 ENDDO 462 463 WRITE(lunout,*) 'preind: finished reading', FLOAT(iim) 464 DO it=1,12 465 DO k=1, klev 466 ! a) at the poles, use the zonal mean: 467 DO i=1,iim 468 ! North pole 469 pi_var(1,k,it,id_aero) = & 470 pi_var(1,k,it,id_aero) + pi_var_1(i,jjm+1,klev+1-k,it) 471 ! South pole 472 pi_var(klon,k,it,id_aero) = & 473 pi_var(klon,k,it,id_aero) + pi_var_1(i,1,klev+1-k,it) 474 ENDDO 475 pi_var(1,k,it,id_aero) = pi_var(1,k,it,id_aero)/FLOAT(iim) 476 pi_var(klon,k,it,id_aero) = pi_var(klon,k,it,id_aero)/FLOAT(iim) 477 478 ! b) the values between the poles: 479 ig=1 480 DO j=2,jjm 481 DO i=1,iim 482 ig=ig+1 483 IF (ig.GT.klon) THEN 484 WRITE(lunout,*) 'Attention dans readaerosol_preind', & 485 name_aero(id_aero), 'ig > klon', ig, klon 486 ENDIF 487 pi_var(ig,k,it,id_aero) = & 488 pi_var_1(i,jjm+1+1-j,klev+1-k,it) 489 ENDDO 490 ENDDO 491 IF (ig.NE.klon-1) THEN 492 WRITE(lunout,*) 'Error in readaerosol_preind (', name_aero(id_aero), ' conversion)' 493 CALL abort_gcm('readaerosol_preind','Error in readaerosol_preind 1',1) 494 ENDIF 495 ENDDO ! Loop over k (vertical) 496 ENDDO ! Loop over it (months) 497 498 ENDIF ! Had to read new data? 499 500 501 ! Interpolate to actual day: 502 IF (iday.LT.im*30-15) THEN 503 ! in the first half of the month use month before and actual month 504 im2=im-1 505 day1 = im2*30+15 506 day2 = im2*30-15 507 IF (im2.LE.0) THEN 508 ! the month is january, thus the month before december 509 im2=12 510 ENDIF 511 DO k=1,klev 512 DO i=1,klon 513 pi_massvar(i,k) = & 514 pi_var(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * & 515 (pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero)) 516 517 IF (pi_massvar(i,k).LT.0.) THEN 518 IF (iday-day2.LT.0.) WRITE(lunout,*) 'iday-day2',iday-day2 519 IF (pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero).LT.0.) THEN 520 WRITE(lunout,*) 'pi_',name_aero(id_aero),'(i,k,im2) - pi_', & 521 name_aero(id_aero),'(i,k,im)', & 522 pi_var(i,k,im2,id_aero) -pi_var(i,k,im,id_aero) 523 ENDIF 524 IF (day1-day2.LT.0.)WRITE(lunout,*) 'day1-day2',day1-day2 525 PRINT *, 'pi_',name_aero(id_aero) 526 CALL abort_gcm('readaerosol_preind','Error in readaerosol_preind 2',1) 527 ENDIF 528 ENDDO 529 ENDDO 530 ELSE 531 532 ! the second half of the month 533 im2=im+1 534 day1 = im*30+15 535 IF (im2.GT.12) THEN 536 ! the month is december, the following thus january 537 im2=1 538 ENDIF 539 day2 = im*30-15 540 541 DO k=1,klev 542 DO i=1,klon 543 pi_massvar(i,k) = & 544 pi_var(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * & 545 (pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero)) 546 IF (pi_massvar(i,k).LT.0.) THEN 547 IF (iday-day2.LT.0.) WRITE(lunout,*) 'iday-day2',iday-day2 548 IF (pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero).LT.0.) THEN 549 WRITE(lunout,*) 'pi_', name_aero(id_aero),'(i,k,im2) - pi_',& 550 name_aero(id_aero), '(i,k,im)',& 551 pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero) 552 ENDIF 553 IF (day1-day2.LT.0.) & 554 WRITE(lunout,*) 'day1-day2',day1-day2 555 PRINT *, 'pi_',name_aero(id_aero) 556 CALL abort_gcm('readaerosol_preind','Error in readaerosol_preind 3',1) 557 ENDIF 558 ENDDO 559 ENDDO 560 ENDIF 561 410 562 411 im = iday/30 +1 ! the actual month 412 413 ! annee_ref is the initial year (defined in temps.h) 414 iyr = iyr + annee_ref 415 416 417 IF (first) THEN 418 cyear='.nat' 419 CALL getso4fromfile(cyear,pi_so4_1) 420 421 ! Transform the horizontal 2D-field into the physics-field 422 ! (Also the levels and the latitudes have to be inversed) 423 424 ! Initialize field 425 DO it=1,12 426 DO k=1,klev 427 DO i=1,klon 428 pi_so4(i,k,it)=0. 429 ENDDO 563 ! The massvar concentration [molec cm-3] is read in. 564 ! Convert it into mass [ug VAR/m3] 565 ! masse_var in [g/mol], n_avogadro in [molec/mol] 566 DO k=1,klev 567 DO i=1,klon 568 pi_var_out(i,k,id_aero) = pi_massvar(i,k) 430 569 ENDDO 431 ENDDO 570 ENDDO 571 572 ELSE ! If no new day, use old data: 573 DO k=1,klev 574 DO i=1,klon 575 pi_massvar(i,k) = pi_var_out(i,k,id_aero) 576 ENDDO 577 ENDDO 578 ENDIF ! Was this the beginning of a new day? 579 ENDIF ! is_mpi_root 580 581 !$OMP END MASTER 582 CALL Scatter(pi_massvar,pi_massvar_p) 583 584 END SUBROUTINE readaerosol_preind 585 586 587 588 !----------------------------------------------------------------------------- 589 ! Routine for reading VAR data from files 590 !----------------------------------------------------------------------------- 591 592 593 SUBROUTINE get_aero_fromfile (cyr, var, varname) 594 USE dimphy 595 IMPLICIT NONE 596 597 INCLUDE "netcdf.inc" 598 INCLUDE "dimensions.h" 599 INCLUDE "iniprint.h" 600 601 CHARACTER(len=4), INTENT(in) :: cyr 602 REAL, DIMENSION(iim, jjm+1, klev, 12), INTENT(out) :: var 603 CHARACTER(len=*), INTENT(in) :: varname 604 605 606 CHARACTER(len=30) :: fname 607 CHARACTER(len=30) :: cvar 608 INTEGER, DIMENSION(3) :: START, COUNT 609 INTEGER :: STATUS, NCID, VARID 610 INTEGER :: imth, i, j, k 611 INTEGER, PARAMETER :: ny=jjm+1 612 REAL, DIMENSION(iim, ny, klev) :: varmth 613 ! 614 ! 615 ! 616 fname = TRIM(varname)//'.run'//cyr//'.cdf' 617 618 WRITE(lunout,*) 'reading ', TRIM(fname) 619 STATUS = NF_OPEN (TRIM(fname), NF_NOWRITE, NCID) 620 IF (STATUS .NE. NF_NOERR) THEN 621 WRITE(lunout,*) 'err in open ',status 622 CALL abort_gcm('get_aero_fromfile','Error in opening file',1) 623 ENDIF 624 DO imth=1, 12 625 IF (imth.EQ.1) THEN 626 cvar=TRIM(varname)//'JAN' 627 ELSEIF (imth.EQ.2) THEN 628 cvar=TRIM(varname)//'FEB' 629 ELSEIF (imth.EQ.3) THEN 630 cvar=TRIM(varname)//'MAR' 631 ELSEIF (imth.EQ.4) THEN 632 cvar=TRIM(varname)//'APR' 633 ELSEIF (imth.EQ.5) THEN 634 cvar=TRIM(varname)//'MAY' 635 ELSEIF (imth.EQ.6) THEN 636 cvar=TRIM(varname)//'JUN' 637 ELSEIF (imth.EQ.7) THEN 638 cvar=TRIM(varname)//'JUL' 639 ELSEIF (imth.EQ.8) THEN 640 cvar=TRIM(varname)//'AUG' 641 ELSEIF (imth.EQ.9) THEN 642 cvar=TRIM(varname)//'SEP' 643 ELSEIF (imth.EQ.10) THEN 644 cvar=TRIM(varname)//'OCT' 645 ELSEIF (imth.EQ.11) THEN 646 cvar=TRIM(varname)//'NOV' 647 ELSEIF (imth.EQ.12) THEN 648 cvar=TRIM(varname)//'DEC' 649 ENDIF 650 start(1)=1 651 start(2)=1 652 start(3)=1 653 COUNT(1)=iim 654 COUNT(2)=ny 655 COUNT(3)=klev 656 STATUS = NF_INQ_VARID (NCID, TRIM(cvar), VARID) 657 WRITE(lunout,*) ncid,imth,TRIM(cvar), varid 658 659 IF (STATUS .NE. NF_NOERR) WRITE(lunout,*) 'err in read ',status 660 661 #ifdef NC_DOUBLE 662 status = NF_GET_VARA_DOUBLE(NCID, VARID, START, COUNT,varmth) 663 #else 664 status = NF_GET_VARA_REAL(NCID, VARID, START, COUNT, varmth) 665 #endif 666 IF (STATUS .NE. NF_NOERR) WRITE(lunout,*) 'err in read data',status 432 667 433 write (*,*) 'preind: finished reading', FLOAT(iim) 434 DO it=1,12 435 DO k=1, klev 436 ! a) at the poles, use the zonal mean: 437 DO i=1,iim 438 ! North pole 439 pi_so4(1,k,it)=pi_so4(1,k,it)+pi_so4_1(i,jjm+1,klev+1-k,it) 440 ! South pole 441 pi_so4(klon,k,it)=pi_so4(klon,k,it)+pi_so4_1(i,1,klev+1-k,it) 442 ENDDO 443 pi_so4(1,k,it)=pi_so4(1,k,it)/FLOAT(iim) 444 pi_so4(klon,k,it)=pi_so4(klon,k,it)/FLOAT(iim) 445 446 ! b) the values between the poles: 447 ig=1 448 DO j=2,jjm 449 DO i=1,iim 450 ig=ig+1 451 if (ig.gt.klon) write (*,*) 'shit' 452 pi_so4(ig,k,it) = pi_so4_1(i,jjm+1+1-j,klev+1-k,it) 453 ENDDO 454 ENDDO 455 IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)' 456 ENDDO ! Loop over k (vertical) 457 ENDDO ! Loop over it (months) 458 459 ENDIF ! Had to read new data? 460 461 462 ! Interpolate to actual day: 463 IF (iday.LT.im*30-15) THEN 464 ! in the first half of the month use month before and actual month 465 im2=im-1 466 day1 = im2*30+15 467 day2 = im2*30-15 468 IF (im2.LE.0) THEN 469 ! the month is january, thus the month before december 470 im2=12 471 ENDIF 472 DO k=1,klev 473 DO i=1,klon 474 pi_sulfate(i,k) = pi_so4(i,k,im2) 475 . - FLOAT(iday-day2)/FLOAT(day1-day2) 476 . * (pi_so4(i,k,im2) - pi_so4(i,k,im)) 477 IF (pi_sulfate(i,k).LT.0.) THEN 478 IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2 479 IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.) 480 . write(*,*) 'pi_so4(i,k,im2) - pi_so4(i,k,im)', 481 . pi_so4(i,k,im2) - pi_so4(i,k,im) 482 IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2 483 stop 'pi_sulfate' 484 endif 485 ENDDO 486 ENDDO 487 ELSE 488 ! the second half of the month 489 im2=im+1 490 day1 = im*30+15 491 IF (im2.GT.12) THEN 492 ! the month is december, the following thus january 493 im2=1 494 ENDIF 495 day2 = im*30-15 496 497 DO k=1,klev 498 DO i=1,klon 499 pi_sulfate(i,k) = pi_so4(i,k,im2) 500 . - FLOAT(iday-day2)/FLOAT(day1-day2) 501 . * (pi_so4(i,k,im2) - pi_so4(i,k,im)) 502 IF (pi_sulfate(i,k).LT.0.) THEN 503 IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2 504 IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.) 505 . write(*,*) 'pi_so4(i,k,im2) - pi_so4(i,k,im)', 506 . pi_so4(i,k,im2) - pi_so4(i,k,im) 507 IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2 508 stop 'pi_sulfate' 509 endif 510 ENDDO 511 ENDDO 512 ENDIF 513 514 515 CJLD ! The sulfate concentration [molec cm-3] is read in. 516 CJLD ! Convert it into mass [ug SO4/m3] 517 CJLD ! masse_so4 in [g/mol], n_avogadro in [molec/mol] 518 DO k=1,klev 519 DO i=1,klon 520 CJLD pi_sulfate(i,k) = pi_sulfate(i,k)*masse_so4 521 CJLD . /n_avogadro*1.e+12 522 pi_so4_out(i,k) = pi_sulfate(i,k) 523 ENDDO 668 DO k=1,klev 669 DO j=1,jjm+1 670 DO i=1,iim 671 IF (varmth(i,j,k).LT.0.) THEN 672 WRITE(lunout,*) 'this is shit' 673 WRITE(lunout,*) varname,'(',i,j,k,') =',varmth(i,j,k) 674 ENDIF 675 var(i,j,k,imth)=varmth(i,j,k) 676 ENDDO 524 677 ENDDO 525 526 ELSE ! If no new day, use old data: 527 DO k=1,klev 528 DO i=1,klon 529 pi_sulfate(i,k) = pi_so4_out(i,k) 530 ENDDO 531 ENDDO 532 533 534 ENDIF ! Was this the beginning of a new day? 535 536 endif ! is_mpi_root==0 537 538 c$OMP END MASTER 539 call Scatter(pi_sulfate,pi_sulfate_p) 540 541 RETURN 542 END 543 544 545 546 547 548 549 550 551 552 553 c----------------------------------------------------------------------------- 554 c Routine for reading SO4 data from files 555 c----------------------------------------------------------------------------- 556 557 558 SUBROUTINE getso4fromfile (cyr, so4) 559 USE dimphy 560 #include "netcdf.inc" 561 #include "dimensions.h" 562 CHARACTER*15 fname 563 CHARACTER*4 cyr 564 565 CHARACTER*6 cvar 566 INTEGER START(3), COUNT(3) 567 INTEGER STATUS, NCID, VARID 568 INTEGER imth, i, j, k, ny 569 PARAMETER (ny=jjm+1) 570 571 572 REAL so4mth(iim, ny, klev) 573 REAL so4(iim, ny, klev, 12) 574 575 576 fname = 'so4.run'//cyr//'.cdf' 577 578 write (*,*) 'reading ', fname 579 STATUS = NF_OPEN (fname, NF_NOWRITE, NCID) 580 IF (STATUS .NE. NF_NOERR) write (*,*) 'err in open ',status 581 582 DO imth=1, 12 583 IF (imth.eq.1) THEN 584 cvar='SO4JAN' 585 ELSEIF (imth.eq.2) THEN 586 cvar='SO4FEB' 587 ELSEIF (imth.eq.3) THEN 588 cvar='SO4MAR' 589 ELSEIF (imth.eq.4) THEN 590 cvar='SO4APR' 591 ELSEIF (imth.eq.5) THEN 592 cvar='SO4MAY' 593 ELSEIF (imth.eq.6) THEN 594 cvar='SO4JUN' 595 ELSEIF (imth.eq.7) THEN 596 cvar='SO4JUL' 597 ELSEIF (imth.eq.8) THEN 598 cvar='SO4AUG' 599 ELSEIF (imth.eq.9) THEN 600 cvar='SO4SEP' 601 ELSEIF (imth.eq.10) THEN 602 cvar='SO4OCT' 603 ELSEIF (imth.eq.11) THEN 604 cvar='SO4NOV' 605 ELSEIF (imth.eq.12) THEN 606 cvar='SO4DEC' 607 ENDIF 608 start(1)=1 609 start(2)=1 610 start(3)=1 611 count(1)=iim 612 count(2)=ny 613 count(3)=klev 614 c write(*,*) 'here i am' 615 STATUS = NF_INQ_VARID (NCID, cvar, VARID) 616 write (*,*) ncid,imth,cvar, varid 617 618 IF (STATUS .NE. NF_NOERR) write (*,*) 'err in read ',status 619 620 #ifdef NC_DOUBLE 621 status = NF_GET_VARA_DOUBLE(NCID, VARID, START, COUNT, so4mth) 622 #else 623 status = NF_GET_VARA_REAL(NCID, VARID, START, COUNT, so4mth) 624 #endif 625 IF (STATUS .NE. NF_NOERR) write (*,*) 'err in read data',status 626 627 DO k=1,klev 628 DO j=1,jjm+1 629 DO i=1,iim 630 IF (so4mth(i,j,k).LT.0.) then 631 write(*,*) 'this is shit' 632 write(*,*) 'so4(',i,j,k,') =',so4mth(i,j,k) 633 endif 634 so4(i,j,k,imth)=so4mth(i,j,k) 635 ENDDO 636 ENDDO 637 ENDDO 638 ENDDO 639 640 STATUS = NF_CLOSE(NCID) 641 IF (STATUS .NE. NF_NOERR) write (*,*) 'err in closing file',status 642 643 644 END ! subroutine getso4fromfile 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 678 ENDDO 679 ENDDO 680 681 STATUS = NF_CLOSE(NCID) 682 IF (STATUS .NE. NF_NOERR) WRITE(lunout,*) 'err in closing file',status 683 684 685 END SUBROUTINE get_aero_fromfile 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701
Note: See TracChangeset
for help on using the changeset viewer.
