| 1 | SUBROUTINE flusv(KDLON,nsf,n,omega,g,tau,emis,bh,bsol,fah,fdh) |
|---|
| 2 | use dimradmars_mod, only: ndlo2, ndlon |
|---|
| 3 | IMPLICIT NONE |
|---|
| 4 | c....................................................................... |
|---|
| 5 | c |
|---|
| 6 | c calcul des flux ascendant et descendant aux interfaces entre n couches |
|---|
| 7 | c * dans l'infrarouge |
|---|
| 8 | c * B est une fonction lineaire de $\tau$ a l'interieur de chaque couche |
|---|
| 9 | c * le B du sol peut etre different de celui qui correspond au profil |
|---|
| 10 | c de la n-ieme couche |
|---|
| 11 | c * l'hypothese est une hypothese a deux flux isotropes sur chaque |
|---|
| 12 | c hemisphere ("hemispheric constant") + "source function technique" |
|---|
| 13 | c (voir Toon et al. 1988) |
|---|
| 14 | c * le flux descendant en haut de l'atmosphere est nul |
|---|
| 15 | c * les couches sont numerotees du haut de l'atmosphere vers le sol |
|---|
| 16 | c |
|---|
| 17 | c in : * KDLON ---> dimension de vectorisation |
|---|
| 18 | c * nsf ---> nsf=0 ==> "hemispheric constant" |
|---|
| 19 | c nsf>0 ==> "hemispheric constant" + "source function" |
|---|
| 20 | c * n ---> nombre de couches |
|---|
| 21 | c * omega(i) ---> single scattering albedo pour la i-eme couche |
|---|
| 22 | c * g(i) ---> asymmetry parameter pour la i-eme couche |
|---|
| 23 | c * tau(i) ---> epaisseur optique de la i-eme couche |
|---|
| 24 | c * emis ---> emissivite du sol |
|---|
| 25 | c * bh(i) ---> luminance du corps noir en haut de la i-eme |
|---|
| 26 | c couche, bh(n+1) pour la valeur au sol qui |
|---|
| 27 | c correspond au profil de la n-ieme couche |
|---|
| 28 | c * bsol ---> luminance du corps noir au sol |
|---|
| 29 | c |
|---|
| 30 | c out : * fah(i) ---> flux ascendant en haut de la i-eme couche, |
|---|
| 31 | c fah(n+1) pour le sol |
|---|
| 32 | c * fdh(i) ---> flux descendant en haut de la i-eme couche, |
|---|
| 33 | c fdh(n+1) pour le sol |
|---|
| 34 | c |
|---|
| 35 | c....................................................................... |
|---|
| 36 | c include des dimensions locales |
|---|
| 37 | c |
|---|
| 38 | #include "dimensions.h" |
|---|
| 39 | #include "dimphys.h" |
|---|
| 40 | !#include "dimradmars.h" |
|---|
| 41 | c....................................................................... |
|---|
| 42 | c declaration des arguments |
|---|
| 43 | c |
|---|
| 44 | INTEGER KDLON,nsf,n |
|---|
| 45 | REAL omega(NDLO2,n),g(NDLO2,n),tau(NDLO2,n),emis(NDLO2) |
|---|
| 46 | &,bh(NDLO2,n+1),bsol(NDLO2),fah(NDLO2,n+1),fdh(NDLO2,n+1) |
|---|
| 47 | c....................................................................... |
|---|
| 48 | c declaration des variables locales |
|---|
| 49 | c |
|---|
| 50 | REAL pi |
|---|
| 51 | PARAMETER (pi=3.141592653589793E+0) |
|---|
| 52 | INTEGER iv,i,j |
|---|
| 53 | REAL beta,gama1,gama2,amu1,grgama,b0,b1 |
|---|
| 54 | REAL a(NDLON,4*nlayermx),b(NDLON,4*nlayermx) |
|---|
| 55 | & ,d(NDLON,4*nlayermx),e(NDLON,4*nlayermx) |
|---|
| 56 | & ,y(NDLON,4*nlayermx) |
|---|
| 57 | & ,alambda(NDLON,2*nlayermx) |
|---|
| 58 | & ,e1(NDLON,2*nlayermx),e2(NDLON,2*nlayermx) |
|---|
| 59 | & ,e3(NDLON,2*nlayermx),e4(NDLON,2*nlayermx) |
|---|
| 60 | & ,cah(NDLON,2*nlayermx),cab(NDLON,2*nlayermx) |
|---|
| 61 | & ,cdh(NDLON,2*nlayermx),cdb(NDLON,2*nlayermx) |
|---|
| 62 | REAL grg(NDLON,2*nlayermx),grh(NDLON,2*nlayermx) |
|---|
| 63 | & ,grj(NDLON,2*nlayermx),grk(NDLON,2*nlayermx) |
|---|
| 64 | & ,alpha1(NDLON,2*nlayermx),alpha2(NDLON,2*nlayermx) |
|---|
| 65 | & ,sigma1(NDLON,2*nlayermx),sigma2(NDLON,2*nlayermx) |
|---|
| 66 | INTEGER nq |
|---|
| 67 | PARAMETER (nq=8) |
|---|
| 68 | REAL x(nq),w(nq),gri(NDLON,nq) |
|---|
| 69 | DATA x/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 | DATA w/5.0614268145185310E-2 , 0.1111905172266872E+0 |
|---|
| 74 | & , 0.1568533229389437E+0 , 0.1813418916891810E+0 |
|---|
| 75 | & , 0.1813418916891810E+0 , 0.1568533229389437E+0 |
|---|
| 76 | & , 0.1111905172266872E+0 , 5.0614268145185310E-2/ |
|---|
| 77 | c....................................................................... |
|---|
| 78 | c |
|---|
| 79 | c....................................................................... |
|---|
| 80 | do 10001 i=1,n |
|---|
| 81 | do 99999 iv=1,KDLON |
|---|
| 82 | beta=(1.E+0-g(iv,i))/2.E+0 |
|---|
| 83 | gama1=2.E+0*(1.E+0-omega(iv,i)*(1.E+0-beta)) |
|---|
| 84 | gama2=2.E+0*omega(iv,i)*beta |
|---|
| 85 | amu1=5.E-1 |
|---|
| 86 | alambda(iv,i)=sqrt(gama1**2-gama2**2) |
|---|
| 87 | grgama=(gama1-alambda(iv,i))/gama2 |
|---|
| 88 | c |
|---|
| 89 | c ici une petite bidouille : si l'epaisseur optique d'une couche |
|---|
| 90 | c est trop faible, $dB \over d\tau$ devient tres grand et le schema |
|---|
| 91 | c s'ecroule. dans ces cas, on fait l'hypothese de couche isotherme. |
|---|
| 92 | c |
|---|
| 93 | if (tau(iv,i).gt.1.E-3) then |
|---|
| 94 | b0=bh(iv,i) |
|---|
| 95 | b1=(bh(iv,i+1)-b0)/tau(iv,i) |
|---|
| 96 | else |
|---|
| 97 | b0=(bh(iv,i)+bh(iv,i+1))/2.E+0 |
|---|
| 98 | b1=0.E+0 |
|---|
| 99 | endif |
|---|
| 100 | c |
|---|
| 101 | e1(iv,i)=1.E+0+grgama*exp(-alambda(iv,i)*tau(iv,i)) |
|---|
| 102 | e2(iv,i)=1.E+0-grgama*exp(-alambda(iv,i)*tau(iv,i)) |
|---|
| 103 | e3(iv,i)=grgama+exp(-alambda(iv,i)*tau(iv,i)) |
|---|
| 104 | e4(iv,i)=grgama-exp(-alambda(iv,i)*tau(iv,i)) |
|---|
| 105 | cah(iv,i)=2.E+0*pi*amu1*(b0+b1/(gama1+gama2)) |
|---|
| 106 | cab(iv,i)=2.E+0*pi*amu1*(b0+b1*(tau(iv,i)+1.E+0/(gama1+gama2))) |
|---|
| 107 | cdh(iv,i)=2.E+0*pi*amu1*(b0-b1/(gama1+gama2)) |
|---|
| 108 | cdb(iv,i)=2.E+0*pi*amu1*(b0+b1*(tau(iv,i)-1.E+0/(gama1+gama2))) |
|---|
| 109 | c |
|---|
| 110 | grg(iv,i)=(1.E+0/amu1-alambda(iv,i)) |
|---|
| 111 | grh(iv,i)=grgama*(alambda(iv,i)+1.E+0/amu1) |
|---|
| 112 | grj(iv,i)=grh(iv,i) |
|---|
| 113 | grk(iv,i)=grg(iv,i) |
|---|
| 114 | alpha1(iv,i)=2.E+0*pi*(b0+b1*(1.E+0/(gama1+gama2)-amu1)) |
|---|
| 115 | alpha2(iv,i)=2.E+0*pi*b1 |
|---|
| 116 | sigma1(iv,i)=2.E+0*pi*(b0-b1*(1.E+0/(gama1+gama2)-amu1)) |
|---|
| 117 | sigma2(iv,i)=alpha2(iv,i) |
|---|
| 118 | c |
|---|
| 119 | 99999 continue |
|---|
| 120 | 10001 continue |
|---|
| 121 | c....................................................................... |
|---|
| 122 | do 99998 iv=1,KDLON |
|---|
| 123 | a(iv,1)=0.E+0 |
|---|
| 124 | b(iv,1)=e1(iv,1) |
|---|
| 125 | d(iv,1)=-e2(iv,1) |
|---|
| 126 | e(iv,1)=-cdh(iv,1) |
|---|
| 127 | 99998 continue |
|---|
| 128 | c |
|---|
| 129 | do 10002 i=1,n-1 |
|---|
| 130 | j=2*i+1 |
|---|
| 131 | do 99997 iv=1,KDLON |
|---|
| 132 | a(iv,j)=e2(iv,i)*e3(iv,i)-e4(iv,i)*e1(iv,i) |
|---|
| 133 | b(iv,j)=e1(iv,i)*e1(iv,i+1)-e3(iv,i)*e3(iv,i+1) |
|---|
| 134 | d(iv,j)=e3(iv,i)*e4(iv,i+1)-e1(iv,i)*e2(iv,i+1) |
|---|
| 135 | e(iv,j)=e3(iv,i)*(cah(iv,i+1)-cab(iv,i)) |
|---|
| 136 | & +e1(iv,i)*(cdb(iv,i)-cdh(iv,i+1)) |
|---|
| 137 | 99997 continue |
|---|
| 138 | 10002 continue |
|---|
| 139 | c |
|---|
| 140 | do 10003 i=1,n-1 |
|---|
| 141 | j=2*i |
|---|
| 142 | do 99996 iv=1,KDLON |
|---|
| 143 | a(iv,j)=e2(iv,i+1)*e1(iv,i)-e3(iv,i)*e4(iv,i+1) |
|---|
| 144 | b(iv,j)=e2(iv,i)*e2(iv,i+1)-e4(iv,i)*e4(iv,i+1) |
|---|
| 145 | d(iv,j)=e1(iv,i+1)*e4(iv,i+1)-e2(iv,i+1)*e3(iv,i+1) |
|---|
| 146 | e(iv,j)=e2(iv,i+1)*(cah(iv,i+1)-cab(iv,i)) |
|---|
| 147 | & +e4(iv,i+1)*(cdb(iv,i)-cdh(iv,i+1)) |
|---|
| 148 | 99996 continue |
|---|
| 149 | 10003 continue |
|---|
| 150 | c |
|---|
| 151 | j=2*n |
|---|
| 152 | do 99995 iv=1,KDLON |
|---|
| 153 | a(iv,j)=e1(iv,n)-(1.E+0-emis(iv))*e3(iv,n) |
|---|
| 154 | b(iv,j)=e2(iv,n)-(1.E+0-emis(iv))*e4(iv,n) |
|---|
| 155 | d(iv,j)=0.E+0 |
|---|
| 156 | e(iv,j)=emis(iv)*pi*bsol(iv)-cab(iv,n)+(1.E+0-emis(iv))*cdb(iv,n) |
|---|
| 157 | 99995 continue |
|---|
| 158 | c....................................................................... |
|---|
| 159 | call sys3v(KDLON,2*n,a,b,d,e,y) |
|---|
| 160 | c....................................................................... |
|---|
| 161 | do 10004 i=1,n |
|---|
| 162 | do 99994 iv=1,KDLON |
|---|
| 163 | grg(iv,i)=grg(iv,i)*(y(iv,2*i-1)+y(iv,2*i)) |
|---|
| 164 | grh(iv,i)=grh(iv,i)*(y(iv,2*i-1)-y(iv,2*i)) |
|---|
| 165 | grj(iv,i)=grj(iv,i)*(y(iv,2*i-1)+y(iv,2*i)) |
|---|
| 166 | grk(iv,i)=grk(iv,i)*(y(iv,2*i-1)-y(iv,2*i)) |
|---|
| 167 | 99994 continue |
|---|
| 168 | 10004 continue |
|---|
| 169 | c....................................................................... |
|---|
| 170 | c les valeurs de flux "hemispheric constant" |
|---|
| 171 | c |
|---|
| 172 | IF (nsf.eq.0) THEN |
|---|
| 173 | do 10005 i=1,n |
|---|
| 174 | do 99993 iv=1,KDLON |
|---|
| 175 | fah(iv,i)=e3(iv,i)*y(iv,2*i-1)-e4(iv,i)*y(iv,2*i)+cah(iv,i) |
|---|
| 176 | fdh(iv,i)=e1(iv,i)*y(iv,2*i-1)-e2(iv,i)*y(iv,2*i)+cdh(iv,i) |
|---|
| 177 | 99993 continue |
|---|
| 178 | 10005 continue |
|---|
| 179 | do 99992 iv=1,KDLON |
|---|
| 180 | fah(iv,n+1)=e1(iv,n)*y(iv,2*n-1)+e2(iv,n)*y(iv,2*n)+cab(iv,n) |
|---|
| 181 | fdh(iv,n+1)=e3(iv,n)*y(iv,2*n-1)+e4(iv,n)*y(iv,2*n)+cdb(iv,n) |
|---|
| 182 | 99992 continue |
|---|
| 183 | GOTO 10100 |
|---|
| 184 | ENDIF |
|---|
| 185 | c....................................................................... |
|---|
| 186 | c passage a "source function" |
|---|
| 187 | c |
|---|
| 188 | c on applique une quadrature de dimension nq |
|---|
| 189 | c x est le vecteur des \mu de la quadrature |
|---|
| 190 | c w est le vecteur des poids correspondants |
|---|
| 191 | c nq est fixe en parameter |
|---|
| 192 | c x et w sont fixes en data |
|---|
| 193 | c |
|---|
| 194 | c....................................................................... |
|---|
| 195 | c on part d'en haut et on descent selon les nq angles pour calculer |
|---|
| 196 | c tous les flux descendants |
|---|
| 197 | c |
|---|
| 198 | do 10006 j=1,nq |
|---|
| 199 | do 99991 iv=1,KDLON |
|---|
| 200 | gri(iv,j)=0.E+0 |
|---|
| 201 | 99991 continue |
|---|
| 202 | 10006 continue |
|---|
| 203 | do 99990 iv=1,KDLON |
|---|
| 204 | fdh(iv,1)=0.E+0 |
|---|
| 205 | 99990 continue |
|---|
| 206 | do 10007 i=1,n |
|---|
| 207 | do 10008 j=1,nq |
|---|
| 208 | do 99989 iv=1,KDLON |
|---|
| 209 | gri(iv,j)=gri(iv,j)*exp(-tau(iv,i)/x(j)) |
|---|
| 210 | &+grj(iv,i)/(alambda(iv,i)*x(j)+1.E+0) |
|---|
| 211 | &*(1.E+0-exp(-tau(iv,i)*(alambda(iv,i)+1.E+0/x(j)))) |
|---|
| 212 | &+grk(iv,i)/(alambda(iv,i)*x(j)-1.E+0) |
|---|
| 213 | &*(exp(-tau(iv,i)/x(j))-exp(-tau(iv,i)*alambda(iv,i))) |
|---|
| 214 | &+sigma1(iv,i)*(1.E+0-exp(-tau(iv,i)/x(j))) |
|---|
| 215 | &+sigma2(iv,i)*(x(j)*exp(-tau(iv,i)/x(j))+tau(iv,i)-x(j)) |
|---|
| 216 | 99989 continue |
|---|
| 217 | 10008 continue |
|---|
| 218 | do 99988 iv=1,KDLON |
|---|
| 219 | fdh(iv,i+1)=0.E+0 |
|---|
| 220 | 99988 continue |
|---|
| 221 | do 10009 j=1,nq |
|---|
| 222 | do 99987 iv=1,KDLON |
|---|
| 223 | fdh(iv,i+1)=fdh(iv,i+1)+w(j)*x(j)*gri(iv,j) |
|---|
| 224 | 99987 continue |
|---|
| 225 | 10009 continue |
|---|
| 226 | 10007 continue |
|---|
| 227 | c....................................................................... |
|---|
| 228 | c on applique la condition de reflexion a sol |
|---|
| 229 | c |
|---|
| 230 | do 99986 iv=1,KDLON |
|---|
| 231 | fah(iv,n+1)=(1.E+0-emis(iv))*fdh(iv,n+1)+pi*emis(iv)*bsol(iv) |
|---|
| 232 | 99986 continue |
|---|
| 233 | do 10010 j=1,nq |
|---|
| 234 | do 99985 iv=1,KDLON |
|---|
| 235 | gri(iv,j)=2.E+0*fah(iv,n+1) |
|---|
| 236 | 99985 continue |
|---|
| 237 | 10010 continue |
|---|
| 238 | c....................................................................... |
|---|
| 239 | c on remonte pour calculer tous les flux ascendants |
|---|
| 240 | c |
|---|
| 241 | do 10011 i=n,1,-1 |
|---|
| 242 | do 10012 j=1,nq |
|---|
| 243 | do 99984 iv=1,KDLON |
|---|
| 244 | gri(iv,j)=gri(iv,j)*exp(-tau(iv,i)/x(j)) |
|---|
| 245 | &+grg(iv,i)/(alambda(iv,i)*x(j)-1.E+0) |
|---|
| 246 | &*(exp(-tau(iv,i)/x(j))-exp(-tau(iv,i)*alambda(iv,i))) |
|---|
| 247 | &+grh(iv,i)/(alambda(iv,i)*x(j)+1.E+0) |
|---|
| 248 | &*(1.E+0-exp(-tau(iv,i)*(alambda(iv,i)+1.E+0/x(j)))) |
|---|
| 249 | &+alpha1(iv,i)*(1.E+0-exp(-tau(iv,i)/x(j))) |
|---|
| 250 | &+alpha2(iv,i)*(x(j)-(tau(iv,i)+x(j))*exp(-tau(iv,i)/x(j))) |
|---|
| 251 | 99984 continue |
|---|
| 252 | 10012 continue |
|---|
| 253 | do 99983 iv=1,KDLON |
|---|
| 254 | fah(iv,i)=0.E+0 |
|---|
| 255 | 99983 continue |
|---|
| 256 | do 10013 j=1,nq |
|---|
| 257 | do 99982 iv=1,KDLON |
|---|
| 258 | fah(iv,i)=fah(iv,i)+w(j)*x(j)*gri(iv,j) |
|---|
| 259 | 99982 continue |
|---|
| 260 | 10013 continue |
|---|
| 261 | 10011 continue |
|---|
| 262 | c....................................................................... |
|---|
| 263 | 10100 continue |
|---|
| 264 | c....................................................................... |
|---|
| 265 | c |
|---|
| 266 | c....................................................................... |
|---|
| 267 | RETURN |
|---|
| 268 | END |
|---|
| 269 | |
|---|
| 270 | c *************************************************************** |
|---|
| 271 | |
|---|
| 272 | |
|---|
| 273 | SUBROUTINE sys3v(KDLON,n,a,b,d,e,y) |
|---|
| 274 | use dimradmars_mod, only: ndlon, ndlo2 |
|---|
| 275 | IMPLICIT NONE |
|---|
| 276 | c....................................................................... |
|---|
| 277 | c |
|---|
| 278 | c inversion d'un systeme lineaire tridiagonal |
|---|
| 279 | c |
|---|
| 280 | c | b1 d1 | | y1 | | e1 | |
|---|
| 281 | c | a2 b2 d2 | | y2 | | e2 | |
|---|
| 282 | c | a3 b3 d3 | * | y3 | = | e3 | |
|---|
| 283 | c | .... | | | | | |
|---|
| 284 | c | an bn | | yn | | en | |
|---|
| 285 | c |
|---|
| 286 | c in : * KDLON --> dimension de vectorisation |
|---|
| 287 | c * n --> dimension du systeme |
|---|
| 288 | c * a,b,d,e --> voir dessin |
|---|
| 289 | c |
|---|
| 290 | c out : * y --> voir dessin |
|---|
| 291 | c |
|---|
| 292 | c....................................................................... |
|---|
| 293 | c include des dimensions locales |
|---|
| 294 | c |
|---|
| 295 | #include "dimensions.h" |
|---|
| 296 | #include "dimphys.h" |
|---|
| 297 | !#include "dimradmars.h" |
|---|
| 298 | c....................................................................... |
|---|
| 299 | c declaration des arguments |
|---|
| 300 | c |
|---|
| 301 | INTEGER KDLON,n |
|---|
| 302 | REAL a(NDLO2,n),b(NDLO2,n),d(NDLO2,n),e(NDLO2,n),y(NDLO2,n) |
|---|
| 303 | c....................................................................... |
|---|
| 304 | c declaration des variables locales |
|---|
| 305 | c |
|---|
| 306 | INTEGER iv,i |
|---|
| 307 | REAL as(NDLON,4*nlayermx),ds(NDLON,4*nlayermx) |
|---|
| 308 | & ,x(NDLON,4*nlayermx) |
|---|
| 309 | c....................................................................... |
|---|
| 310 | c |
|---|
| 311 | c....................................................................... |
|---|
| 312 | do 99999 iv=1,KDLON |
|---|
| 313 | as(iv,n)=a(iv,n)/b(iv,n) |
|---|
| 314 | ds(iv,n)=e(iv,n)/b(iv,n) |
|---|
| 315 | 99999 continue |
|---|
| 316 | do 10001 i=n-1,1,-1 |
|---|
| 317 | do 99998 iv=1,KDLON |
|---|
| 318 | x(iv,i)=1.E+0/(b(iv,i)-d(iv,i)*as(iv,i+1)) |
|---|
| 319 | as(iv,i)=a(iv,i)*x(iv,i) |
|---|
| 320 | ds(iv,i)=(e(iv,i)-d(iv,i)*ds(iv,i+1))*x(iv,i) |
|---|
| 321 | 99998 continue |
|---|
| 322 | 10001 continue |
|---|
| 323 | do 99997 iv=1,KDLON |
|---|
| 324 | y(iv,1)=ds(iv,1) |
|---|
| 325 | 99997 continue |
|---|
| 326 | do 10002 i=2,n |
|---|
| 327 | do 99996 iv=1,KDLON |
|---|
| 328 | y(iv,i)=ds(iv,i)-as(iv,i)*y(iv,i-1) |
|---|
| 329 | 99996 continue |
|---|
| 330 | 10002 continue |
|---|
| 331 | c....................................................................... |
|---|
| 332 | c |
|---|
| 333 | c....................................................................... |
|---|
| 334 | RETURN |
|---|
| 335 | END |
|---|
