| 1 | MODULE flusv_mod |
|---|
| 2 | |
|---|
| 3 | IMPLICIT NONE |
|---|
| 4 | |
|---|
| 5 | CONTAINS |
|---|
| 6 | |
|---|
| 7 | SUBROUTINE flusv(KDLON,nsf,n,omega,g,tau,emis,bh,bsol,fah,fdh) |
|---|
| 8 | use dimradmars_mod, only: ndlo2, ndlon, nflev |
|---|
| 9 | IMPLICIT NONE |
|---|
| 10 | c....................................................................... |
|---|
| 11 | c |
|---|
| 12 | c compute the upward and downward fluxes at the interface between n layers |
|---|
| 13 | c * in the infrared |
|---|
| 14 | c * B is a linear function of $\tau$ in each layer |
|---|
| 15 | c * B at the surface can be different than what corresponds to the profile |
|---|
| 16 | c in the n-th layer |
|---|
| 17 | c * the work hypothes isthat we have two isotropic fluxes for each |
|---|
| 18 | c hemisphere ("hemispheric constant") + "technical source function" |
|---|
| 19 | c (see Toon et al. 1988) |
|---|
| 20 | c * the downward flux at the top of the atmosphere is zero |
|---|
| 21 | c * layers are numbered from the top of the atmosphere to the ground |
|---|
| 22 | c |
|---|
| 23 | c in : * KDLON ---> vectorisation dimension |
|---|
| 24 | c * nsf ---> nsf=0 ==> "hemispheric constant" |
|---|
| 25 | c nsf>0 ==> "hemispheric constant" + "source function" |
|---|
| 26 | c * n ---> number of layers |
|---|
| 27 | c * omega(i) ---> single scattering albedo for the i-th layer |
|---|
| 28 | c * g(i) ---> asymmetry parameter for the i-th layer |
|---|
| 29 | c * tau(i) ---> optical thickness of the i-th layer |
|---|
| 30 | c * emis ---> ground emissivity |
|---|
| 31 | c * bh(i) ---> black body luminance at the top of the i-th layer, |
|---|
| 32 | c bh(n+1) for the ground value which |
|---|
| 33 | c corresponds to the profile for the n-th layer |
|---|
| 34 | c * bsol ---> black body luminance of the ground |
|---|
| 35 | c |
|---|
| 36 | c out : * fah(i) ---> upward flux at the top of the i-th layer, |
|---|
| 37 | c fah(n+1) for the ground |
|---|
| 38 | c * fdh(i) ---> downward flux at the top of the i-th layer, |
|---|
| 39 | c fdh(n+1) for the ground |
|---|
| 40 | c |
|---|
| 41 | c....................................................................... |
|---|
| 42 | c arguments |
|---|
| 43 | c |
|---|
| 44 | INTEGER,INTENT(IN) :: KDLON,nsf,n |
|---|
| 45 | REAL,INTENT(IN) :: omega(NDLO2,n),g(NDLO2,n) |
|---|
| 46 | REAL,INTENT(IN) :: tau(NDLO2,n),emis(NDLO2) |
|---|
| 47 | REAL,INTENT(IN) :: bh(NDLO2,n+1),bsol(NDLO2) |
|---|
| 48 | REAL,INTENT(OUT) :: fah(NDLO2,n+1),fdh(NDLO2,n+1) |
|---|
| 49 | c....................................................................... |
|---|
| 50 | c local variables |
|---|
| 51 | c |
|---|
| 52 | REAL,PARAMETER :: pi=3.141592653589793E+0 |
|---|
| 53 | INTEGER iv,i,j |
|---|
| 54 | REAL beta,gama1,gama2,amu1,grgama,b0,b1 |
|---|
| 55 | REAL a(NDLON,4*nflev),b(NDLON,4*nflev) |
|---|
| 56 | & ,d(NDLON,4*nflev),e(NDLON,4*nflev) |
|---|
| 57 | & ,y(NDLON,4*nflev) |
|---|
| 58 | & ,alambda(NDLON,2*nflev) |
|---|
| 59 | & ,e1(NDLON,2*nflev),e2(NDLON,2*nflev) |
|---|
| 60 | & ,e3(NDLON,2*nflev),e4(NDLON,2*nflev) |
|---|
| 61 | & ,cah(NDLON,2*nflev),cab(NDLON,2*nflev) |
|---|
| 62 | & ,cdh(NDLON,2*nflev),cdb(NDLON,2*nflev) |
|---|
| 63 | REAL grg(NDLON,2*nflev),grh(NDLON,2*nflev) |
|---|
| 64 | & ,grj(NDLON,2*nflev),grk(NDLON,2*nflev) |
|---|
| 65 | & ,alpha1(NDLON,2*nflev),alpha2(NDLON,2*nflev) |
|---|
| 66 | & ,sigma1(NDLON,2*nflev),sigma2(NDLON,2*nflev) |
|---|
| 67 | INTEGER,PARAMETER :: nq=8 |
|---|
| 68 | REAL,PARAMETER :: x(nq) = |
|---|
| 69 | & (/1.9855071751231860E-2 , 0.1016667612931866E+0 |
|---|
| 70 | & , 0.2372337950418355E+0 , 0.4082826787521751E+0 |
|---|
| 71 | & , 0.5917173212478250E+0 , 0.7627662049581645E+0 |
|---|
| 72 | & , 0.8983332387068134E+0 , 0.9801449282487682E+0/) |
|---|
| 73 | REAL,PARAMETER :: w(nq) = |
|---|
| 74 | & (/5.0614268145185310E-2 , 0.1111905172266872E+0 |
|---|
| 75 | & , 0.1568533229389437E+0 , 0.1813418916891810E+0 |
|---|
| 76 | & , 0.1813418916891810E+0 , 0.1568533229389437E+0 |
|---|
| 77 | & , 0.1111905172266872E+0 , 5.0614268145185310E-2/) |
|---|
| 78 | REAL :: gri(NDLON,nq) |
|---|
| 79 | c....................................................................... |
|---|
| 80 | c |
|---|
| 81 | c....................................................................... |
|---|
| 82 | do i=1,n |
|---|
| 83 | do iv=1,KDLON |
|---|
| 84 | beta=(1.E+0-g(iv,i))/2.E+0 |
|---|
| 85 | gama1=2.E+0*(1.E+0-omega(iv,i)*(1.E+0-beta)) |
|---|
| 86 | gama2=2.E+0*omega(iv,i)*beta |
|---|
| 87 | amu1=5.E-1 |
|---|
| 88 | alambda(iv,i)=sqrt(gama1**2-gama2**2) |
|---|
| 89 | grgama=(gama1-alambda(iv,i))/gama2 |
|---|
| 90 | c |
|---|
| 91 | c small hack here : if the optical depth of a layer is too small, |
|---|
| 92 | c $dB \over d\tau$ becomes very large and the scheme fails. |
|---|
| 93 | c In those cases we assume an isothermal layer. |
|---|
| 94 | c |
|---|
| 95 | if (tau(iv,i).gt.1.E-3) then |
|---|
| 96 | b0=bh(iv,i) |
|---|
| 97 | b1=(bh(iv,i+1)-b0)/tau(iv,i) |
|---|
| 98 | else |
|---|
| 99 | b0=(bh(iv,i)+bh(iv,i+1))/2.E+0 |
|---|
| 100 | b1=0.E+0 |
|---|
| 101 | endif |
|---|
| 102 | c |
|---|
| 103 | e1(iv,i)=1.E+0+grgama*exp(-alambda(iv,i)*tau(iv,i)) |
|---|
| 104 | e2(iv,i)=1.E+0-grgama*exp(-alambda(iv,i)*tau(iv,i)) |
|---|
| 105 | e3(iv,i)=grgama+exp(-alambda(iv,i)*tau(iv,i)) |
|---|
| 106 | e4(iv,i)=grgama-exp(-alambda(iv,i)*tau(iv,i)) |
|---|
| 107 | cah(iv,i)=2.E+0*pi*amu1*(b0+b1/(gama1+gama2)) |
|---|
| 108 | cab(iv,i)=2.E+0*pi*amu1*(b0+b1*(tau(iv,i)+1.E+0/(gama1+gama2))) |
|---|
| 109 | cdh(iv,i)=2.E+0*pi*amu1*(b0-b1/(gama1+gama2)) |
|---|
| 110 | cdb(iv,i)=2.E+0*pi*amu1*(b0+b1*(tau(iv,i)-1.E+0/(gama1+gama2))) |
|---|
| 111 | c |
|---|
| 112 | grg(iv,i)=(1.E+0/amu1-alambda(iv,i)) |
|---|
| 113 | grh(iv,i)=grgama*(alambda(iv,i)+1.E+0/amu1) |
|---|
| 114 | grj(iv,i)=grh(iv,i) |
|---|
| 115 | grk(iv,i)=grg(iv,i) |
|---|
| 116 | alpha1(iv,i)=2.E+0*pi*(b0+b1*(1.E+0/(gama1+gama2)-amu1)) |
|---|
| 117 | alpha2(iv,i)=2.E+0*pi*b1 |
|---|
| 118 | sigma1(iv,i)=2.E+0*pi*(b0-b1*(1.E+0/(gama1+gama2)-amu1)) |
|---|
| 119 | sigma2(iv,i)=alpha2(iv,i) |
|---|
| 120 | c |
|---|
| 121 | enddo ! of do iv=1,KDLON |
|---|
| 122 | enddo ! of do i=1,n |
|---|
| 123 | c....................................................................... |
|---|
| 124 | do iv=1,KDLON |
|---|
| 125 | a(iv,1)=0.E+0 |
|---|
| 126 | b(iv,1)=e1(iv,1) |
|---|
| 127 | d(iv,1)=-e2(iv,1) |
|---|
| 128 | e(iv,1)=-cdh(iv,1) |
|---|
| 129 | enddo |
|---|
| 130 | c |
|---|
| 131 | do i=1,n-1 |
|---|
| 132 | j=2*i+1 |
|---|
| 133 | do iv=1,KDLON |
|---|
| 134 | a(iv,j)=e2(iv,i)*e3(iv,i)-e4(iv,i)*e1(iv,i) |
|---|
| 135 | b(iv,j)=e1(iv,i)*e1(iv,i+1)-e3(iv,i)*e3(iv,i+1) |
|---|
| 136 | d(iv,j)=e3(iv,i)*e4(iv,i+1)-e1(iv,i)*e2(iv,i+1) |
|---|
| 137 | e(iv,j)=e3(iv,i)*(cah(iv,i+1)-cab(iv,i)) |
|---|
| 138 | & +e1(iv,i)*(cdb(iv,i)-cdh(iv,i+1)) |
|---|
| 139 | enddo |
|---|
| 140 | enddo ! of do i=1,n-1 |
|---|
| 141 | c |
|---|
| 142 | do i=1,n-1 |
|---|
| 143 | j=2*i |
|---|
| 144 | do iv=1,KDLON |
|---|
| 145 | a(iv,j)=e2(iv,i+1)*e1(iv,i)-e3(iv,i)*e4(iv,i+1) |
|---|
| 146 | b(iv,j)=e2(iv,i)*e2(iv,i+1)-e4(iv,i)*e4(iv,i+1) |
|---|
| 147 | d(iv,j)=e1(iv,i+1)*e4(iv,i+1)-e2(iv,i+1)*e3(iv,i+1) |
|---|
| 148 | e(iv,j)=e2(iv,i+1)*(cah(iv,i+1)-cab(iv,i)) |
|---|
| 149 | & +e4(iv,i+1)*(cdb(iv,i)-cdh(iv,i+1)) |
|---|
| 150 | enddo |
|---|
| 151 | enddo ! of do i=1,n-1 |
|---|
| 152 | c |
|---|
| 153 | j=2*n |
|---|
| 154 | do iv=1,KDLON |
|---|
| 155 | a(iv,j)=e1(iv,n)-(1.E+0-emis(iv))*e3(iv,n) |
|---|
| 156 | b(iv,j)=e2(iv,n)-(1.E+0-emis(iv))*e4(iv,n) |
|---|
| 157 | d(iv,j)=0.E+0 |
|---|
| 158 | e(iv,j)=emis(iv)*pi*bsol(iv)-cab(iv,n) |
|---|
| 159 | & +(1.E+0-emis(iv))*cdb(iv,n) |
|---|
| 160 | enddo |
|---|
| 161 | c....................................................................... |
|---|
| 162 | call sys3v(KDLON,2*n,a,b,d,e,y) |
|---|
| 163 | c....................................................................... |
|---|
| 164 | do i=1,n |
|---|
| 165 | do iv=1,KDLON |
|---|
| 166 | grg(iv,i)=grg(iv,i)*(y(iv,2*i-1)+y(iv,2*i)) |
|---|
| 167 | grh(iv,i)=grh(iv,i)*(y(iv,2*i-1)-y(iv,2*i)) |
|---|
| 168 | grj(iv,i)=grj(iv,i)*(y(iv,2*i-1)+y(iv,2*i)) |
|---|
| 169 | grk(iv,i)=grk(iv,i)*(y(iv,2*i-1)-y(iv,2*i)) |
|---|
| 170 | enddo |
|---|
| 171 | enddo |
|---|
| 172 | c....................................................................... |
|---|
| 173 | c values of "hemispheric constant" fluxes |
|---|
| 174 | c |
|---|
| 175 | IF (nsf.eq.0) THEN |
|---|
| 176 | do i=1,n |
|---|
| 177 | do iv=1,KDLON |
|---|
| 178 | fah(iv,i)=e3(iv,i)*y(iv,2*i-1)-e4(iv,i)*y(iv,2*i)+cah(iv,i) |
|---|
| 179 | fdh(iv,i)=e1(iv,i)*y(iv,2*i-1)-e2(iv,i)*y(iv,2*i)+cdh(iv,i) |
|---|
| 180 | enddo |
|---|
| 181 | enddo |
|---|
| 182 | do iv=1,KDLON |
|---|
| 183 | fah(iv,n+1)=e1(iv,n)*y(iv,2*n-1)+e2(iv,n)*y(iv,2*n)+cab(iv,n) |
|---|
| 184 | fdh(iv,n+1)=e3(iv,n)*y(iv,2*n-1)+e4(iv,n)*y(iv,2*n)+cdb(iv,n) |
|---|
| 185 | enddo |
|---|
| 186 | ELSE |
|---|
| 187 | c....................................................................... |
|---|
| 188 | c going to the "source function" |
|---|
| 189 | c |
|---|
| 190 | c apply a quadrature over nq (fixed parameter) points |
|---|
| 191 | c x is the vector of the \mu of the quadrature |
|---|
| 192 | c w is the vector of corresponding weights |
|---|
| 193 | c x() et w() are fixed parameters |
|---|
| 194 | c |
|---|
| 195 | c....................................................................... |
|---|
| 196 | c start from the top and go down along the nq angles to compute all |
|---|
| 197 | c downward fluxes |
|---|
| 198 | c |
|---|
| 199 | do j=1,nq |
|---|
| 200 | do iv=1,KDLON |
|---|
| 201 | gri(iv,j)=0.E+0 |
|---|
| 202 | enddo |
|---|
| 203 | enddo |
|---|
| 204 | do iv=1,KDLON |
|---|
| 205 | fdh(iv,1)=0.E+0 |
|---|
| 206 | enddo |
|---|
| 207 | do i=1,n |
|---|
| 208 | do j=1,nq |
|---|
| 209 | do iv=1,KDLON |
|---|
| 210 | gri(iv,j)=gri(iv,j)*exp(-tau(iv,i)/x(j)) |
|---|
| 211 | & +grj(iv,i)/(alambda(iv,i)*x(j)+1.E+0) |
|---|
| 212 | & *(1.E+0-exp(-tau(iv,i)*(alambda(iv,i)+1.E+0/x(j)))) |
|---|
| 213 | & +grk(iv,i)/(alambda(iv,i)*x(j)-1.E+0) |
|---|
| 214 | & *(exp(-tau(iv,i)/x(j))-exp(-tau(iv,i)*alambda(iv,i))) |
|---|
| 215 | & +sigma1(iv,i)*(1.E+0-exp(-tau(iv,i)/x(j))) |
|---|
| 216 | & +sigma2(iv,i)*(x(j)*exp(-tau(iv,i)/x(j))+tau(iv,i)-x(j)) |
|---|
| 217 | enddo ! of do iv=1,KDLON |
|---|
| 218 | enddo ! of do j=1,nq |
|---|
| 219 | do iv=1,KDLON |
|---|
| 220 | fdh(iv,i+1)=0.E+0 |
|---|
| 221 | enddo |
|---|
| 222 | do j=1,nq |
|---|
| 223 | do iv=1,KDLON |
|---|
| 224 | fdh(iv,i+1)=fdh(iv,i+1)+w(j)*x(j)*gri(iv,j) |
|---|
| 225 | enddo ! of do iv=1,KDLON |
|---|
| 226 | enddo ! of do j=1,nq |
|---|
| 227 | enddo ! of do i=1,n |
|---|
| 228 | c....................................................................... |
|---|
| 229 | c apply the reflexion condition on the ground |
|---|
| 230 | c |
|---|
| 231 | do iv=1,KDLON |
|---|
| 232 | fah(iv,n+1)=(1.E+0-emis(iv))*fdh(iv,n+1)+pi*emis(iv)*bsol(iv) |
|---|
| 233 | enddo |
|---|
| 234 | do j=1,nq |
|---|
| 235 | do iv=1,KDLON |
|---|
| 236 | gri(iv,j)=2.E+0*fah(iv,n+1) |
|---|
| 237 | enddo |
|---|
| 238 | enddo |
|---|
| 239 | c....................................................................... |
|---|
| 240 | c going back up to compute all the upward fluxes |
|---|
| 241 | c |
|---|
| 242 | do i=n,1,-1 |
|---|
| 243 | do j=1,nq |
|---|
| 244 | do iv=1,KDLON |
|---|
| 245 | gri(iv,j)=gri(iv,j)*exp(-tau(iv,i)/x(j)) |
|---|
| 246 | & +grg(iv,i)/(alambda(iv,i)*x(j)-1.E+0) |
|---|
| 247 | & *(exp(-tau(iv,i)/x(j))-exp(-tau(iv,i)*alambda(iv,i))) |
|---|
| 248 | & +grh(iv,i)/(alambda(iv,i)*x(j)+1.E+0) |
|---|
| 249 | & *(1.E+0-exp(-tau(iv,i)*(alambda(iv,i)+1.E+0/x(j)))) |
|---|
| 250 | & +alpha1(iv,i)*(1.E+0-exp(-tau(iv,i)/x(j))) |
|---|
| 251 | & +alpha2(iv,i)*(x(j)-(tau(iv,i)+x(j))*exp(-tau(iv,i)/x(j))) |
|---|
| 252 | enddo ! of do iv=1,KDLON |
|---|
| 253 | enddo ! of do j=1,nq |
|---|
| 254 | do iv=1,KDLON |
|---|
| 255 | fah(iv,i)=0.E+0 |
|---|
| 256 | enddo |
|---|
| 257 | do j=1,nq |
|---|
| 258 | do iv=1,KDLON |
|---|
| 259 | fah(iv,i)=fah(iv,i)+w(j)*x(j)*gri(iv,j) |
|---|
| 260 | enddo |
|---|
| 261 | enddo ! of do j=1,nq |
|---|
| 262 | enddo ! of do i=n,1,-1 |
|---|
| 263 | c....................................................................... |
|---|
| 264 | ENDIF ! of IF (nsf.eq.0) |
|---|
| 265 | c....................................................................... |
|---|
| 266 | c |
|---|
| 267 | c....................................................................... |
|---|
| 268 | |
|---|
| 269 | END SUBROUTINE flusv |
|---|
| 270 | |
|---|
| 271 | c *************************************************************** |
|---|
| 272 | |
|---|
| 273 | |
|---|
| 274 | SUBROUTINE sys3v(KDLON,n,a,b,d,e,y) |
|---|
| 275 | use dimradmars_mod, only: ndlon, ndlo2, nflev |
|---|
| 276 | IMPLICIT NONE |
|---|
| 277 | c....................................................................... |
|---|
| 278 | c |
|---|
| 279 | c solve a tridiagonal linear system such that: |
|---|
| 280 | c |
|---|
| 281 | c | b1 d1 | | y1 | | e1 | |
|---|
| 282 | c | a2 b2 d2 | | y2 | | e2 | |
|---|
| 283 | c | a3 b3 d3 | * | y3 | = | e3 | |
|---|
| 284 | c | .... | | | | | |
|---|
| 285 | c | an bn | | yn | | en | |
|---|
| 286 | c |
|---|
| 287 | c in : * KDLON --> vectorisation dimension |
|---|
| 288 | c * n --> system size |
|---|
| 289 | c * a,b,d,e --> coefficients as shown above |
|---|
| 290 | c |
|---|
| 291 | c out : * y --> see above |
|---|
| 292 | c |
|---|
| 293 | c....................................................................... |
|---|
| 294 | c arguments |
|---|
| 295 | c |
|---|
| 296 | INTEGER,INTENT(IN) :: KDLON,n |
|---|
| 297 | REAL,INTENT(IN) :: a(NDLO2,n),b(NDLO2,n),d(NDLO2,n),e(NDLO2,n) |
|---|
| 298 | REAL,INTENT(OUT) :: y(NDLO2,n) |
|---|
| 299 | c....................................................................... |
|---|
| 300 | c local variables |
|---|
| 301 | c |
|---|
| 302 | INTEGER :: iv,i |
|---|
| 303 | REAL :: as(NDLON,4*nflev),ds(NDLON,4*nflev) |
|---|
| 304 | & ,x(NDLON,4*nflev) |
|---|
| 305 | c....................................................................... |
|---|
| 306 | c |
|---|
| 307 | c....................................................................... |
|---|
| 308 | do iv=1,KDLON |
|---|
| 309 | as(iv,n)=a(iv,n)/b(iv,n) |
|---|
| 310 | ds(iv,n)=e(iv,n)/b(iv,n) |
|---|
| 311 | enddo |
|---|
| 312 | do i=n-1,1,-1 |
|---|
| 313 | do iv=1,KDLON |
|---|
| 314 | x(iv,i)=1.E+0/(b(iv,i)-d(iv,i)*as(iv,i+1)) |
|---|
| 315 | as(iv,i)=a(iv,i)*x(iv,i) |
|---|
| 316 | ds(iv,i)=(e(iv,i)-d(iv,i)*ds(iv,i+1))*x(iv,i) |
|---|
| 317 | enddo |
|---|
| 318 | enddo |
|---|
| 319 | do iv=1,KDLON |
|---|
| 320 | y(iv,1)=ds(iv,1) |
|---|
| 321 | enddo |
|---|
| 322 | do i=2,n |
|---|
| 323 | do iv=1,KDLON |
|---|
| 324 | y(iv,i)=ds(iv,i)-as(iv,i)*y(iv,i-1) |
|---|
| 325 | enddo |
|---|
| 326 | enddo |
|---|
| 327 | c....................................................................... |
|---|
| 328 | c |
|---|
| 329 | c....................................................................... |
|---|
| 330 | |
|---|
| 331 | END SUBROUTINE sys3v |
|---|
| 332 | |
|---|
| 333 | END MODULE flusv_mod |
|---|
