*IF DEF,A03_7A SFOROG7A.2
C *****************************COPYRIGHT****************************** SFOROG7A.3
C (c) CROWN COPYRIGHT 1997, METEOROLOGICAL OFFICE, All Rights Reserved. SFOROG7A.4
C SFOROG7A.5
C Use, duplication or disclosure of this code is subject to the SFOROG7A.6
C restrictions as set forth in the contract. SFOROG7A.7
C SFOROG7A.8
C Meteorological Office SFOROG7A.9
C London Road SFOROG7A.10
C BRACKNELL SFOROG7A.11
C Berkshire UK SFOROG7A.12
C RG12 2SZ SFOROG7A.13
C SFOROG7A.14
C If no contract has been raised with this copy of the code, the use, SFOROG7A.15
C duplication or disclosure of it is strictly prohibited. Permission SFOROG7A.16
C to do so must first be obtained in writing from the Head of Numerical SFOROG7A.17
C Modelling at the above address. SFOROG7A.18
C ******************************COPYRIGHT****************************** SFOROG7A.19
! SFOROG7A.20
!!! SUBROUTINE SF_OROG------------------------------------------------ SFOROG7A.21
!!! SFOROG7A.22
!!! Purpose: Calculate roughness lengths, blending height and wind SFOROG7A.23
!!! profile factor SFOROG7A.24
!!! SFOROG7A.25
!!! SJ, RE <- programmer of some or all of previous code changes SFOROG7A.26
C Modification History: AJC1F405.63
C Version Date Change AJC1F405.64
C 4.5 Jul. 98 Kill the IBM specific lines (JCThil) AJC1F405.65
!!! SFOROG7A.27
!!!-------------------------------------------------------------------- SFOROG7A.28
SUBROUTINE SF_OROG( 2,2SFOROG7A.29
& P_FIELD,LAND_FIELD,TILE_PTS,LAND_INDEX,TILE_INDEX, SFOROG7A.30
& L_Z0_OROG,LTIMER, SFOROG7A.31
& HO2R2_OROG,RIB,SIL_OROG,Z0M,Z1, SFOROG7A.32
& WIND_PROFILE_FACTOR,Z0M_EFF SFOROG7A.33
& ) SFOROG7A.34
SFOROG7A.35
IMPLICIT NONE SFOROG7A.36
SFOROG7A.37
INTEGER SFOROG7A.38
& P_FIELD ! IN Number of points on P-grid. SFOROG7A.39
&,LAND_FIELD ! IN Number of land points. SFOROG7A.40
&,TILE_PTS ! IN Number of tile points. SFOROG7A.41
&,LAND_INDEX(P_FIELD) ! IN Index of land points. SFOROG7A.42
&,TILE_INDEX(LAND_FIELD)! IN Index of tile points. SFOROG7A.43
SFOROG7A.44
LOGICAL SFOROG7A.45
& LTIMER ! IN .TRUE. for timer diagnostics SFOROG7A.46
&,L_Z0_OROG ! IN .TRUE. to use orographic roughness. SFOROG7A.47
SFOROG7A.48
REAL SFOROG7A.49
& HO2R2_OROG(LAND_FIELD)!IN Peak to trough height of unresolved SFOROG7A.50
! ! orography divided by 2SQRT(2) (m). SFOROG7A.51
&,RIB(LAND_FIELD) ! IN Bulk Richardson number for lowest layer SFOROG7A.52
&,SIL_OROG(LAND_FIELD) ! IN Silhouette area of unresolved orography SFOROG7A.53
! ! per unit horizontal area SFOROG7A.54
&,Z0M(LAND_FIELD) ! IN Roughness length for momentum (m). SFOROG7A.55
&,Z1(P_FIELD) ! IN Height of lowest atmospheric level (m). SFOROG7A.56
SFOROG7A.57
! Output variables. SFOROG7A.58
SFOROG7A.59
REAL SFOROG7A.60
& WIND_PROFILE_FACTOR(LAND_FIELD) SFOROG7A.61
! ! OUT For transforming effective surface SFOROG7A.62
! ! transfer coefficients to those SFOROG7A.63
! ! excluding form drag. SFOROG7A.64
&,Z0M_EFF(LAND_FIELD) ! OUT Effective roughness length for SFOROG7A.65
! momentum (m) SFOROG7A.66
SFOROG7A.67
! Work Variables SFOROG7A.68
SFOROG7A.69
INTEGER SFOROG7A.70
& I ! Horizontal field index SFOROG7A.71
&,J ! Tile field index SFOROG7A.72
&,L ! Land field index SFOROG7A.73
SFOROG7A.74
REAL SFOROG7A.75
& H_BLEND_OROG ! Blending height SFOROG7A.76
&,RIB_FN ! Interpolation coefficient for 0 < RIB < RI_CRIT SFOROG7A.77
&,ZETA1 ! Work space SFOROG7A.78
&,ZETA2 ! More work space SFOROG7A.79
&,ZETA3 ! Even more work space SFOROG7A.80
SFOROG7A.81
! Common parameters SFOROG7A.82
SFOROG7A.83
*CALL C_MDI
SFOROG7A.84
*CALL C_VKMAN
SFOROG7A.85
*CALL C_SURF
SFOROG7A.86
SFOROG7A.87
! Local parameters SFOROG7A.88
SFOROG7A.89
REAL H_BLEND_MIN,H_BLEND_MAX SFOROG7A.90
PARAMETER ( SFOROG7A.91
& H_BLEND_MIN=0.0 ! Minimun value of blending height SFOROG7A.92
&,H_BLEND_MAX=1000.0 ! Maximum value of blending height SFOROG7A.93
& ) SFOROG7A.94
SFOROG7A.95
EXTERNAL TIMER SFOROG7A.96
SFOROG7A.97
IF (LTIMER) THEN SFOROG7A.98
CALL TIMER
('SF_OROG ',3) SFOROG7A.99
ENDIF SFOROG7A.100
SFOROG7A.101
! Set blending height, effective roughness length and SFOROG7A.102
! wind profile factor at land points. SFOROG7A.103
DO J=1,TILE_PTS SFOROG7A.104
L = TILE_INDEX(J) SFOROG7A.105
I = LAND_INDEX(L) SFOROG7A.106
SFOROG7A.107
WIND_PROFILE_FACTOR(L) = 1.0 SFOROG7A.108
Z0M_EFF(L) = Z0M(L) SFOROG7A.109
SFOROG7A.110
IF (L_Z0_OROG) THEN SFOROG7A.111
SFOROG7A.112
ZETA1 = 0.5 * OROG_DRAG_PARAM * SIL_OROG(L) SFOROG7A.113
ZETA2 = LOG ( 1.0 + Z1(I)/Z0M(L) ) SFOROG7A.114
ZETA3 = 1.0 / SQRT ( ZETA1/(VKMAN*VKMAN) + 1.0/(ZETA2*ZETA2) ) SFOROG7A.115
ZETA2 = 1.0 / EXP(ZETA3) SFOROG7A.116
H_BLEND_OROG = MAX ( Z1(I) / (1.0 - ZETA2) , SFOROG7A.117
& HO2R2_OROG(L) * SQRT(2.0) ) SFOROG7A.118
H_BLEND_OROG = MIN ( H_BLEND_MAX, H_BLEND_OROG ) SFOROG7A.119
SFOROG7A.120
! Apply simple stability correction to form drag if RIB is stable SFOROG7A.121
SFOROG7A.122
IF (SIL_OROG(L) .EQ. RMDI) THEN SFOROG7A.123
ZETA1 = 0.0 SFOROG7A.124
ELSE SFOROG7A.125
RIB_FN = ( 1.0 - RIB(L) / RI_CRIT ) SFOROG7A.126
IF (RIB_FN.GT.1.0) RIB_FN = 1.0 SFOROG7A.127
IF (RIB_FN.LT.0.0) RIB_FN = 0.0 SFOROG7A.128
ZETA1 = 0.5 * OROG_DRAG_PARAM * SIL_OROG(L) * RIB_FN SFOROG7A.129
ENDIF SFOROG7A.130
SFOROG7A.131
Z0M_EFF(L) = H_BLEND_OROG / EXP ( VKMAN / SQRT ( ZETA1 + SFOROG7A.132
& (VKMAN / LOG ( H_BLEND_OROG / Z0M(L) ) ) * SFOROG7A.133
& (VKMAN / LOG ( H_BLEND_OROG / Z0M(L) ) ) ) ) SFOROG7A.134
IF (RIB(L).GT.RI_CRIT) Z0M_EFF(L)=Z0M(L) SFOROG7A.135
SFOROG7A.136
WIND_PROFILE_FACTOR(L) = LOG( H_BLEND_OROG / Z0M_EFF(L) ) / SFOROG7A.137
& LOG( H_BLEND_OROG / Z0M(L) ) SFOROG7A.138
SFOROG7A.139
ENDIF SFOROG7A.140
SFOROG7A.141
ENDDO SFOROG7A.142
SFOROG7A.143
IF (LTIMER) THEN SFOROG7A.144
CALL TIMER
('SF_OROG ',4) SFOROG7A.145
ENDIF SFOROG7A.146
SFOROG7A.147
RETURN SFOROG7A.148
END SFOROG7A.149
SFOROG7A.150
!!! SUBROUTINE SF_OROG_GB -------------------------------------------- SFOROG7A.151
!!! SFOROG7A.152
!!! Purpose: Calculate effective roughness length and blending height SFOROG7A.153
!!! SFOROG7A.154
!!! SJ, RE <- programmer of some or all of previous code changes SFOROG7A.155
!!! SFOROG7A.156
!!!-------------------------------------------------------------------- SFOROG7A.157
SUBROUTINE SF_OROG_GB( 1,2SFOROG7A.158
& P_FIELD,P1,P_POINTS,LAND_FIELD,LAND1,LAND_PTS,LAND_INDEX, SFOROG7A.159
& LAND_MASK,L_Z0_OROG,HO2R2_OROG,RIB,SIL_OROG,Z0M,Z1, SFOROG7A.160
& H_BLEND_OROG,Z0M_EFF,LTIMER SFOROG7A.161
& ) SFOROG7A.162
SFOROG7A.163
IMPLICIT NONE SFOROG7A.164
SFOROG7A.165
INTEGER SFOROG7A.166
& P_FIELD ! IN Number of points on P-grid. SFOROG7A.167
&,P1 ! IN First P-point to be processed. SFOROG7A.168
&,P_POINTS ! IN Number of P-grid points to be processed SFOROG7A.169
&,LAND_FIELD ! IN Number of land points. SFOROG7A.170
&,LAND1 ! IN First land point to be processed. SFOROG7A.171
&,LAND_PTS ! IN Number of land points to be processed. SFOROG7A.172
&,LAND_INDEX(P_FIELD) ! IN Index of land points. SFOROG7A.173
SFOROG7A.174
LOGICAL SFOROG7A.175
& LAND_MASK(P_FIELD) ! IN .TRUE. for land; .FALSE. elsewhere. SFOROG7A.176
&,LTIMER ! IN .TRUE. for timer diagnostics SFOROG7A.177
&,L_Z0_OROG ! IN .TRUE. to use orographic roughness. SFOROG7A.178
SFOROG7A.179
SFOROG7A.180
REAL SFOROG7A.181
& HO2R2_OROG(LAND_FIELD)!IN Peak to trough height of unresolved SFOROG7A.182
! ! orography divided by 2SQRT(2) (m). SFOROG7A.183
&,RIB(P_FIELD) ! IN GBM Bulk Richardson number for lowest SFOROG7A.184
! ! layer SFOROG7A.185
&,SIL_OROG(LAND_FIELD) ! IN Silhouette area of unresolved orography SFOROG7A.186
! ! per unit horizontal area SFOROG7A.187
&,Z0M(LAND_FIELD) ! IN GBM Roughness length for momentum (m). SFOROG7A.188
&,Z1(P_FIELD) ! IN Height of lowest atmospheric level (m). SFOROG7A.189
SFOROG7A.190
! Output variables. SFOROG7A.191
SFOROG7A.192
REAL SFOROG7A.193
& H_BLEND_OROG(P_FIELD)! OUT Blending height SFOROG7A.194
&,Z0M_EFF(P_FIELD) ! OUT Effective roughness length for SFOROG7A.195
! momentum (m) SFOROG7A.196
SFOROG7A.197
! Work Variables SFOROG7A.198
SFOROG7A.199
INTEGER SFOROG7A.200
& I ! Horizontal field index SFOROG7A.201
&,L ! Land field index SFOROG7A.202
SFOROG7A.203
REAL SFOROG7A.204
& RIB_FN ! Interpolation coefficient for 0 < RIB < RI_CRIT SFOROG7A.205
&,ZETA1 ! Work space SFOROG7A.206
&,ZETA2 ! More work space SFOROG7A.207
&,ZETA3 ! Even more work space SFOROG7A.208
SFOROG7A.209
! Common parameters SFOROG7A.210
SFOROG7A.211
*CALL C_MDI
SFOROG7A.212
*CALL C_VKMAN
SFOROG7A.213
*CALL C_SURF
SFOROG7A.214
SFOROG7A.215
! Local parameters SFOROG7A.216
SFOROG7A.217
REAL H_BLEND_MIN,H_BLEND_MAX SFOROG7A.218
PARAMETER ( SFOROG7A.219
& H_BLEND_MIN=0.0 ! Minimun value of blending height SFOROG7A.220
&,H_BLEND_MAX=1000.0 ! Maximum value of blending height SFOROG7A.221
& ) SFOROG7A.222
SFOROG7A.223
EXTERNAL TIMER SFOROG7A.224
SFOROG7A.225
IF (LTIMER) THEN SFOROG7A.226
CALL TIMER
('SF_OROG ',3) SFOROG7A.227
ENDIF SFOROG7A.228
SFOROG7A.229
! Set blending height, effective roughness length and SFOROG7A.230
! wind profile factor at land points. SFOROG7A.231
SFOROG7A.232
DO L = LAND1,LAND1+LAND_PTS-1 SFOROG7A.238
I = LAND_INDEX(L) SFOROG7A.239
SFOROG7A.241
H_BLEND_OROG(I) = H_BLEND_MIN SFOROG7A.242
Z0M_EFF(I) = Z0M(L) SFOROG7A.243
SFOROG7A.244
IF (L_Z0_OROG) THEN SFOROG7A.245
SFOROG7A.246
ZETA1 = 0.5 * OROG_DRAG_PARAM * SIL_OROG(L) SFOROG7A.247
ZETA2 = LOG ( 1.0 + Z1(I)/Z0M(L) ) SFOROG7A.248
ZETA3 = 1.0 / SQRT ( ZETA1/(VKMAN*VKMAN) + 1.0/(ZETA2*ZETA2) ) SFOROG7A.249
ZETA2 = 1.0 / EXP(ZETA3) SFOROG7A.250
H_BLEND_OROG(I) = MAX ( Z1(I) / (1.0 - ZETA2) , ABX1F405.927
& HO2R2_OROG(L) * SQRT(2.0) ) SFOROG7A.252
H_BLEND_OROG(I) = MIN ( H_BLEND_MAX, H_BLEND_OROG(I) ) ABX1F405.928
SFOROG7A.254
! Apply simple stability correction to form drag if RIB is stable SFOROG7A.255
SFOROG7A.256
IF (SIL_OROG(L) .EQ. RMDI) THEN SFOROG7A.257
ZETA1 = 0.0 SFOROG7A.258
ELSE SFOROG7A.259
RIB_FN = ( 1.0 - RIB(I) / RI_CRIT ) SFOROG7A.260
IF (RIB_FN.GT.1.0) RIB_FN = 1.0 SFOROG7A.261
IF (RIB_FN.LT.0.0) RIB_FN = 0.0 SFOROG7A.262
ZETA1 = 0.5 * OROG_DRAG_PARAM * SIL_OROG(L) * RIB_FN SFOROG7A.263
ENDIF SFOROG7A.264
SFOROG7A.265
Z0M_EFF(I) = H_BLEND_OROG(I) / EXP ( VKMAN / SQRT ( ZETA1 + ABX1F405.929
& (VKMAN / LOG ( H_BLEND_OROG(I) / Z0M(L) ) ) * ABX1F405.930
& (VKMAN / LOG ( H_BLEND_OROG(I) / Z0M(L) ) ) ) ) ABX1F405.931
IF (RIB(I).GT.RI_CRIT) Z0M_EFF(I) = Z0M(L) SFOROG7A.269
SFOROG7A.270
ENDIF SFOROG7A.271
SFOROG7A.272
ENDDO SFOROG7A.276
SFOROG7A.277
IF (LTIMER) THEN SFOROG7A.278
CALL TIMER
('SF_OROG ',4) SFOROG7A.279
ENDIF SFOROG7A.280
SFOROG7A.281
RETURN SFOROG7A.282
END SFOROG7A.283
*ENDIF SFOROG7A.284