SUBROUTINE FLOTT_GWD_RAN(NLON,NLEV,DTIME, pp, pn2, I tt,uu,vv, O zustr,zvstr, O d_t, d_u, d_v) C ######################################################################## C Parametrization of the momentum flux deposition due to a discrete C numbe of gravity waves. C F. Lott (2008) C LMDz model stand alone version C ######################################################################## C C use dimphy implicit none #include "dimensions.h" #include "paramet.h" #include "YOEGWD.h" #include "YOMCST.h" c VENUS ATTENTION: CP VARIABLE real cp(NLON,NLEV) C DECLARATIONS: C INPUTS INTEGER NLON,NLEV REAL DTIME REAL pp(NLON,NLEV) ! Pressure at full levels c VENUS ATTENTION: CP VARIABLE PN2 CALCULE EN AMONT DES PARAMETRISATIONS REAL pn2(NLON,NLEV) ! N2 (BV^2) at 1/2 levels REAL TT(NLON,NLEV) ! Temp at full levels REAL UU(NLON,NLEV),VV(NLON,NLEV) ! Temp, and winds at full levels C OUTPUTS REAL zustr(NLON),zvstr(NLON) ! Surface Stresses REAL d_t(NLON,NLEV) ! Tendency on Temp. REAL d_u(NLON,NLEV), d_v(NLON,NLEV) ! Tendencies on winds c ON CONSERVE LA MEMOIRE un certain temps AVEC UN SAVE ET UN COEF real coef parameter (coef = 0.986) real,save,allocatable :: d_u_sav(:,:),d_v_sav(:,:) LOGICAL firstcall SAVE firstcall DATA firstcall/.true./ C GRAVITY-WAVES SPECIFICATIONS INTEGER JK,NK,JP,NP,JO,NO,JW,NW,II,LL PARAMETER(NK=4,NP=4,NO=4,NW=NK*NP*NO) ! nombres d'ondes envoyées REAL KMIN,KMAX ! Min and Max horizontal wavenumbers REAL OMIN,OMAX ! Max frequency and security for divisions! PARAMETER (KMIN=2.E-5,KMAX=2.E-4,OMIN=2.E-5,OMAX=2.E-3) REAL RUWMAX,SAT ! MAX EP FLUX AT LAUNCH LEVEL, SATURATION PARAMETER PARAMETER (RUWMAX=1.E-3,SAT=0.1) REAL ZK(NW,KLON),ZP(NW),ZO(NW,KLON) ! Waves absolute specifications! REAL ZOM(NW,KLON),ZOP(NW,KLON) ! Intrisic frequency at 1/2 levels REAL PHM(NW,KLON),PHP(NW,KLON) ! Phi at 1/2 levels REAL RUW0(NW,KLON) ! Fluxes at lauch level REAL RUWP(NW,KLON),RVWP(NW,KLON)! Fluxes X and Y for each wave INTEGER LAUNCH REAL ALEAS EXTERNAL ALEAS C OUTPUT FOR RAW PROFILES: INTEGER IU C PARAMETERS OF WAVES DISSPATIONS REAL RDISS,ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ PARAMETER (RDISS=10.,ZOISEC=1.E-5) C INTERNAL ARRAYS real PR,RH(KLON,KLEV+1) ! Ref Press. intervient dans amplitude ? PARAMETER (PR=9.2e6) ! VENUS!! REAL UH(KLON,KLEV+1),VH(KLON,KLEV+1) ! Winds at 1/2 levels REAL PH(KLON,KLEV+1) ! Pressure at 1/2 levels REAL PHM1(KLON,KLEV+1) ! 1/Press at 1/2 levels REAL BV(KLON,KLEV+1),BVSEC ! BVF at 1/2 levels REAL PSEC ! SECURITY TO AVOID DIVISION PARAMETER (PSEC=1.E-6) ! BY O PRESSURE VALUE REAL RUW(KLON,KLEV+1) ! Flux x at semi levels REAL RVW(KLON,KLEV+1) ! Flux y at semi levels C INITIALISATION IF (firstcall) THEN allocate(d_u_sav(klon,klev),d_v_sav(klon,klev)) firstcall=.false. ENDIF BVSEC=1.E-6 c PRINT *,ALEAS(0.) C WAVES CHARACTERISTTICS JW=0 DO 11 JP=1,NP DO 11 JK=1,NK DO 11 JO=1,NO JW=JW+1 ZP(JW)=2.*RPI*FLOAT(JP-1)/FLOAT(NP) ! Angle DO 11 II=1,NLON ZK(JW,II)=KMIN+(KMAX-KMIN)*ALEAS(0.) ! Hor. Wavenbr ZO(JW,II)=OMIN+(OMAX-OMIN)*ALEAS(0.) ! Absolute frequency RUW0(JW,II)=RUWMAX/FLOAT(NW)*ALEAS(0.) ! Momentum flux at launch lev 11 CONTINUE c print*,"cos(ZP)=",cos(ZP) c PRINT *,'Catch 22' C 2. SETUP c VENUS: CP(T) at ph levels do II=1,NLON do LL=2,KLEV cp(II,LL)=cpdet(0.5*(TT(II,LL)+TT(II,LL-1))) enddo enddo C PRESSURE AT HALPH LEVELS II=KLON/2 DO 20 LL=2,KLEV PH(:,LL)=EXP((ALOG(PP(:,LL))+ALOG(PP(:,LL-1)))/2.) PHM1(:,LL)=1./PH(:,LL) 20 CONTINUE PH(:,KLEV+1)=0. PHM1(:,KLEV+1)=1./PSEC PH(:,1)=2.*PP(:,1)-PH(:,2) C LOG-PRESSURE ALTITUDE AND DENSITY DO 22 LL=1,KLEV+1 RH(:,LL)=PH(:,LL)/PR ! Reference density 22 CONTINUE C LAUNCHING ALTITUDE DO LL=1,NLEV IF(PH(NLON/2,LL)/PH(NLON/2,1).GT.0.8)LAUNCH=LL ENDDO c PRINT *,'LAUNCH:',LAUNCH C WINDS AND BV FREQUENCY AT 1/2 LEVELS DO 23 LL=2,KLEV UH(:,LL)=0.5*(UU(:,LL)+UU(:,LL-1)) VH(:,LL)=0.5*(VV(:,LL)+VV(:,LL-1)) c VENUS ATTENTION: CP VARIABLE PSTAB CALCULE EN AMONT DES PARAMETRISATIONS BV(:,LL)=SQRT(AMAX1(BVSEC,pn2(:,LL))) 23 CONTINUE c print *,'catch 23' BV(:,1)=BV(:,2) UH(:,1)=0. VH(:,1)=0. BV(:,KLEV+1)=BV(:,KLEV) UH(:,KLEV+1)=UU(:,KLEV) VH(:,KLEV+1)=VV(:,KLEV) C 3. COMPUTE THE FLUXES!! DO 3 JW=1,NW ZOP(JW,:)=ZO(JW,:)-ZK(JW,:)*COS(ZP(JW))*UH(:,LAUNCH) ! Intrinsic C -ZK(JW,:)*SIN(ZP(JW))*VH(:,LAUNCH) ! Frequency PHP(JW,:)= C AMIN1( ! Phi at launch level C ENSURES THE IMPOSED MOM FLUX: C SQRT(ABS(ZOP(JW,:))*BV(:,LAUNCH) C *RUW0(JW,:)/RH(:,LAUNCH))/ZK(JW,:) RUWP(JW,:)=COS(ZP(JW))*SIGN(1.,ZOP(JW,:))*RUW0(JW,:) RVWP(JW,:)=SIN(ZP(JW))*SIGN(1.,ZOP(JW,:))*RUW0(JW,:) 3 CONTINUE do LL=1,LAUNCH RUW(:,LL)=0 RVW(:,LL)=0 do JW=1,NW RUW(:,LL)=RUW(:,LL)+RUWP(JW,:) RVW(:,LL)=RVW(:,LL)+RVWP(JW,:) enddo enddo DO 33 LL=LAUNCH,KLEV-1 DO 34 JW=1,NW ZOM(JW,:)=ZOP(JW,:) PHM(JW,:)=PHP(JW,:) ZOP(JW,:)=ZO(JW,:)-ZK(JW,:)*COS(ZP(JW))*UH(:,LL+1) ! Intrinsic C -ZK(JW,:)*SIN(ZP(JW))*VH(:,LL+1) ! Frequency C PHP(JW,:)= C AMIN1( C NO BREAKING C SQRT(PH(:,LL)*PHM1(:,LL+1))*PHM(JW,:) C *SQRT(BV(:,LL+1)/BV(:,LL) C *ABS(ZOP(JW,:))/AMAX1(ABS(ZOM(JW,:)),ZOISEC)) C CRITICAL LEVELS C , C AMAX1(0.,SIGN(1.,ZOP(JW,:)*ZOM(JW,:))) C SATURATION C *ZOP(JW,:)**2/ZK(JW,:)**2/SQRT(FLOAT(NK*NO))*SAT C ) 34 CONTINUE DO 35 JW=1,NW RUWP(JW,:)=ZOP(JW,:)*COS(ZP(JW))*ZK(JW,:)**2 C *PH(:,LL+1)/PR*PHP(JW,:)**2/BV(:,LL+1) C /AMAX1(ABS(ZOP(JW,:)),ZOISEC)**2 RVWP(JW,:)=ZOP(JW,:)*SIN(ZP(JW))*ZK(JW,:)**2 C *PH(:,LL+1)/PR*PHP(JW,:)**2/BV(:,LL+1) C /AMAX1(ABS(ZOP(JW,:)),ZOISEC)**2 35 CONTINUE RUW(:,LL+1)=0. RVW(:,LL+1)=0. DO 36 JW=1,NW RUW(:,LL+1)=RUW(:,LL+1)+RUWP(JW,:) RVW(:,LL+1)=RVW(:,LL+1)+RVWP(JW,:) 36 CONTINUE 33 CONTINUE C 4 CALCUL DES TENDANCES: C RECTIFICATION AU SOMMET ET DANS LES BASSES COUCHES: RUW(:,1)=RUW(:,LAUNCH) RVW(:,1)=RVW(:,LAUNCH) DO 4 LL=2,LAUNCH RUW(:,LL)=RUW(:,LL-1)+(RUW(:,LAUNCH+1)-RUW(:,1))* C (PH(:,LL)-PH(:,LL-1))/(PH(:,LAUNCH+1)-PH(:,1)) RVW(:,LL)=RVW(:,LL-1)+(RVW(:,LAUNCH+1)-RVW(:,1))* C (PH(:,LL)-PH(:,LL-1))/(PH(:,LAUNCH+1)-PH(:,1)) 4 CONTINUE DO 41 LL=1,KLEV D_U(:,LL)=RG*(RUW(:,LL+1)-RUW(:,LL)) C /(PH(:,LL+1)-PH(:,LL))*DTIME D_V(:,LL)=RG*(RVW(:,LL+1)-RVW(:,LL)) C /(PH(:,LL+1)-PH(:,LL))*DTIME 41 CONTINUE c ON CONSERVE LA MEMOIRE un certain temps AVEC UN SAVE ET UN COEF d_u = d_u + coef * d_u_sav d_v = d_v + coef * d_v_sav d_u_sav = d_u d_v_sav = d_v c PRINT *,'OUT LOTT, LONG?' c DO JW=1,NW c PRINT *,'ONDE ',COS(ZP(JW))*JW,' KM:',2*RPI/ZK(JW)/1000. c C ,' HR:',2*RPI/ZO(JW)/3600.,' C:',ZO(JW)/ZK(JW) c ENDDO RETURN END c=================================================================== c=================================================================== c=================================================================== c=================================================================== FUNCTION ALEAS (R) C***BEGIN PROLOGUE ALEAS C***PURPOSE Generate a uniformly distributed random number. C***LIBRARY SLATEC (FNLIB) C***CATEGORY L6A21 C***TYPE SINGLE PRECISION (ALEAS-S) C***KEYWORDS FNLIB, ALEAS NUMBER, SPECIAL FUNCTIONS, UNIFORM C***AUTHOR Fullerton, W., (LANL) C***DESCRIPTION C C This pseudo-random number generator is portable among a wide C variety of computers. RAND(R) undoubtedly is not as good as many C readily available installation dependent versions, and so this C routine is not recommended for widespread usage. Its redeeming C feature is that the exact same random numbers (to within final round- C off error) can be generated from machine to machine. Thus, programs C that make use of random numbers can be easily transported to and C checked in a new environment. C C The random numbers are generated by the linear congruential C method described, e.g., by Knuth in Seminumerical Methods (p.9), C Addison-Wesley, 1969. Given the I-th number of a pseudo-random C sequence, the I+1 -st number is generated from C X(I+1) = (A*X(I) + C) MOD M, C where here M = 2**22 = 4194304, C = 1731 and several suitable values C of the multiplier A are discussed below. Both the multiplier A and C random number X are represented in double precision as two 11-bit C words. The constants are chosen so that the period is the maximum C possible, 4194304. C C In order that the same numbers be generated from machine to C machine, it is necessary that 23-bit integers be reducible modulo C 2**11 exactly, that 23-bit integers be added exactly, and that 11-bit C integers be multiplied exactly. Furthermore, if the restart option C is used (where R is between 0 and 1), then the product R*2**22 = C R*4194304 must be correct to the nearest integer. C C The first four random numbers should be .0004127026, C .6750836372, .1614754200, and .9086198807. The tenth random number C is .5527787209, and the hundredth is .3600893021 . The thousandth C number should be .2176990509 . C C In order to generate several effectively independent sequences C with the same generator, it is necessary to know the random number C for several widely spaced calls. The I-th random number times 2**22, C where I=K*P/8 and P is the period of the sequence (P = 2**22), is C still of the form L*P/8. In particular we find the I-th random C number multiplied by 2**22 is given by 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 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 C Thus the 4*P/8 = 2097152 random number is 2097152/2**22. C C Several multipliers have been subjected to the spectral test C (see Knuth, p. 82). Four suitable multipliers roughly in order of C goodness according to the spectral test are C 3146757 = 1536*2048 + 1029 = 2**21 + 2**20 + 2**10 + 5 C 2098181 = 1024*2048 + 1029 = 2**21 + 2**10 + 5 C 3146245 = 1536*2048 + 517 = 2**21 + 2**20 + 2**9 + 5 C 2776669 = 1355*2048 + 1629 = 5**9 + 7**7 + 1 C C In the table below LOG10(NU(I)) gives roughly the number of C random decimal digits in the random numbers considered I at a time. C C is the primary measure of goodness. In both cases bigger is better. C C LOG10 NU(I) C(I) C A I=2 I=3 I=4 I=5 I=2 I=3 I=4 I=5 C C 3146757 3.3 2.0 1.6 1.3 3.1 1.3 4.6 2.6 C 2098181 3.3 2.0 1.6 1.2 3.2 1.3 4.6 1.7 C 3146245 3.3 2.2 1.5 1.1 3.2 4.2 1.1 0.4 C 2776669 3.3 2.1 1.6 1.3 2.5 2.0 1.9 2.6 C Best C Possible 3.3 2.3 1.7 1.4 3.6 5.9 9.7 14.9 C C Input Argument -- C R If R=0., the next random number of the sequence is generated. C If R .LT. 0., the last generated number will be returned for C possible use in a restart procedure. C If R .GT. 0., the sequence of random numbers will start with C the seed R mod 1. This seed is also returned as the value of C RAND provided the arithmetic is done exactly. C C Output Value -- C RAND a pseudo-random number between 0. and 1. C C***REFERENCES (NONE) C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 770401 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890531 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE RAND SAVE IA1, IA0, IA1MA0, IC, IX1, IX0 DATA IA1, IA0, IA1MA0 /1536, 1029, 507/ DATA IC /1731/ DATA IX1, IX0 /0, 0/ C***FIRST EXECUTABLE STATEMENT RAND IF (R.LT.0.) GO TO 10 IF (R.GT.0.) GO TO 20 C C A*X = 2**22*IA1*IX1 + 2**11*(IA1*IX1 + (IA1-IA0)*(IX0-IX1) C + IA0*IX0) + IA0*IX0 C IY0 = IA0*IX0 IY1 = IA1*IX1 + IA1MA0*(IX0-IX1) + IY0 IY0 = IY0 + IC IX0 = MOD (IY0, 2048) IY1 = IY1 + (IY0-IX0)/2048 IX1 = MOD (IY1, 2048) C 10 ALEAS = IX1*2048 + IX0 ALEAS = ALEAS / 4194304. RETURN C 20 IX1 = MOD(R,1.)*4194304. + 0.5 IX0 = MOD (IX1, 2048) IX1 = (IX1-IX0)/2048 GO TO 10 C END