| 1 | SUBROUTINE FLOTT_GWD_RAN(NLON,NLEV,DTIME, pp, pn2, |
|---|
| 2 | I tt,uu,vv, |
|---|
| 3 | O zustr,zvstr, |
|---|
| 4 | O d_t, d_u, d_v) |
|---|
| 5 | |
|---|
| 6 | C ######################################################################## |
|---|
| 7 | C Parametrization of the momentum flux deposition due to a discrete |
|---|
| 8 | C numbe of gravity waves. |
|---|
| 9 | C F. Lott (2008) |
|---|
| 10 | C LMDz model stand alone version |
|---|
| 11 | C ######################################################################## |
|---|
| 12 | C |
|---|
| 13 | C |
|---|
| 14 | use dimphy |
|---|
| 15 | implicit none |
|---|
| 16 | |
|---|
| 17 | #include "dimensions.h" |
|---|
| 18 | #include "paramet.h" |
|---|
| 19 | |
|---|
| 20 | #include "YOEGWD.h" |
|---|
| 21 | #include "YOMCST.h" |
|---|
| 22 | |
|---|
| 23 | c VENUS ATTENTION: CP VARIABLE |
|---|
| 24 | real cp(NLON,NLEV) |
|---|
| 25 | |
|---|
| 26 | C DECLARATIONS: |
|---|
| 27 | C INPUTS |
|---|
| 28 | INTEGER NLON,NLEV |
|---|
| 29 | REAL DTIME |
|---|
| 30 | REAL pp(NLON,NLEV) ! Pressure at full levels |
|---|
| 31 | c VENUS ATTENTION: CP VARIABLE PN2 CALCULE EN AMONT DES PARAMETRISATIONS |
|---|
| 32 | REAL pn2(NLON,NLEV) ! N2 (BV^2) at 1/2 levels |
|---|
| 33 | REAL TT(NLON,NLEV) ! Temp at full levels |
|---|
| 34 | REAL UU(NLON,NLEV),VV(NLON,NLEV) ! Temp, and winds at full levels |
|---|
| 35 | |
|---|
| 36 | C OUTPUTS |
|---|
| 37 | REAL zustr(NLON),zvstr(NLON) ! Surface Stresses |
|---|
| 38 | REAL d_t(NLON,NLEV) ! Tendency on Temp. |
|---|
| 39 | REAL d_u(NLON,NLEV), d_v(NLON,NLEV) ! Tendencies on winds |
|---|
| 40 | |
|---|
| 41 | c ON CONSERVE LA MEMOIRE un certain temps AVEC UN SAVE ET UN COEF |
|---|
| 42 | real coef |
|---|
| 43 | parameter (coef = 0.986) |
|---|
| 44 | real,save,allocatable :: d_u_sav(:,:),d_v_sav(:,:) |
|---|
| 45 | LOGICAL firstcall |
|---|
| 46 | SAVE firstcall |
|---|
| 47 | DATA firstcall/.true./ |
|---|
| 48 | |
|---|
| 49 | C GRAVITY-WAVES SPECIFICATIONS |
|---|
| 50 | |
|---|
| 51 | INTEGER JK,NK,JP,NP,JO,NO,JW,NW,II,LL |
|---|
| 52 | PARAMETER(NK=4,NP=4,NO=4,NW=NK*NP*NO) ! nombres d'ondes envoyées |
|---|
| 53 | REAL KMIN,KMAX ! Min and Max horizontal wavenumbers |
|---|
| 54 | REAL OMIN,OMAX ! Max frequency and security for divisions! |
|---|
| 55 | PARAMETER (KMIN=2.E-5,KMAX=2.E-4,OMIN=2.E-5,OMAX=2.E-3) |
|---|
| 56 | REAL RUWMAX,SAT ! MAX EP FLUX AT LAUNCH LEVEL, SATURATION PARAMETER |
|---|
| 57 | PARAMETER (RUWMAX=1.E-3,SAT=0.1) |
|---|
| 58 | REAL ZK(NW,KLON),ZP(NW),ZO(NW,KLON) ! Waves absolute specifications! |
|---|
| 59 | REAL ZOM(NW,KLON),ZOP(NW,KLON) ! Intrisic frequency at 1/2 levels |
|---|
| 60 | REAL PHM(NW,KLON),PHP(NW,KLON) ! Phi at 1/2 levels |
|---|
| 61 | REAL RUW0(NW,KLON) ! Fluxes at lauch level |
|---|
| 62 | REAL RUWP(NW,KLON),RVWP(NW,KLON)! Fluxes X and Y for each wave |
|---|
| 63 | INTEGER LAUNCH |
|---|
| 64 | REAL ALEAS |
|---|
| 65 | EXTERNAL ALEAS |
|---|
| 66 | |
|---|
| 67 | C OUTPUT FOR RAW PROFILES: |
|---|
| 68 | |
|---|
| 69 | INTEGER IU |
|---|
| 70 | |
|---|
| 71 | C PARAMETERS OF WAVES DISSPATIONS |
|---|
| 72 | REAL RDISS,ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ |
|---|
| 73 | PARAMETER (RDISS=10.,ZOISEC=1.E-5) |
|---|
| 74 | |
|---|
| 75 | C INTERNAL ARRAYS |
|---|
| 76 | real PR,RH(KLON,KLEV+1) ! Ref Press. intervient dans amplitude ? |
|---|
| 77 | PARAMETER (PR=9.2e6) ! VENUS!! |
|---|
| 78 | |
|---|
| 79 | |
|---|
| 80 | REAL UH(KLON,KLEV+1),VH(KLON,KLEV+1) ! Winds at 1/2 levels |
|---|
| 81 | REAL PH(KLON,KLEV+1) ! Pressure at 1/2 levels |
|---|
| 82 | REAL PHM1(KLON,KLEV+1) ! 1/Press at 1/2 levels |
|---|
| 83 | REAL BV(KLON,KLEV+1),BVSEC ! BVF at 1/2 levels |
|---|
| 84 | REAL PSEC ! SECURITY TO AVOID DIVISION |
|---|
| 85 | PARAMETER (PSEC=1.E-6) ! BY O PRESSURE VALUE |
|---|
| 86 | REAL RUW(KLON,KLEV+1) ! Flux x at semi levels |
|---|
| 87 | REAL RVW(KLON,KLEV+1) ! Flux y at semi levels |
|---|
| 88 | |
|---|
| 89 | C INITIALISATION |
|---|
| 90 | |
|---|
| 91 | IF (firstcall) THEN |
|---|
| 92 | allocate(d_u_sav(klon,klev),d_v_sav(klon,klev)) |
|---|
| 93 | firstcall=.false. |
|---|
| 94 | ENDIF |
|---|
| 95 | |
|---|
| 96 | BVSEC=1.E-6 |
|---|
| 97 | |
|---|
| 98 | c PRINT *,ALEAS(0.) |
|---|
| 99 | |
|---|
| 100 | C WAVES CHARACTERISTTICS |
|---|
| 101 | |
|---|
| 102 | JW=0 |
|---|
| 103 | DO 11 JP=1,NP |
|---|
| 104 | DO 11 JK=1,NK |
|---|
| 105 | DO 11 JO=1,NO |
|---|
| 106 | JW=JW+1 |
|---|
| 107 | ZP(JW)=2.*RPI*FLOAT(JP-1)/FLOAT(NP) ! Angle |
|---|
| 108 | DO 11 II=1,NLON |
|---|
| 109 | ZK(JW,II)=KMIN+(KMAX-KMIN)*ALEAS(0.) ! Hor. Wavenbr |
|---|
| 110 | ZO(JW,II)=OMIN+(OMAX-OMIN)*ALEAS(0.) ! Absolute frequency |
|---|
| 111 | RUW0(JW,II)=RUWMAX/FLOAT(NW)*ALEAS(0.) ! Momentum flux at launch lev |
|---|
| 112 | 11 CONTINUE |
|---|
| 113 | |
|---|
| 114 | c print*,"cos(ZP)=",cos(ZP) |
|---|
| 115 | |
|---|
| 116 | c PRINT *,'Catch 22' |
|---|
| 117 | |
|---|
| 118 | C 2. SETUP |
|---|
| 119 | |
|---|
| 120 | c VENUS: CP(T) at ph levels |
|---|
| 121 | do II=1,NLON |
|---|
| 122 | do LL=2,KLEV |
|---|
| 123 | cp(II,LL)=cpdet(0.5*(TT(II,LL)+TT(II,LL-1))) |
|---|
| 124 | enddo |
|---|
| 125 | enddo |
|---|
| 126 | C PRESSURE AT HALPH LEVELS |
|---|
| 127 | |
|---|
| 128 | II=KLON/2 |
|---|
| 129 | |
|---|
| 130 | DO 20 LL=2,KLEV |
|---|
| 131 | PH(:,LL)=EXP((ALOG(PP(:,LL))+ALOG(PP(:,LL-1)))/2.) |
|---|
| 132 | PHM1(:,LL)=1./PH(:,LL) |
|---|
| 133 | 20 CONTINUE |
|---|
| 134 | |
|---|
| 135 | PH(:,KLEV+1)=0. |
|---|
| 136 | PHM1(:,KLEV+1)=1./PSEC |
|---|
| 137 | |
|---|
| 138 | PH(:,1)=2.*PP(:,1)-PH(:,2) |
|---|
| 139 | |
|---|
| 140 | C LOG-PRESSURE ALTITUDE AND DENSITY |
|---|
| 141 | |
|---|
| 142 | |
|---|
| 143 | DO 22 LL=1,KLEV+1 |
|---|
| 144 | RH(:,LL)=PH(:,LL)/PR ! Reference density |
|---|
| 145 | 22 CONTINUE |
|---|
| 146 | |
|---|
| 147 | C LAUNCHING ALTITUDE |
|---|
| 148 | |
|---|
| 149 | DO LL=1,NLEV |
|---|
| 150 | IF(PH(NLON/2,LL)/PH(NLON/2,1).GT.0.8)LAUNCH=LL |
|---|
| 151 | ENDDO |
|---|
| 152 | c PRINT *,'LAUNCH:',LAUNCH |
|---|
| 153 | |
|---|
| 154 | C WINDS AND BV FREQUENCY AT 1/2 LEVELS |
|---|
| 155 | |
|---|
| 156 | DO 23 LL=2,KLEV |
|---|
| 157 | UH(:,LL)=0.5*(UU(:,LL)+UU(:,LL-1)) |
|---|
| 158 | VH(:,LL)=0.5*(VV(:,LL)+VV(:,LL-1)) |
|---|
| 159 | |
|---|
| 160 | c VENUS ATTENTION: CP VARIABLE PSTAB CALCULE EN AMONT DES PARAMETRISATIONS |
|---|
| 161 | BV(:,LL)=SQRT(AMAX1(BVSEC,pn2(:,LL))) |
|---|
| 162 | |
|---|
| 163 | 23 CONTINUE |
|---|
| 164 | |
|---|
| 165 | c print *,'catch 23' |
|---|
| 166 | |
|---|
| 167 | BV(:,1)=BV(:,2) |
|---|
| 168 | UH(:,1)=0. |
|---|
| 169 | VH(:,1)=0. |
|---|
| 170 | BV(:,KLEV+1)=BV(:,KLEV) |
|---|
| 171 | UH(:,KLEV+1)=UU(:,KLEV) |
|---|
| 172 | VH(:,KLEV+1)=VV(:,KLEV) |
|---|
| 173 | |
|---|
| 174 | C 3. COMPUTE THE FLUXES!! |
|---|
| 175 | |
|---|
| 176 | DO 3 JW=1,NW |
|---|
| 177 | |
|---|
| 178 | ZOP(JW,:)=ZO(JW,:)-ZK(JW,:)*COS(ZP(JW))*UH(:,LAUNCH) ! Intrinsic |
|---|
| 179 | C -ZK(JW,:)*SIN(ZP(JW))*VH(:,LAUNCH) ! Frequency |
|---|
| 180 | PHP(JW,:)= |
|---|
| 181 | C AMIN1( ! Phi at launch level |
|---|
| 182 | C ENSURES THE IMPOSED MOM FLUX: |
|---|
| 183 | C SQRT(ABS(ZOP(JW,:))*BV(:,LAUNCH) |
|---|
| 184 | C *RUW0(JW,:)/RH(:,LAUNCH))/ZK(JW,:) |
|---|
| 185 | RUWP(JW,:)=COS(ZP(JW))*SIGN(1.,ZOP(JW,:))*RUW0(JW,:) |
|---|
| 186 | RVWP(JW,:)=SIN(ZP(JW))*SIGN(1.,ZOP(JW,:))*RUW0(JW,:) |
|---|
| 187 | |
|---|
| 188 | 3 CONTINUE |
|---|
| 189 | |
|---|
| 190 | do LL=1,LAUNCH |
|---|
| 191 | RUW(:,LL)=0 |
|---|
| 192 | RVW(:,LL)=0 |
|---|
| 193 | do JW=1,NW |
|---|
| 194 | RUW(:,LL)=RUW(:,LL)+RUWP(JW,:) |
|---|
| 195 | RVW(:,LL)=RVW(:,LL)+RVWP(JW,:) |
|---|
| 196 | enddo |
|---|
| 197 | enddo |
|---|
| 198 | |
|---|
| 199 | DO 33 LL=LAUNCH,KLEV-1 |
|---|
| 200 | |
|---|
| 201 | DO 34 JW=1,NW |
|---|
| 202 | ZOM(JW,:)=ZOP(JW,:) |
|---|
| 203 | PHM(JW,:)=PHP(JW,:) |
|---|
| 204 | |
|---|
| 205 | ZOP(JW,:)=ZO(JW,:)-ZK(JW,:)*COS(ZP(JW))*UH(:,LL+1) ! Intrinsic |
|---|
| 206 | C -ZK(JW,:)*SIN(ZP(JW))*VH(:,LL+1) ! Frequency |
|---|
| 207 | C |
|---|
| 208 | PHP(JW,:)= |
|---|
| 209 | C AMIN1( |
|---|
| 210 | C NO BREAKING |
|---|
| 211 | C SQRT(PH(:,LL)*PHM1(:,LL+1))*PHM(JW,:) |
|---|
| 212 | C *SQRT(BV(:,LL+1)/BV(:,LL) |
|---|
| 213 | C *ABS(ZOP(JW,:))/AMAX1(ABS(ZOM(JW,:)),ZOISEC)) |
|---|
| 214 | C CRITICAL LEVELS |
|---|
| 215 | C , |
|---|
| 216 | C AMAX1(0.,SIGN(1.,ZOP(JW,:)*ZOM(JW,:))) |
|---|
| 217 | C SATURATION |
|---|
| 218 | C *ZOP(JW,:)**2/ZK(JW,:)**2/SQRT(FLOAT(NK*NO))*SAT |
|---|
| 219 | C ) |
|---|
| 220 | |
|---|
| 221 | 34 CONTINUE |
|---|
| 222 | |
|---|
| 223 | DO 35 JW=1,NW |
|---|
| 224 | RUWP(JW,:)=ZOP(JW,:)*COS(ZP(JW))*ZK(JW,:)**2 |
|---|
| 225 | C *PH(:,LL+1)/PR*PHP(JW,:)**2/BV(:,LL+1) |
|---|
| 226 | C /AMAX1(ABS(ZOP(JW,:)),ZOISEC)**2 |
|---|
| 227 | RVWP(JW,:)=ZOP(JW,:)*SIN(ZP(JW))*ZK(JW,:)**2 |
|---|
| 228 | C *PH(:,LL+1)/PR*PHP(JW,:)**2/BV(:,LL+1) |
|---|
| 229 | C /AMAX1(ABS(ZOP(JW,:)),ZOISEC)**2 |
|---|
| 230 | 35 CONTINUE |
|---|
| 231 | |
|---|
| 232 | RUW(:,LL+1)=0. |
|---|
| 233 | RVW(:,LL+1)=0. |
|---|
| 234 | |
|---|
| 235 | |
|---|
| 236 | DO 36 JW=1,NW |
|---|
| 237 | RUW(:,LL+1)=RUW(:,LL+1)+RUWP(JW,:) |
|---|
| 238 | RVW(:,LL+1)=RVW(:,LL+1)+RVWP(JW,:) |
|---|
| 239 | 36 CONTINUE |
|---|
| 240 | |
|---|
| 241 | |
|---|
| 242 | 33 CONTINUE |
|---|
| 243 | |
|---|
| 244 | C 4 CALCUL DES TENDANCES: |
|---|
| 245 | |
|---|
| 246 | C RECTIFICATION AU SOMMET ET DANS LES BASSES COUCHES: |
|---|
| 247 | |
|---|
| 248 | RUW(:,1)=RUW(:,LAUNCH) |
|---|
| 249 | RVW(:,1)=RVW(:,LAUNCH) |
|---|
| 250 | DO 4 LL=2,LAUNCH |
|---|
| 251 | RUW(:,LL)=RUW(:,LL-1)+(RUW(:,LAUNCH+1)-RUW(:,1))* |
|---|
| 252 | C (PH(:,LL)-PH(:,LL-1))/(PH(:,LAUNCH+1)-PH(:,1)) |
|---|
| 253 | RVW(:,LL)=RVW(:,LL-1)+(RVW(:,LAUNCH+1)-RVW(:,1))* |
|---|
| 254 | C (PH(:,LL)-PH(:,LL-1))/(PH(:,LAUNCH+1)-PH(:,1)) |
|---|
| 255 | 4 CONTINUE |
|---|
| 256 | |
|---|
| 257 | DO 41 LL=1,KLEV |
|---|
| 258 | D_U(:,LL)=RG*(RUW(:,LL+1)-RUW(:,LL)) |
|---|
| 259 | C /(PH(:,LL+1)-PH(:,LL))*DTIME |
|---|
| 260 | D_V(:,LL)=RG*(RVW(:,LL+1)-RVW(:,LL)) |
|---|
| 261 | C /(PH(:,LL+1)-PH(:,LL))*DTIME |
|---|
| 262 | 41 CONTINUE |
|---|
| 263 | |
|---|
| 264 | c ON CONSERVE LA MEMOIRE un certain temps AVEC UN SAVE ET UN COEF |
|---|
| 265 | d_u = d_u + coef * d_u_sav |
|---|
| 266 | d_v = d_v + coef * d_v_sav |
|---|
| 267 | d_u_sav = d_u |
|---|
| 268 | d_v_sav = d_v |
|---|
| 269 | |
|---|
| 270 | c PRINT *,'OUT LOTT, LONG?' |
|---|
| 271 | c DO JW=1,NW |
|---|
| 272 | c PRINT *,'ONDE ',COS(ZP(JW))*JW,' KM:',2*RPI/ZK(JW)/1000. |
|---|
| 273 | c C ,' HR:',2*RPI/ZO(JW)/3600.,' C:',ZO(JW)/ZK(JW) |
|---|
| 274 | c ENDDO |
|---|
| 275 | |
|---|
| 276 | RETURN |
|---|
| 277 | END |
|---|
| 278 | |
|---|
| 279 | |
|---|
| 280 | c=================================================================== |
|---|
| 281 | c=================================================================== |
|---|
| 282 | c=================================================================== |
|---|
| 283 | c=================================================================== |
|---|
| 284 | |
|---|
| 285 | FUNCTION ALEAS (R) |
|---|
| 286 | C***BEGIN PROLOGUE ALEAS |
|---|
| 287 | C***PURPOSE Generate a uniformly distributed random number. |
|---|
| 288 | C***LIBRARY SLATEC (FNLIB) |
|---|
| 289 | C***CATEGORY L6A21 |
|---|
| 290 | C***TYPE SINGLE PRECISION (ALEAS-S) |
|---|
| 291 | C***KEYWORDS FNLIB, ALEAS NUMBER, SPECIAL FUNCTIONS, UNIFORM |
|---|
| 292 | C***AUTHOR Fullerton, W., (LANL) |
|---|
| 293 | C***DESCRIPTION |
|---|
| 294 | C |
|---|
| 295 | C This pseudo-random number generator is portable among a wide |
|---|
| 296 | C variety of computers. RAND(R) undoubtedly is not as good as many |
|---|
| 297 | C readily available installation dependent versions, and so this |
|---|
| 298 | C routine is not recommended for widespread usage. Its redeeming |
|---|
| 299 | C feature is that the exact same random numbers (to within final round- |
|---|
| 300 | C off error) can be generated from machine to machine. Thus, programs |
|---|
| 301 | C that make use of random numbers can be easily transported to and |
|---|
| 302 | C checked in a new environment. |
|---|
| 303 | C |
|---|
| 304 | C The random numbers are generated by the linear congruential |
|---|
| 305 | C method described, e.g., by Knuth in Seminumerical Methods (p.9), |
|---|
| 306 | C Addison-Wesley, 1969. Given the I-th number of a pseudo-random |
|---|
| 307 | C sequence, the I+1 -st number is generated from |
|---|
| 308 | C X(I+1) = (A*X(I) + C) MOD M, |
|---|
| 309 | C where here M = 2**22 = 4194304, C = 1731 and several suitable values |
|---|
| 310 | C of the multiplier A are discussed below. Both the multiplier A and |
|---|
| 311 | C random number X are represented in double precision as two 11-bit |
|---|
| 312 | C words. The constants are chosen so that the period is the maximum |
|---|
| 313 | C possible, 4194304. |
|---|
| 314 | C |
|---|
| 315 | C In order that the same numbers be generated from machine to |
|---|
| 316 | C machine, it is necessary that 23-bit integers be reducible modulo |
|---|
| 317 | C 2**11 exactly, that 23-bit integers be added exactly, and that 11-bit |
|---|
| 318 | C integers be multiplied exactly. Furthermore, if the restart option |
|---|
| 319 | C is used (where R is between 0 and 1), then the product R*2**22 = |
|---|
| 320 | C R*4194304 must be correct to the nearest integer. |
|---|
| 321 | C |
|---|
| 322 | C The first four random numbers should be .0004127026, |
|---|
| 323 | C .6750836372, .1614754200, and .9086198807. The tenth random number |
|---|
| 324 | C is .5527787209, and the hundredth is .3600893021 . The thousandth |
|---|
| 325 | C number should be .2176990509 . |
|---|
| 326 | C |
|---|
| 327 | C In order to generate several effectively independent sequences |
|---|
| 328 | C with the same generator, it is necessary to know the random number |
|---|
| 329 | C for several widely spaced calls. The I-th random number times 2**22, |
|---|
| 330 | C where I=K*P/8 and P is the period of the sequence (P = 2**22), is |
|---|
| 331 | C still of the form L*P/8. In particular we find the I-th random |
|---|
| 332 | C number multiplied by 2**22 is given by |
|---|
| 333 | C I = 0 1*P/8 2*P/8 3*P/8 4*P/8 5*P/8 6*P/8 7*P/8 8*P/8 |
|---|
| 334 | C RAND= 0 5*P/8 2*P/8 7*P/8 4*P/8 1*P/8 6*P/8 3*P/8 0 |
|---|
| 335 | C Thus the 4*P/8 = 2097152 random number is 2097152/2**22. |
|---|
| 336 | C |
|---|
| 337 | C Several multipliers have been subjected to the spectral test |
|---|
| 338 | C (see Knuth, p. 82). Four suitable multipliers roughly in order of |
|---|
| 339 | C goodness according to the spectral test are |
|---|
| 340 | C 3146757 = 1536*2048 + 1029 = 2**21 + 2**20 + 2**10 + 5 |
|---|
| 341 | C 2098181 = 1024*2048 + 1029 = 2**21 + 2**10 + 5 |
|---|
| 342 | C 3146245 = 1536*2048 + 517 = 2**21 + 2**20 + 2**9 + 5 |
|---|
| 343 | C 2776669 = 1355*2048 + 1629 = 5**9 + 7**7 + 1 |
|---|
| 344 | C |
|---|
| 345 | C In the table below LOG10(NU(I)) gives roughly the number of |
|---|
| 346 | C random decimal digits in the random numbers considered I at a time. |
|---|
| 347 | C C is the primary measure of goodness. In both cases bigger is better. |
|---|
| 348 | C |
|---|
| 349 | C LOG10 NU(I) C(I) |
|---|
| 350 | C A I=2 I=3 I=4 I=5 I=2 I=3 I=4 I=5 |
|---|
| 351 | C |
|---|
| 352 | C 3146757 3.3 2.0 1.6 1.3 3.1 1.3 4.6 2.6 |
|---|
| 353 | C 2098181 3.3 2.0 1.6 1.2 3.2 1.3 4.6 1.7 |
|---|
| 354 | C 3146245 3.3 2.2 1.5 1.1 3.2 4.2 1.1 0.4 |
|---|
| 355 | C 2776669 3.3 2.1 1.6 1.3 2.5 2.0 1.9 2.6 |
|---|
| 356 | C Best |
|---|
| 357 | C Possible 3.3 2.3 1.7 1.4 3.6 5.9 9.7 14.9 |
|---|
| 358 | C |
|---|
| 359 | C Input Argument -- |
|---|
| 360 | C R If R=0., the next random number of the sequence is generated. |
|---|
| 361 | C If R .LT. 0., the last generated number will be returned for |
|---|
| 362 | C possible use in a restart procedure. |
|---|
| 363 | C If R .GT. 0., the sequence of random numbers will start with |
|---|
| 364 | C the seed R mod 1. This seed is also returned as the value of |
|---|
| 365 | C RAND provided the arithmetic is done exactly. |
|---|
| 366 | C |
|---|
| 367 | C Output Value -- |
|---|
| 368 | C RAND a pseudo-random number between 0. and 1. |
|---|
| 369 | C |
|---|
| 370 | C***REFERENCES (NONE) |
|---|
| 371 | C***ROUTINES CALLED (NONE) |
|---|
| 372 | C***REVISION HISTORY (YYMMDD) |
|---|
| 373 | C 770401 DATE WRITTEN |
|---|
| 374 | C 890531 Changed all specific intrinsics to generic. (WRB) |
|---|
| 375 | C 890531 REVISION DATE from Version 3.2 |
|---|
| 376 | C 891214 Prologue converted to Version 4.0 format. (BAB) |
|---|
| 377 | C***END PROLOGUE RAND |
|---|
| 378 | SAVE IA1, IA0, IA1MA0, IC, IX1, IX0 |
|---|
| 379 | DATA IA1, IA0, IA1MA0 /1536, 1029, 507/ |
|---|
| 380 | DATA IC /1731/ |
|---|
| 381 | DATA IX1, IX0 /0, 0/ |
|---|
| 382 | C***FIRST EXECUTABLE STATEMENT RAND |
|---|
| 383 | IF (R.LT.0.) GO TO 10 |
|---|
| 384 | IF (R.GT.0.) GO TO 20 |
|---|
| 385 | C |
|---|
| 386 | C A*X = 2**22*IA1*IX1 + 2**11*(IA1*IX1 + (IA1-IA0)*(IX0-IX1) |
|---|
| 387 | C + IA0*IX0) + IA0*IX0 |
|---|
| 388 | C |
|---|
| 389 | IY0 = IA0*IX0 |
|---|
| 390 | IY1 = IA1*IX1 + IA1MA0*(IX0-IX1) + IY0 |
|---|
| 391 | IY0 = IY0 + IC |
|---|
| 392 | IX0 = MOD (IY0, 2048) |
|---|
| 393 | IY1 = IY1 + (IY0-IX0)/2048 |
|---|
| 394 | IX1 = MOD (IY1, 2048) |
|---|
| 395 | C |
|---|
| 396 | 10 ALEAS = IX1*2048 + IX0 |
|---|
| 397 | ALEAS = ALEAS / 4194304. |
|---|
| 398 | RETURN |
|---|
| 399 | C |
|---|
| 400 | 20 IX1 = MOD(R,1.)*4194304. + 0.5 |
|---|
| 401 | IX0 = MOD (IX1, 2048) |
|---|
| 402 | IX1 = (IX1-IX0)/2048 |
|---|
| 403 | GO TO 10 |
|---|
| 404 | C |
|---|
| 405 | END |
|---|
| 406 | |
|---|
| 407 | |
|---|
