*IF DEF,A06_1A GWRICH1A.2
C ******************************COPYRIGHT****************************** GTS2F400.3637
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.3638
C GTS2F400.3639
C Use, duplication or disclosure of this code is subject to the GTS2F400.3640
C restrictions as set forth in the contract. GTS2F400.3641
C GTS2F400.3642
C Meteorological Office GTS2F400.3643
C London Road GTS2F400.3644
C BRACKNELL GTS2F400.3645
C Berkshire UK GTS2F400.3646
C RG12 2SZ GTS2F400.3647
C GTS2F400.3648
C If no contract has been raised with this copy of the code, the use, GTS2F400.3649
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.3650
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.3651
C Modelling at the above address. GTS2F400.3652
C ******************************COPYRIGHT****************************** GTS2F400.3653
C GTS2F400.3654
CLL SUBROUTINE GW_RICH------------------------------------------ GWRICH1A.3
CLL GWRICH1A.4
CLL PURPOSE: TO CALCULATE STRESS PROFILE DUE TO SUBGRID-SCALE GWRICH1A.5
CLL OROGRAPHIC GRAVITY WAVES AND HENCE DRAG ON MEAN FLOW. GWRICH1A.6
CLL THE WAVES PROPOGATE VERTICALLY WITH STRESS INDEPENDENT GWRICH1A.7
CLL HEIGHT UNLESS A CRITICAL LEVEL OR WAVE BREAKING IS GWRICH1A.8
CLL DIAGNOSED. A MINIMUM RICHARDSON NUMBER CRITERION GWRICH1A.9
CLL DETERMINES WHERE THIS OCCURS AND THE WAVE AMPLITUDE GWRICH1A.10
CLL AND STRESS REDUCED SO THAT THE WAVES ARE MAINTAINED GWRICH1A.11
CLL AT MARGINAL STABILITY. GWRICH1A.12
CLL SUITABLE FOR SINGLE COLUMN USE GWRICH1A.13
CLL GWRICH1A.14
CLL SUITABLE FOR ROTATED GRIDS GWRICH1A.15
CLL GWRICH1A.16
CLL ORIGINAL VERSION FOR CRAY Y-MP GWRICH1A.17
CLL WRITTEN BY C. WILSON GWRICH1A.18
CLL FURTHER ALTERATIONS MAY BE REQUIRED FOR AUTOTASKING EFFICIENCY GWRICH1A.19
CLL GWRICH1A.20
CLL MODEL MODIFICATION HISTORY FROM MODEL VERSION 3.0: GWRICH1A.21
CLL VERSION DATE GWRICH1A.22
CLL GWRICH1A.23
CLL 3.3 25/10/93 Removal of DIAG06 directive. New arguments to DR251093.89
CLL dimension diagnostic arrays. D. Robinson. DR251093.90
CLL 3.4 24/09/94 Test correct diagnostic switch for GW Stress ADR4F304.1
CLL v component. D. Robinson. ADR4F304.2
CLL DR251093.91
CLL 4.4 19/09/97 Remove *IF -DEF,CRAY compile options. S.Webster ASW1F404.3
CLL ASW1F404.4
CLL PROGRAMMING STANDARD: UNIFIED MODEL DOCUMENTATION PAPER NO. 4, GWRICH1A.24
CLL VERSION 1, DATED 12/09/89 GWRICH1A.25
CLL GWRICH1A.26
CLL logical components covered: P22 GWRICH1A.27
CLL GWRICH1A.28
CLL SYSTEM TASK: PART OF P22 GWRICH1A.29
CLL GWRICH1A.30
CLL DOCUMENTATION: THE EQUATIONS USED ARE (2.6) TO (2.10) GWRICH1A.31
CLL IN UNIFIED MODEL DOCUMENTATION PAPER NO 22 GWRICH1A.32
CLL C. A. WILSON AND R. SWINBANK GWRICH1A.33
CLL VERSION 1,DATED 15/12/89 GWRICH1A.34
CLLEND------------------------------------------------------------- GWRICH1A.35
GWRICH1A.36
C GWRICH1A.37
C*L ARGUMENTS:--------------------------------------------------- GWRICH1A.38
SUBROUTINE GW_RICH 1GWRICH1A.39
1 (PSTAR,PEXNER,THETA,U,V,S_STRESS,START_L,LEVELS,POINTS, GWRICH1A.40
2 AKH,BKH,DELTA_AK,DELTA_BK,KAY,SIN_A,COS_A, GWRICH1A.41
3 DU_DT,DV_DT, DR251093.92
4 STRESS_UD_LAND,LAND_POINTS_UD,STRESS_UD_ON, DR251093.93
5 STRESS_VD_LAND,LAND_POINTS_VD,STRESS_VD_ON) DR251093.94
GWRICH1A.53
IMPLICIT NONE GWRICH1A.54
GWRICH1A.55
INTEGER GWRICH1A.56
* LEVELS !IN NUMBER OF MODEL LEVELS GWRICH1A.57
*,START_L !IN START LEVEL FOR WAVE-BREAKING TEST GWRICH1A.58
*,POINTS !IN NUMBER OF POINTS GWRICH1A.59
*,LAND_POINTS_UD !IN ) No of land points in diagnostic DR251093.95
*,LAND_POINTS_VD !IN ) arrays for GW stress - u and v. DR251093.96
GWRICH1A.61
REAL GWRICH1A.62
* PSTAR(POINTS) !IN PSTAR FIELD GWRICH1A.63
*,PEXNER(POINTS,LEVELS+1) !IN PEXNER GWRICH1A.64
*,THETA(POINTS,LEVELS) !IN THETA FIELD GWRICH1A.65
*,U(POINTS,LEVELS) !IN U FIELD GWRICH1A.66
*,V(POINTS,LEVELS) !IN V FIELD GWRICH1A.67
*,S_STRESS(POINTS) !IN 'SURFACE' STRESS GWRICH1A.68
C AKH,BKH DEFINE HYBRID VERTICAL COORDINATES P=A+BP*-LAYER EDGES, GWRICH1A.69
C DELTA_AK,DELTA_BK DEFINE PRESSURE DIFFERENCES ACROSS LAYERS GWRICH1A.70
REAL GWRICH1A.71
* AKH(LEVELS+1) !IN VALUE AT LAYER BOUNDARY GWRICH1A.72
*,BKH(LEVELS+1) !IN VALUE AT LAYER BOUMDARY GWRICH1A.73
*,DELTA_AK (LEVELS) !IN DIFFERENCE ACROSS LAYER GWRICH1A.74
*,DELTA_BK (LEVELS) !IN DIFFERENCE ACROSS LAYER GWRICH1A.75
*,KAY !IN stress constant (m-1) GWRICH1A.76
*,SIN_A(POINTS) !IN SIN (STRESS DIRECTION FROM NORTH) GWRICH1A.77
*,COS_A(POINTS) !IN COS (STRESS DIRECTION FROM NORTH) GWRICH1A.78
*,DU_DT(POINTS,LEVELS) !OUT U-ACCELERATION GWRICH1A.79
*,DV_DT(POINTS,LEVELS) !OUT V-ACCELERATION GWRICH1A.80
*,STRESS_UD_LAND(LAND_POINTS_UD,LEVELS+1) !U STRESS DIAGNOSTIC DR251093.97
*,STRESS_VD_LAND(LAND_POINTS_VD,LEVELS+1) !V STRESS DIAGNOSTIC DR251093.98
GWRICH1A.81
LOGICAL DR251093.99
+ STRESS_UD_ON ! ) Diagnostic switches for GW stress - DR251093.100
+ ,STRESS_VD_ON ! ) u and v (Item Nos 201 and 202) DR251093.101
GWRICH1A.89
C*--------------------------------------------------------------------- GWRICH1A.90
GWRICH1A.91
C*L WORKSPACE USAGE:------------------------------------------------- GWRICH1A.92
C DEFINE LOCAL WORKSPACE ARRAYS: GWRICH1A.93
C 10 REAL ARRAYS OF FULL FIELD LENGTH REQUIRED GWRICH1A.94
C GWRICH1A.95
GWRICH1A.101
REAL GWRICH1A.102
* DZ(POINTS,3) ! HEIGHT DIFFERENCES IN EACH HALF LAYER GWRICH1A.103
*,T(POINTS,2) ! TEMPERATURES (LEVELS) GWRICH1A.104
*,SPEED(POINTS,2) ! WIND SPEEDS (LEVELS) GWRICH1A.105
*,STRESS(POINTS,2) ! STRESSES (LAYER BOUNDARIES) GWRICH1A.106
*,DRAG(POINTS) ! DRAG EXERTED ON LAYER(IN DIRECTION OF GWRICH1A.107
* ! SURFACE STRESS GWRICH1A.108
GWRICH1A.110
C*--------------------------------------------------------------------- GWRICH1A.111
C GWRICH1A.112
C*L EXTERNAL SUBROUTINES CALLED--------------------------------------- GWRICH1A.113
C NONE GWRICH1A.114
C*------------------------------------------------------------------ GWRICH1A.115
CL MAXIMUM VECTOR LENGTH ASSUMED = POINTS GWRICH1A.116
CL--------------------------------------------------------------------- GWRICH1A.117
C---------------------------------------------------------------------- GWRICH1A.118
C GWRICH1A.119
C DEFINE LOCAL VARIABLES GWRICH1A.120
C LOCAL VARIABLES: GWRICH1A.121
C GWRICH1A.122
REAL GWRICH1A.123
* RHO ! DENSITY AT LAYER BOUNDARY GWRICH1A.124
*,TB ! TEMPERATURE AT LAYER BOUNDARY GWRICH1A.125
*,SPEEDB ! WIND SPEEDS AT LAYER BOUNDARY GWRICH1A.126
*,DZB ! HEIGHT DIFFERENCE ACROSS LAYER BOUNDARY GWRICH1A.127
*,N ! BRUNT_VAISALA FREQUENCY GWRICH1A.128
*,N_SQ ! SQUARE OF BRUNT_VAISALA FREQUENCY GWRICH1A.129
*,AMPL ! WAVE-AMPLITUDE GWRICH1A.130
*,DV_DZ ! MAGNITUDE OF WIND SHEAR GWRICH1A.131
*,RI ! MINIMUM RICHARDSON NUMBER GWRICH1A.132
*,EPSILON ! MAX WAVE-AMPLITUDE*N**2/WIND SPEED GWRICH1A.133
*,DELTA_P ! DIFFERENCE IN PRESSURE ACROSS LAYER GWRICH1A.134
REAL GWRICH1A.135
* DELTA_AK_SUM ! DELTA_AK SUMMED OVER LOWEST LAYERS UP TO START_L GWRICH1A.136
*,DELTA_BK_SUM ! DELTA_BK SUMMED OVER LOWEST LAYERS UP TO START_L GWRICH1A.137
INTEGER I,K ! LOOP COUNTER IN ROUTINE GWRICH1A.138
INTEGER KK,KL,KU,KT ! LEVEL COUNTERS IN ROUTINE GWRICH1A.139
C GWRICH1A.140
C INCLUDE PHYSICAL CONSTANTS GWRICH1A.141
*CALL C_G 
GWRICH1A.142
*CALL C_R_CP GWRICH1A.143
GWRICH1A.144
C LOCAL CONSTANTS GWRICH1A.145
*CALL C_GWAVE GWRICH1A.146
GWRICH1A.147
REAL CPBYG GWRICH1A.148
PARAMETER(CPBYG=CP/G) GWRICH1A.149
GWRICH1A.150
REAL GWRICH1A.151
& PU,PL,P_EXNER_CENTRE GWRICH1A.152
*CALL P_EXNERC GWRICH1A.153
GWRICH1A.154
C------------------------------------------------------------------- GWRICH1A.155
CL INTERNAL STRUCTURE INCLUDING SUBROUTINE CALLS: GWRICH1A.156
CL 1. START LEVEL PRELIMINARIES GWRICH1A.157
C------------------------------------------ GWRICH1A.158
GWRICH1A.159
CFPP$ NOCONCUR L GWRICH1A.160
C TREAT LAYERS BELOW AND INCLUDING START_L AS ONE LAYER GWRICH1A.161
DELTA_AK_SUM = 0.0 GWRICH1A.162
DELTA_BK_SUM = 0.0 GWRICH1A.163
DO K=1,START_L GWRICH1A.164
DELTA_AK_SUM = DELTA_AK_SUM + DELTA_AK(K) GWRICH1A.165
DELTA_BK_SUM = DELTA_BK_SUM + DELTA_BK(K) GWRICH1A.166
END DO GWRICH1A.167
CFPP$ CONCUR GWRICH1A.168
GWRICH1A.169
KL=1 GWRICH1A.170
KU=2 GWRICH1A.171
KT=3 GWRICH1A.172
GWRICH1A.173
DO I=1,POINTS GWRICH1A.174
GWRICH1A.175
SPEED(I,KL) = U(I,START_L)*SIN_A(I) + V(I,START_L)*COS_A(I) GWRICH1A.176
STRESS(I,KL) = S_STRESS(I) GWRICH1A.177
GWRICH1A.178
PU=PSTAR(I)*BKH(START_L+1) + AKH(START_L+1) GWRICH1A.179
PL=PSTAR(I)*BKH(START_L) + AKH(START_L) GWRICH1A.180
P_EXNER_CENTRE= GWRICH1A.181
& P_EXNER_C( PEXNER(I,START_L+1),PEXNER(I,START_L),PU,PL,KAPPA) GWRICH1A.182
GWRICH1A.183
DZ(I,KL) = (P_EXNER_CENTRE - PEXNER(I,START_L+1)) GWRICH1A.184
* *THETA(I,START_L)*CPBYG GWRICH1A.185
T(I,KL) = P_EXNER_CENTRE*THETA(I,START_L) GWRICH1A.186
GWRICH1A.193
END DO GWRICH1A.194
DR251093.102
IF (STRESS_UD_ON) THEN DR251093.103
DO I=1,POINTS DR251093.104
STRESS_UD_LAND(I,START_L) = STRESS(I,KL)*SIN_A(I) DR251093.105
ENDDO DR251093.106
ENDIF DR251093.107
IF (STRESS_VD_ON) THEN ADR4F304.3
DO I=1,POINTS DR251093.109
STRESS_VD_LAND(I,START_L) = STRESS(I,KL)*COS_A(I) DR251093.110
ENDDO DR251093.111
ENDIF DR251093.112
GWRICH1A.195
C------------------------------------------------------------------ GWRICH1A.196
CL 2 LOOP LEVELS GWRICH1A.197
C------------------------------------------------------------------ GWRICH1A.198
GWRICH1A.199
DO K=START_L+1,LEVELS GWRICH1A.200
GWRICH1A.201
GWRICH1A.202
DO I=1,POINTS GWRICH1A.203
STRESS(I,KU) = STRESS(I,KL) GWRICH1A.204
IF(STRESS(I,KL) .GT. 0.0) THEN GWRICH1A.205
SPEED(I,KU) = U(I,K)*SIN_A(I) + V(I,K)*COS_A(I) GWRICH1A.206
GWRICH1A.207
PU=PSTAR(I)*BKH(K+1) + AKH(K+1) GWRICH1A.208
PL=PSTAR(I)*BKH(K) + AKH(K) GWRICH1A.209
P_EXNER_CENTRE= GWRICH1A.210
& P_EXNER_C( PEXNER(I,K+1),PEXNER(I,K),PU,PL,KAPPA) GWRICH1A.211
GWRICH1A.212
C lower half height of upper layer GWRICH1A.213
DZ(I,KU) = (PEXNER(I,K) - P_EXNER_CENTRE)*THETA(I,K) GWRICH1A.214
* *CPBYG GWRICH1A.215
C upper half height of upper layer GWRICH1A.216
DZ(I,KT) = (P_EXNER_CENTRE - PEXNER(I,K+1))*THETA(I,K) GWRICH1A.217
* *CPBYG GWRICH1A.218
C model level height difference GWRICH1A.219
DZB = DZ(I,KU) + DZ(I,KL) GWRICH1A.220
SPEEDB = (DZ(I,KU)*SPEED(I,KL)+DZ(I,KL)*SPEED(I,KU))/ GWRICH1A.221
* DZB GWRICH1A.222
GWRICH1A.223
C------------------------------------------------------------------ GWRICH1A.224
CL 2.1 TEST FOR CRITICAL LEVEL V < or = 0 GWRICH1A.225
C------------------------------------------------------------------ GWRICH1A.226
GWRICH1A.227
IF( SPEEDB .LE. 0.0) THEN GWRICH1A.228
STRESS(I,KU) = 0.0 GWRICH1A.229
ELSE GWRICH1A.230
T(I,KU) = P_EXNER_CENTRE*THETA(I,K) GWRICH1A.231
TB = (DZ(I,KU)*T(I,KL) + DZ(I,KL)*T(I,KU))/DZB GWRICH1A.232
RHO = ( AKH(K) + BKH(K)*PSTAR(I) )/(R*TB) GWRICH1A.233
GWRICH1A.234
C------------------------------------------------------------------ GWRICH1A.235
CL 2.2 CALCULATE BRUNT-VAISALA FREQUENCY GWRICH1A.236
C------------------------------------------------------------------ GWRICH1A.237
GWRICH1A.238
N_SQ = G*( THETA(I,K) - THETA(I,K-1) )*PEXNER(I,K)/ GWRICH1A.239
* ( TB*DZB ) GWRICH1A.240
GWRICH1A.241
IF( N_SQ .LE.0.0 ) THEN GWRICH1A.242
CL SET STRESS TO ZERO IF UNSTABLE GWRICH1A.243
N_SQ = 0.0 GWRICH1A.244
STRESS(I,KU) = 0.0 GWRICH1A.245
ELSE GWRICH1A.246
N = SQRT( N_SQ ) GWRICH1A.247
GWRICH1A.248
C------------------------------------------------------------------ GWRICH1A.249
CL 2.2 CALCULATE MINIMUM RICHARDSON NO. DUE TO GRAVITY GWRICH1A.250
CL WAVES EQNS 2.6 AND 2.8 GWRICH1A.251
C------------------------------------------------------------------ GWRICH1A.252
GWRICH1A.253
AMPL = GWRICH1A.254
* SQRT(STRESS(I,KU)/( KAY*RHO*SPEEDB*N ) ) GWRICH1A.255
DV_DZ = SQRT( (U(I,K)-U(I,K-1))*(U(I,K)-U(I,K-1)) + GWRICH1A.256
* (V(I,K)-V(I,K-1))*(V(I,K)-V(I,K-1)) )/DZB GWRICH1A.257
RI = N_SQ*(1. - N*AMPL/SPEEDB)/ GWRICH1A.258
* ((DV_DZ + N_SQ*AMPL/SPEEDB)*(DV_DZ + N_SQ*AMPL/SPEEDB)) GWRICH1A.259
GWRICH1A.260
C------------------------------------------------------------------ GWRICH1A.261
CL 2.3 TEST FOR WAVE-BREAKING (RI < RIC ) GWRICH1A.262
CL AND MODIFY STRESS AT UPPER LAYER BOUNDARY GWRICH1A.263
CL EQNS 2.9 , 2.10 AND 2.6 GWRICH1A.264
C------------------------------------------------------------------ GWRICH1A.265
GWRICH1A.266
IF( RI.LT.RIC ) THEN GWRICH1A.267
EPSILON = GWRICH1A.268
* ( SQRT(4.*RIC*DV_DZ*N + N_SQ*(1.+4.*RIC)) GWRICH1A.269
* -(2.*RIC*DV_DZ + N ) )/(2.*RIC) GWRICH1A.270
AMPL = EPSILON*SPEEDB/ N_SQ GWRICH1A.271
IF( AMPL.LT.0.0 ) THEN GWRICH1A.272
STRESS(I,KU) = 0.0 GWRICH1A.273
ELSE GWRICH1A.274
STRESS(I,KU) = KAY*RHO*N*SPEEDB*AMPL*AMPL GWRICH1A.275
GWRICH1A.276
END IF ! AMPL < 0 ELSE AMPL > 0 GWRICH1A.277
GWRICH1A.278
END IF ! RI < RIC GWRICH1A.279
GWRICH1A.280
END IF ! N_SQ < 0 ELSE N_SQ > 0 GWRICH1A.281
GWRICH1A.282
END IF ! SPEED < 0 ELSE SPEED >0 GWRICH1A.283
GWRICH1A.284
END IF ! STRESS >0 GWRICH1A.285
GWRICH1A.286
IF( K .EQ. START_L+1 ) THEN GWRICH1A.287
DELTA_P = DELTA_AK_SUM+DELTA_BK_SUM*PSTAR(I) GWRICH1A.288
ELSE GWRICH1A.289
DELTA_P = DELTA_AK(K-1)+DELTA_BK(K-1)*PSTAR(I) GWRICH1A.290
END IF GWRICH1A.291
GWRICH1A.292
C------------------------------------------------------------------ GWRICH1A.293
CL 2.4 CALCULATE DRAG FROM VERTICAL STRESS CONVERGENCE GWRICH1A.294
CL AND ACCELERATIONS FOR WIND COMPONENTS GWRICH1A.295
CL EQN 2.1 GWRICH1A.296
C------------------------------------------------------------------ GWRICH1A.297
GWRICH1A.298
DRAG(I) = G*(STRESS(I,KU) - STRESS(I,KL) )/DELTA_P GWRICH1A.299
DU_DT(I,K-1) = -DRAG(I)*SIN_A(I) GWRICH1A.300
DV_DT(I,K-1) = -DRAG(I)*COS_A(I) GWRICH1A.301
GWRICH1A.302
END DO GWRICH1A.309
GWRICH1A.310
IF (STRESS_UD_ON) THEN DR251093.113
DO I=1,POINTS DR251093.114
STRESS_UD_LAND(I,K) = STRESS(I,KU)*SIN_A(I) DR251093.115
ENDDO DR251093.116
ENDIF DR251093.117
IF (STRESS_VD_ON) THEN DR251093.118
DO I=1,POINTS DR251093.119
STRESS_VD_LAND(I,K) = STRESS(I,KU)*COS_A(I) DR251093.120
ENDDO DR251093.121
ENDIF DR251093.122
GWRICH1A.311
IF( K .EQ. START_L+1 .AND. START_L .GT.1 ) THEN GWRICH1A.312
C SET ACCELERATION SAME IN ALL LAYERS UP TO START_L GWRICH1A.313
DO KK=1,START_L-1 GWRICH1A.314
DO I=1,POINTS GWRICH1A.315
DU_DT(I,KK) = DU_DT(I,START_L) GWRICH1A.316
DV_DT(I,KK) = DV_DT(I,START_L) GWRICH1A.317
END DO GWRICH1A.318
END DO GWRICH1A.319
END IF GWRICH1A.320
GWRICH1A.321
C Swap storage for lower and upper layers GWRICH1A.322
KK=KL GWRICH1A.323
KL=KU GWRICH1A.324
KU=KK GWRICH1A.325
GWRICH1A.326
C Replace top half height of lower layer ready for next pass GWRICH1A.327
DO I=1,POINTS GWRICH1A.328
DZ(I,KL)=DZ(I,KT) GWRICH1A.329
ENDDO GWRICH1A.330
GWRICH1A.331
END DO GWRICH1A.332
CL END LOOP LEVELS GWRICH1A.333
GWRICH1A.334
C------------------------------------------------------------------ GWRICH1A.335
CL 3 TOP OF MODEL. SET ACCELERATION SAME AS PENULTIMATE LAYER GWRICH1A.336
CL WITH PROVISO THAT STRESS >= 0 GWRICH1A.337
C------------------------------------------------------------------ GWRICH1A.338
GWRICH1A.339
DO I=1,POINTS GWRICH1A.340
DELTA_P = DELTA_AK(LEVELS) + DELTA_BK(LEVELS)*PSTAR(I) GWRICH1A.341
STRESS(I,KU) = STRESS(I,KL) + DRAG(I)*DELTA_P/G GWRICH1A.342
IF( STRESS(I,KU) .LT. 0.0) STRESS(I,KU) = 0.0 GWRICH1A.343
DRAG(I) = G*(STRESS(I,KU) - STRESS(I,KL) )/DELTA_P GWRICH1A.344
DU_DT(I,LEVELS) = -DRAG(I)*SIN_A(I) GWRICH1A.345
DV_DT(I,LEVELS) = -DRAG(I)*COS_A(I) GWRICH1A.346
END DO DR251093.123
GWRICH1A.347
IF (STRESS_UD_ON) THEN DR251093.124
DO I=1,POINTS DR251093.125
STRESS_UD_LAND(I,LEVELS+1) = STRESS(I,KU)*SIN_A(I) DR251093.126
END DO DR251093.127
ENDIF DR251093.128
IF (STRESS_VD_ON) THEN DR251093.129
DO I=1,POINTS DR251093.130
STRESS_VD_LAND(I,LEVELS+1) = STRESS(I,KU)*COS_A(I) DR251093.131
END DO DR251093.132
ENDIF DR251093.133
GWRICH1A.355
RETURN GWRICH1A.356
END GWRICH1A.357
GWRICH1A.358
*ENDIF GWRICH1A.359