| 1 | SUBROUTINE SWR ( KDLON, KFLEV, KNU |
|---|
| 2 | S , aerosol,albedo,PDSIG,PPSOL,PRMU,PSEC |
|---|
| 3 | S , PFD,PFU ) |
|---|
| 4 | |
|---|
| 5 | IMPLICIT NONE |
|---|
| 6 | C |
|---|
| 7 | #include "dimensions.h" |
|---|
| 8 | #include "dimphys.h" |
|---|
| 9 | #include "dimradmars.h" |
|---|
| 10 | #include "callkeys.h" |
|---|
| 11 | |
|---|
| 12 | #include "yomaer.h" |
|---|
| 13 | #include "yomlw.h" |
|---|
| 14 | |
|---|
| 15 | #include "fisice.h" |
|---|
| 16 | |
|---|
| 17 | #include "aerice.h" |
|---|
| 18 | C |
|---|
| 19 | C SWR - Continuum scattering computations |
|---|
| 20 | C |
|---|
| 21 | C PURPOSE. |
|---|
| 22 | C -------- |
|---|
| 23 | C Computes the reflectivity and transmissivity in case oF |
|---|
| 24 | C Continuum scattering |
|---|
| 25 | c F. Forget (1999) |
|---|
| 26 | c |
|---|
| 27 | c BASED ON MORCRETTE EARTH MODEL |
|---|
| 28 | C (See radiation's part of the ecmwf research department |
|---|
| 29 | C documentation, and Fouquart and BonneL (1980) |
|---|
| 30 | C |
|---|
| 31 | C IMPLICIT ARGUMENTS : |
|---|
| 32 | C -------------------- |
|---|
| 33 | C |
|---|
| 34 | C ==== INPUTS === |
|---|
| 35 | c |
|---|
| 36 | c KDLON : number of horizontal grid points |
|---|
| 37 | c KFLEV : number of vertical layers |
|---|
| 38 | c KNU : Solar band # (1 or 2) |
|---|
| 39 | c aerosol aerosol extinction optical depth |
|---|
| 40 | c at reference wavelength "longrefvis" set |
|---|
| 41 | c in dimradmars.h , in each layer, for one of |
|---|
| 42 | c the "naerkind" kind of aerosol optical properties. |
|---|
| 43 | c albedo hemispheric surface albedo |
|---|
| 44 | c albedo (i,1) : mean albedo for solar band#1 |
|---|
| 45 | c (see below) |
|---|
| 46 | c albedo (i,2) : mean albedo for solar band#2 |
|---|
| 47 | c (see below) |
|---|
| 48 | c PDSIG layer thickness in sigma coordinates |
|---|
| 49 | c PPSOL Surface pressure (Pa) |
|---|
| 50 | c PRMU: cos of solar zenith angle (=1 when sun at zenith) |
|---|
| 51 | c (CORRECTED for high zenith angle (atmosphere), unlike mu0) |
|---|
| 52 | c PSEC =1./PRMU |
|---|
| 53 | |
|---|
| 54 | C ==== OUTPUTS === |
|---|
| 55 | c |
|---|
| 56 | c PFD : downward flux in spectral band #INU in a given mesh |
|---|
| 57 | c (normalized to the total incident flux at the top of the atmosphere) |
|---|
| 58 | c PFU : upward flux in specatral band #INU in a given mesh |
|---|
| 59 | c (normalized to the total incident flux at the top of the atmosphere) |
|---|
| 60 | C |
|---|
| 61 | C |
|---|
| 62 | C METHOD. |
|---|
| 63 | C ------- |
|---|
| 64 | C |
|---|
| 65 | C Computes continuum fluxes corresponding to aerosoL |
|---|
| 66 | C Or/and rayleigh scattering (no molecular gas absorption) |
|---|
| 67 | C |
|---|
| 68 | C----------------------------------------------------------------------- |
|---|
| 69 | C |
|---|
| 70 | C |
|---|
| 71 | C----------------------------------------------------------------------- |
|---|
| 72 | C |
|---|
| 73 | |
|---|
| 74 | C ARGUMENTS |
|---|
| 75 | C --------- |
|---|
| 76 | INTEGER KDLON, KFLEV, KNU |
|---|
| 77 | REAL aerosol(NDLO2,KFLEV,naerkind), albedo(NDLO2,2), |
|---|
| 78 | S PDSIG(NDLO2,KFLEV),PSEC(NDLO2) |
|---|
| 79 | REAL PPSOL(NDLO2) |
|---|
| 80 | REAL PFD(NDLO2,KFLEV+1),PFU(NDLO2,KFLEV+1) |
|---|
| 81 | REAL PRMU(NDLO2) |
|---|
| 82 | |
|---|
| 83 | C LOCAL ARRAYS |
|---|
| 84 | C ------------ |
|---|
| 85 | |
|---|
| 86 | INTEGER jk,ja,jl,jae, jkl,jklp1,jkm1,jaj |
|---|
| 87 | REAL ZTRAY, ZRATIO,ZGAR, ZFF |
|---|
| 88 | real zfacoa,zcorae |
|---|
| 89 | real ZMUE, zgap,zbmu0, zww,zto,zden,zmu1,zbmu1,zden1,zre11 |
|---|
| 90 | |
|---|
| 91 | REAL ZC1I(NDLON,NFLEV+1), ZGG(NDLON), ZREF(NDLON) |
|---|
| 92 | S , ZRE1(NDLON), ZRE2(NDLON) |
|---|
| 93 | S , ZRMUZ(NDLON), ZRNEB(NDLON), ZR21(NDLON) |
|---|
| 94 | S , ZR23(NDLON), ZSS1(NDLON), ZTO1(NDLON), ZTR(NDLON,2,NFLEV+1) |
|---|
| 95 | S , ZTR1(NDLON), ZTR2(NDLON), ZW(NDLON) |
|---|
| 96 | |
|---|
| 97 | REAL ZRAY1(NDLO2,NFLEV+1), ZRAY2(NDLO2,NFLEV+1) |
|---|
| 98 | s , ZREFZ(NDLO2,2,NFLEV+1) |
|---|
| 99 | S , ZRMUE(NDLO2,NFLEV+1) |
|---|
| 100 | S , ZCGAZ(NDLO2,NFLEV),ZPIZAZ(NDLO2,NFLEV),ZTAUAZ(NDLO2,NFLEV) |
|---|
| 101 | |
|---|
| 102 | REAL ZRAYL(NDLON) |
|---|
| 103 | S , ZRJ(NDLON,6,NFLEV+1) |
|---|
| 104 | S , ZRK(NDLON,6,NFLEV+1) |
|---|
| 105 | S , ZTRA1(NDLON,NFLEV+1), ZTRA2(NDLON,NFLEV+1) |
|---|
| 106 | |
|---|
| 107 | |
|---|
| 108 | real ray,coefsizew |
|---|
| 109 | |
|---|
| 110 | c Function |
|---|
| 111 | c -------- |
|---|
| 112 | real CVMGT |
|---|
| 113 | |
|---|
| 114 | C -------------------------------- |
|---|
| 115 | C OPTICAL PARAMETERS FOR AEROSOLS |
|---|
| 116 | C ------------------------------- |
|---|
| 117 | C |
|---|
| 118 | DO JK = 1 , nlaylte+1 |
|---|
| 119 | DO JA = 1 , 6 |
|---|
| 120 | DO JL = 1 , KDLON |
|---|
| 121 | ZRJ(JL,JA,JK) = 0. |
|---|
| 122 | ZRK(JL,JA,JK) = 0. |
|---|
| 123 | END DO |
|---|
| 124 | END DO |
|---|
| 125 | END DO |
|---|
| 126 | |
|---|
| 127 | c Computing TOTAL single scattering parameters by adding |
|---|
| 128 | c properties of all the NAERKIND kind of aerosols |
|---|
| 129 | |
|---|
| 130 | DO JK = 1 , nlaylte |
|---|
| 131 | DO JL = 1 , KDLON |
|---|
| 132 | ZCGAZ(JL,JK) = 0. |
|---|
| 133 | ZPIZAZ(JL,JK) = 0. |
|---|
| 134 | ZTAUAZ(JL,JK) = 0. |
|---|
| 135 | END DO |
|---|
| 136 | DO 106 JAE=1,naerkind |
|---|
| 137 | DO 105 JL = 1 , KDLON |
|---|
| 138 | c Mean Extinction optical depth in the spectral band |
|---|
| 139 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 140 | ZTAUAZ(JL,JK)=ZTAUAZ(JL,JK) |
|---|
| 141 | S +aerosol(JL,JK,JAE)*QVISsQREF(KNU,JAE) |
|---|
| 142 | c Single scattering albedo |
|---|
| 143 | c ~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 144 | c TEST : to account for the varying w with varying crystal size |
|---|
| 145 | if (activice.and.JAE.eq.naerkind.and.KNU.eq.2) then |
|---|
| 146 | ray=min( max(rice(JL,JK)*1.e+6, 1.),10.) |
|---|
| 147 | coefsizew=(0.0001417*ray**2.-0.00328*ray+0.99667) |
|---|
| 148 | & /omegavis(KNU,JAE) |
|---|
| 149 | else |
|---|
| 150 | coefsizew=1. |
|---|
| 151 | endif |
|---|
| 152 | ZPIZAZ(JL,JK)=ZPIZAZ(JL,JK)+aerosol(JL,JK,JAE) |
|---|
| 153 | S * QVISsQREF(KNU,JAE)*omegavis(KNU,JAE)*coefsizew |
|---|
| 154 | c Assymetry factor |
|---|
| 155 | c ~~~~~~~~~~~~~~~~ |
|---|
| 156 | ZCGAZ(JL,JK) = ZCGAZ(JL,JK) +aerosol(JL,JK,JAE) |
|---|
| 157 | S * QVISsQREF(KNU,JAE)*omegavis(KNU,JAE)*gvis(KNU,JAE) |
|---|
| 158 | 105 CONTINUE |
|---|
| 159 | 106 CONTINUE |
|---|
| 160 | END DO |
|---|
| 161 | C |
|---|
| 162 | DO JK = 1 , nlaylte |
|---|
| 163 | DO JL = 1 , KDLON |
|---|
| 164 | ZCGAZ(JL,JK) = CVMGT( 0., ZCGAZ(JL,JK) / ZPIZAZ(JL,JK), |
|---|
| 165 | S (ZPIZAZ(JL,JK).EQ.0) ) |
|---|
| 166 | ZPIZAZ(JL,JK) = CVMGT( 1., ZPIZAZ(JL,JK) / ZTAUAZ(JL,JK), |
|---|
| 167 | S (ZTAUAZ(JL,JK).EQ.0) ) |
|---|
| 168 | END DO |
|---|
| 169 | END DO |
|---|
| 170 | |
|---|
| 171 | C -------------------------------- |
|---|
| 172 | C INCLUDING RAYLEIGH SCATERRING |
|---|
| 173 | C ------------------------------- |
|---|
| 174 | if (rayleigh) then |
|---|
| 175 | |
|---|
| 176 | call swrayleigh(kdlon,knu,ppsol,prmu,ZRAYL) |
|---|
| 177 | |
|---|
| 178 | c Modifying mean aerosol parameters to account rayleigh scat by gas: |
|---|
| 179 | |
|---|
| 180 | DO JK = 1 , nlaylte |
|---|
| 181 | DO JL = 1 , KDLON |
|---|
| 182 | c Rayleigh opacity in each layer : |
|---|
| 183 | ZTRAY = ZRAYL(JL) * PDSIG(JL,JK) |
|---|
| 184 | c ratio Tau(rayleigh) / Tau (total) |
|---|
| 185 | ZRATIO = ZTRAY / (ZTRAY + ZTAUAZ(JL,JK)) |
|---|
| 186 | ZGAR = ZCGAZ(JL,JK) |
|---|
| 187 | ZFF = ZGAR * ZGAR |
|---|
| 188 | ZTAUAZ(JL,JK)=ZTRAY+ZTAUAZ(JL,JK)*(1.-ZPIZAZ(JL,JK)*ZFF) |
|---|
| 189 | ZCGAZ(JL,JK) = ZGAR * (1. - ZRATIO) / (1. + ZGAR) |
|---|
| 190 | ZPIZAZ(JL,JK) =ZRATIO+(1.-ZRATIO)*ZPIZAZ(JL,JK)*(1.-ZFF) |
|---|
| 191 | S / (1. -ZPIZAZ(JL,JK) * ZFF) |
|---|
| 192 | END DO |
|---|
| 193 | END DO |
|---|
| 194 | end if |
|---|
| 195 | |
|---|
| 196 | |
|---|
| 197 | C ---------------------------------------------- |
|---|
| 198 | C TOTAL EFFECTIVE CLOUDINESS ABOVE A GIVEN LEVEL |
|---|
| 199 | C ---------------------------------------------- |
|---|
| 200 | C |
|---|
| 201 | 200 CONTINUE |
|---|
| 202 | |
|---|
| 203 | DO JL = 1 , KDLON |
|---|
| 204 | ZR23(JL) = 0. |
|---|
| 205 | ZC1I(JL,nlaylte+1) = 0. |
|---|
| 206 | END DO |
|---|
| 207 | |
|---|
| 208 | DO JK = 1 , nlaylte |
|---|
| 209 | JKL = nlaylte+1 - JK |
|---|
| 210 | JKLP1 = JKL + 1 |
|---|
| 211 | DO JL = 1 , KDLON |
|---|
| 212 | ZFACOA = 1.-ZPIZAZ(JL,JKL)*ZCGAZ(JL,JKL)*ZCGAZ(JL,JKL) |
|---|
| 213 | ZCORAE = ZFACOA * ZTAUAZ(JL,JKL) * PSEC(JL) |
|---|
| 214 | ZR21(JL) = EXP(-ZCORAE ) |
|---|
| 215 | ZSS1(JL) = 1.0-ZR21(JL) |
|---|
| 216 | ZC1I(JL,JKL) = 1.0-(1.0-ZSS1(JL))*(1.0-ZC1I(JL,JKLP1)) |
|---|
| 217 | END DO |
|---|
| 218 | END DO |
|---|
| 219 | |
|---|
| 220 | C ----------------------------------------------- |
|---|
| 221 | C REFLECTIVITY/TRANSMISSIVITY FOR PURE SCATTERING |
|---|
| 222 | C ----------------------------------------------- |
|---|
| 223 | C |
|---|
| 224 | DO JL = 1 , KDLON |
|---|
| 225 | ZRAY1(JL,nlaylte+1) = 0. |
|---|
| 226 | ZRAY2(JL,nlaylte+1) = 0. |
|---|
| 227 | ZREFZ(JL,2,1) = albedo(JL,KNU) |
|---|
| 228 | ZREFZ(JL,1,1) = albedo(JL,KNU) |
|---|
| 229 | ZTRA1(JL,nlaylte+1) = 1. |
|---|
| 230 | ZTRA2(JL,nlaylte+1) = 1. |
|---|
| 231 | END DO |
|---|
| 232 | |
|---|
| 233 | DO JK = 2 , nlaylte+1 |
|---|
| 234 | JKM1 = JK-1 |
|---|
| 235 | DO 342 JL = 1 , KDLON |
|---|
| 236 | ZRNEB(JL)= 1.e-5 ! used to be "cloudiness" (PCLDSW in Morcrette) |
|---|
| 237 | |
|---|
| 238 | ZRE1(JL)=0. |
|---|
| 239 | ZTR1(JL)=0. |
|---|
| 240 | ZRE2(JL)=0. |
|---|
| 241 | ZTR2(JL)=0. |
|---|
| 242 | |
|---|
| 243 | C EQUIVALENT ZENITH ANGLE |
|---|
| 244 | c ~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 245 | ZMUE = (1.-ZC1I(JL,JK)) * PSEC(JL) |
|---|
| 246 | S + ZC1I(JL,JK) * 1.66 |
|---|
| 247 | ZRMUE(JL,JK) = 1./ZMUE |
|---|
| 248 | |
|---|
| 249 | C ------------------------------------------------------------------ |
|---|
| 250 | C REFLECT./TRANSMISSIVITY DUE TO AEROSOLS (and rayleigh ?) |
|---|
| 251 | C ------------------------------------------------------------------ |
|---|
| 252 | |
|---|
| 253 | ZGAP = ZCGAZ(JL,JKM1) |
|---|
| 254 | ZBMU0 = 0.5 - 0.75 * ZGAP / ZMUE |
|---|
| 255 | ZWW =ZPIZAZ(JL,JKM1) |
|---|
| 256 | ZTO = ZTAUAZ(JL,JKM1) |
|---|
| 257 | ZDEN = 1. + (1. - ZWW + ZBMU0 * ZWW) * ZTO * ZMUE |
|---|
| 258 | S + (1-ZWW) * (1. - ZWW +2.*ZBMU0*ZWW)*ZTO*ZTO*ZMUE*ZMUE |
|---|
| 259 | ZRAY1(JL,JKM1) = ZBMU0 * ZWW * ZTO * ZMUE / ZDEN |
|---|
| 260 | ZTRA1(JL,JKM1) = 1. / ZDEN |
|---|
| 261 | C |
|---|
| 262 | ZMU1 = 0.5 |
|---|
| 263 | ZBMU1 = 0.5 - 0.75 * ZGAP * ZMU1 |
|---|
| 264 | ZDEN1= 1. + (1. - ZWW + ZBMU1 * ZWW) * ZTO / ZMU1 |
|---|
| 265 | S + (1-ZWW) * (1. - ZWW +2.*ZBMU1*ZWW)*ZTO*ZTO/ZMU1/ZMU1 |
|---|
| 266 | ZRAY2(JL,JKM1) = ZBMU1 * ZWW * ZTO / ZMU1 / ZDEN1 |
|---|
| 267 | ZTRA2(JL,JKM1) = 1. / ZDEN1 |
|---|
| 268 | |
|---|
| 269 | ZGG(JL) = ZCGAZ(JL,JKM1) |
|---|
| 270 | ZW(JL) =ZPIZAZ(JL,JKM1) |
|---|
| 271 | ZREF(JL) = ZREFZ(JL,1,JKM1) |
|---|
| 272 | ZRMUZ(JL) = ZRMUE(JL,JK) |
|---|
| 273 | ZTO1(JL) = ZTAUAZ(JL,JKM1)/ZPIZAZ(JL,JKM1) |
|---|
| 274 | |
|---|
| 275 | 342 CONTINUE |
|---|
| 276 | |
|---|
| 277 | C |
|---|
| 278 | CALL DEDD ( KDLON |
|---|
| 279 | S , ZGG,ZREF,ZRMUZ,ZTO1,ZW |
|---|
| 280 | S , ZRE1,ZRE2,ZTR1,ZTR2 ) |
|---|
| 281 | C |
|---|
| 282 | DO JL = 1 , KDLON |
|---|
| 283 | C |
|---|
| 284 | ZREFZ(JL,1,JK) = (1.-ZRNEB(JL)) * (ZRAY1(JL,JKM1) |
|---|
| 285 | S + ZREFZ(JL,1,JKM1) * ZTRA1(JL,JKM1) |
|---|
| 286 | S * ZTRA2(JL,JKM1) |
|---|
| 287 | S / (1.-ZRAY2(JL,JKM1)*ZREFZ(JL,1,JKM1))) |
|---|
| 288 | S + ZRNEB(JL) * ZRE2(JL) |
|---|
| 289 | C |
|---|
| 290 | ZTR(JL,1,JKM1) = ZRNEB(JL) * ZTR2(JL) + (ZTRA1(JL,JKM1) |
|---|
| 291 | S / (1.-ZRAY2(JL,JKM1)*ZREFZ(JL,1,JKM1))) |
|---|
| 292 | S * (1.-ZRNEB(JL)) |
|---|
| 293 | C |
|---|
| 294 | ZREFZ(JL,2,JK) = (1.-ZRNEB(JL)) * (ZRAY1(JL,JKM1) |
|---|
| 295 | S + ZREFZ(JL,2,JKM1) * ZTRA1(JL,JKM1) |
|---|
| 296 | S * ZTRA2(JL,JKM1) ) |
|---|
| 297 | S + ZRNEB(JL) * ZRE1(JL) |
|---|
| 298 | C |
|---|
| 299 | ZTR(JL,2,JKM1) = ZRNEB(JL) * ZTR1(JL) |
|---|
| 300 | S + ZTRA1(JL,JKM1) * (1.-ZRNEB(JL)) |
|---|
| 301 | C |
|---|
| 302 | END DO |
|---|
| 303 | END DO |
|---|
| 304 | C |
|---|
| 305 | C |
|---|
| 306 | C ------------------------------------------------------------------ |
|---|
| 307 | C |
|---|
| 308 | C * 3.5 REFLECT./TRANSMISSIVITY BETWEEN SURFACE AND LEVEL |
|---|
| 309 | C ------------------------------------------------- |
|---|
| 310 | C |
|---|
| 311 | 350 CONTINUE |
|---|
| 312 | C |
|---|
| 313 | IF (KNU.EQ.1) THEN |
|---|
| 314 | JAJ = 2 |
|---|
| 315 | DO 351 JL = 1 , KDLON |
|---|
| 316 | ZRJ(JL,JAJ,nlaylte+1) = 1. |
|---|
| 317 | ZRK(JL,JAJ,nlaylte+1) = ZREFZ(JL, 1,nlaylte+1) |
|---|
| 318 | 351 CONTINUE |
|---|
| 319 | C |
|---|
| 320 | DO 353 JK = 1 , nlaylte |
|---|
| 321 | JKL = nlaylte+1 - JK |
|---|
| 322 | JKLP1 = JKL + 1 |
|---|
| 323 | DO 352 JL = 1 , KDLON |
|---|
| 324 | ZRE11= ZRJ(JL,JAJ,JKLP1) * ZTR(JL, 1,JKL) |
|---|
| 325 | ZRJ(JL,JAJ,JKL) = ZRE11 |
|---|
| 326 | ZRK(JL,JAJ,JKL) = ZRE11 * ZREFZ(JL, 1,JKL) |
|---|
| 327 | 352 CONTINUE |
|---|
| 328 | 353 CONTINUE |
|---|
| 329 | 354 CONTINUE |
|---|
| 330 | C |
|---|
| 331 | ELSE |
|---|
| 332 | C |
|---|
| 333 | DO 358 JAJ = 1 , 2 |
|---|
| 334 | DO 355 JL = 1 , KDLON |
|---|
| 335 | ZRJ(JL,JAJ,nlaylte+1) = 1. |
|---|
| 336 | ZRK(JL,JAJ,nlaylte+1) = ZREFZ(JL,JAJ,nlaylte+1) |
|---|
| 337 | 355 CONTINUE |
|---|
| 338 | C |
|---|
| 339 | DO 357 JK = 1 , nlaylte |
|---|
| 340 | JKL = nlaylte+1 - JK |
|---|
| 341 | JKLP1 = JKL + 1 |
|---|
| 342 | DO 356 JL = 1 , KDLON |
|---|
| 343 | ZRE11= ZRJ(JL,JAJ,JKLP1) * ZTR(JL,JAJ,JKL) |
|---|
| 344 | ZRJ(JL,JAJ,JKL) = ZRE11 |
|---|
| 345 | ZRK(JL,JAJ,JKL) = ZRE11 * ZREFZ(JL,JAJ,JKL) |
|---|
| 346 | 356 CONTINUE |
|---|
| 347 | 357 CONTINUE |
|---|
| 348 | 358 CONTINUE |
|---|
| 349 | END IF |
|---|
| 350 | |
|---|
| 351 | C |
|---|
| 352 | C |
|---|
| 353 | C |
|---|
| 354 | C ------------------------------------------------------------------ |
|---|
| 355 | C --------------- |
|---|
| 356 | C DOWNWARD FLUXES |
|---|
| 357 | C --------------- |
|---|
| 358 | C |
|---|
| 359 | JAJ = 2 |
|---|
| 360 | |
|---|
| 361 | do JK = 1 , nlaylte+1 |
|---|
| 362 | JKL = nlaylte+1 - JK + 1 |
|---|
| 363 | DO JL = 1 , KDLON |
|---|
| 364 | PFD(JL,JKL) = ZRJ(JL,JAJ,JKL) * sunfr(KNU) |
|---|
| 365 | end do |
|---|
| 366 | end do |
|---|
| 367 | C |
|---|
| 368 | C ------------- |
|---|
| 369 | C UPWARD FLUXES |
|---|
| 370 | C ------------- |
|---|
| 371 | DO JK = 1 , nlaylte+1 |
|---|
| 372 | DO JL = 1 , KDLON |
|---|
| 373 | c ZRK = upward flux / incident top flux |
|---|
| 374 | PFU(JL,JK) = ZRK(JL,JAJ,JK) * sunfr(KNU) |
|---|
| 375 | END DO |
|---|
| 376 | END DO |
|---|
| 377 | |
|---|
| 378 | C |
|---|
| 379 | RETURN |
|---|
| 380 | END |
|---|
