*IF DEF,A03_6A PHIMH6A.2
C ******************************COPYRIGHT****************************** PHIMH6A.3
C (c) CROWN COPYRIGHT 1997, METEOROLOGICAL OFFICE, All Rights Reserved. PHIMH6A.4
C PHIMH6A.5
C Use, duplication or disclosure of this code is subject to the PHIMH6A.6
C restrictions as set forth in the contract. PHIMH6A.7
C PHIMH6A.8
C Meteorological Office PHIMH6A.9
C London Road PHIMH6A.10
C BRACKNELL PHIMH6A.11
C Berkshire UK PHIMH6A.12
C RG12 2SZ PHIMH6A.13
C PHIMH6A.14
C If no contract has been raised with this copy of the code, the use, PHIMH6A.15
C duplication or disclosure of it is strictly prohibited. Permission PHIMH6A.16
C to do so must first be obtained in writing from the Head of Numerical PHIMH6A.17
C Modelling at the above address. PHIMH6A.18
C ******************************COPYRIGHT****************************** PHIMH6A.19
C PHIMH6A.20
!!! SUBROUTINE PHI_M_H ---------------------------------------------- PHIMH6A.21
!!! PHIMH6A.22
!!! Purpose: Calculate the integrated froms of the Monin-Obukhov PHIMH6A.23
!!! stability functions for surface exchanges. PHIMH6A.24
!!! PHIMH6A.25
!!! PHIMH6A.26
!!! Model Modification history: PHIMH6A.27
!!! version Date PHIMH6A.28
!!! PHIMH6A.29
!!! 4.4 12/05/97 *DECK created by R.N.B.Smith PHIMH6A.30
!!! PHIMH6A.31
!!! Programming standard: PHIMH6A.32
!!! PHIMH6A.33
!!! System component covered: Part of P243. PHIMH6A.34
!!! PHIMH6A.35
!!! Documentation: UM Documentation Paper No 24. PHIMH6A.36
!!! PHIMH6A.37
!*L Arguments:--------------------------------------------------------- PHIMH6A.38
SUBROUTINE PHI_M_H( 4,2PHIMH6A.39
& P_POINTS,P_FIELD,P1,L_LAND,LAND_MASK, PHIMH6A.40
& RECIP_L_MO,Z_UV,Z_TQ,Z0M,Z0H,PHI_M,PHI_H,LTIMER PHIMH6A.41
&) PHIMH6A.42
IMPLICIT NONE PHIMH6A.43
PHIMH6A.44
INTEGER PHIMH6A.45
& P_POINTS ! IN Number of gridpoints treated. PHIMH6A.46
&,P_FIELD ! IN Size of field on p-grid. PHIMH6A.47
&,P1 ! IN First p-point to be treated. PHIMH6A.48
PHIMH6A.49
LOGICAL PHIMH6A.50
& LTIMER PHIMH6A.51
&,L_LAND ! IN If true, treat land points only. PHIMH6A.52
&,LAND_MASK(P_FIELD) ! IN Land mask PHIMH6A.53
PHIMH6A.54
PHIMH6A.55
REAL PHIMH6A.56
& RECIP_L_MO(P_FIELD) PHIMH6A.57
! ! IN Reciprocal of the Monin-Obukhov length (m^-1). PHIMH6A.58
&,Z_UV(P_FIELD) ! IN Height of wind level above roughness height (m PHIMH6A.59
&,Z_TQ(P_FIELD) ! IN Height of temperature, moisture and scalar lev PHIMH6A.60
! ! above the roughness height (m). PHIMH6A.61
&,Z0M(P_FIELD) ! IN Roughness length for momentum (m). PHIMH6A.62
&,Z0H(P_FIELD) ! IN Roughness length for heat/moisture/scalars (m) PHIMH6A.63
! PHIMH6A.64
REAL PHIMH6A.65
& PHI_M(P_FIELD)! OUT Stability function for momentum. PHIMH6A.66
&,PHI_H(P_FIELD)! OUT Stability function for heat/moisture/scalars. PHIMH6A.67
! PHIMH6A.68
!*L Workspace usage---------------------------------------------------- PHIMH6A.69
! No work areas are required. PHIMH6A.70
! PHIMH6A.71
!*---------------------------------------------------------------------- PHIMH6A.72
!*L External subprograms called: PHIMH6A.73
PHIMH6A.74
EXTERNAL TIMER PHIMH6A.75
PHIMH6A.76
!*---------------------------------------------------------------------- PHIMH6A.77
! Common and local physical constants. PHIMH6A.78
! PHIMH6A.79
! None. PHIMH6A.80
! PHIMH6A.81
! Define local variables. PHIMH6A.82
! PHIMH6A.83
INTEGER I,L ! Loop counter; horizontal field index. PHIMH6A.84
! PHIMH6A.85
REAL PHIMH6A.86
& PHI_MN ! Neutral value of stability function for momentum PHIMH6A.87
&,PHI_HN ! Neutral value of stability function for scalars. PHIMH6A.88
&,ZETA_UV ! Temporary in calculation of PHI_M. PHIMH6A.89
&,ZETA_0M ! Temporary in calculation of PHI_M. PHIMH6A.90
&,ZETA_TQ ! Temporary in calculation of PHI_H. PHIMH6A.91
&,ZETA_0H ! Temporary in calculation of PHI_H. PHIMH6A.92
&,X_UV_SQ ! Temporary in calculation of PHI_M. PHIMH6A.93
&,X_0M_SQ ! Temporary in calculation of PHI_M. PHIMH6A.94
&,X_UV ! Temporary in calculation of PHI_M. PHIMH6A.95
&,X_0M ! Temporary in calculation of PHI_M. PHIMH6A.96
&,Y_TQ ! Temporary in calculation of PHI_H. PHIMH6A.97
&,Y_0H ! Temporary in calculation of PHI_H. PHIMH6A.98
PHIMH6A.99
IF (LTIMER) THEN PHIMH6A.100
CALL TIMER
('PHI_M_H ',3) PHIMH6A.101
ENDIF PHIMH6A.102
PHIMH6A.103
I=0 ! Reset loop counter PHIMH6A.104
PHIMH6A.105
DO L=P1,P1+P_POINTS-1 PHIMH6A.106
! IF (L_LAND.AND.LAND_MASK(L)) THEN PHIMH6A.107
! I=I+1 PHIMH6A.108
! ELSEIF (.NOT.L_LAND) THEN PHIMH6A.109
I=L PHIMH6A.110
! ENDIF PHIMH6A.111
PHIMH6A.112
! IF(.NOT.L_LAND.OR.(L_LAND.AND.LAND_MASK(L))) THEN PHIMH6A.113
PHIMH6A.114
! PHIMH6A.115
!----------------------------------------------------------------------- PHIMH6A.116
!! 1. Calculate neutral values of PHI_M and PHI_H. PHIMH6A.117
!----------------------------------------------------------------------- PHIMH6A.118
! PHIMH6A.119
PHI_MN = LOG( (Z_UV(I) + Z0M(I)) / Z0M(I) ) PHIMH6A.120
PHI_HN = LOG( (Z_TQ(I) + Z0M(I)) / Z0H(I) ) PHIMH6A.121
! PHIMH6A.122
!----------------------------------------------------------------------- PHIMH6A.123
!! 2. Calculate stability parameters. PHIMH6A.124
!----------------------------------------------------------------------- PHIMH6A.125
! PHIMH6A.126
ZETA_UV = (Z_UV(I) + Z0M(I)) * RECIP_L_MO(I) PHIMH6A.127
ZETA_TQ = (Z_TQ(I) + Z0M(I)) * RECIP_L_MO(I) PHIMH6A.128
ZETA_0M = Z0M(I) * RECIP_L_MO(I) PHIMH6A.129
ZETA_0H = Z0H(I) * RECIP_L_MO(I) PHIMH6A.130
! PHIMH6A.131
!----------------------------------------------------------------------- PHIMH6A.132
!! 3. Calculate PHI_M and PHI_H for neutral and stable conditions. PHIMH6A.133
!----------------------------------------------------------------------- PHIMH6A.134
! PHIMH6A.135
IF (RECIP_L_MO(I) .GE. 0.0) THEN PHIMH6A.136
PHI_M(I) = PHI_MN + 4.0 * (ZETA_UV - ZETA_0M) PHIMH6A.137
PHI_H(I) = PHI_HN + PHIMH6A.138
& (1.0 + 2.0*ZETA_TQ) * (1.0 + 2.0*ZETA_TQ) - PHIMH6A.139
& (1.0 + 2.0*ZETA_0H) * (1.0 + 2.0*ZETA_0H) PHIMH6A.140
! PHIMH6A.141
!----------------------------------------------------------------------- PHIMH6A.142
!! 4. Calculate PHI_M and PHI_H for unstable conditions. PHIMH6A.143
!----------------------------------------------------------------------- PHIMH6A.144
! PHIMH6A.145
ELSE PHIMH6A.146
PHIMH6A.147
X_UV_SQ = SQRT(1.0 - 16.0*ZETA_UV) PHIMH6A.148
X_0M_SQ = SQRT(1.0 - 16.0*ZETA_0M) PHIMH6A.149
X_UV = SQRT(X_UV_SQ) PHIMH6A.150
X_0M = SQRT(X_0M_SQ) PHIMH6A.151
PHI_M(I) = PHI_MN - 2.0*LOG( (1.0+X_UV) / (1.0+X_0M) ) PHIMH6A.152
& - LOG( (1.0+X_UV_SQ) / (1.0+X_0M_SQ) ) PHIMH6A.153
& + 2.0*( ATAN(X_UV) - ATAN(X_0M) ) PHIMH6A.154
PHIMH6A.155
Y_TQ = SQRT(1.0 - 16.0*ZETA_TQ) PHIMH6A.156
Y_0H = SQRT(1.0 - 16.0*ZETA_0H) PHIMH6A.157
PHI_H(I) = PHI_HN - 2.0*LOG( (1.0+Y_TQ) / (1.0+Y_0H) ) PHIMH6A.158
PHIMH6A.159
ENDIF PHIMH6A.160
PHIMH6A.161
! ENDIF ! L_LAND and LAND_MASK PHIMH6A.162
PHIMH6A.163
ENDDO PHIMH6A.164
PHIMH6A.165
IF (LTIMER) THEN PHIMH6A.166
CALL TIMER('PHI_M_H ',4) PHIMH6A.167
ENDIF PHIMH6A.168
PHIMH6A.169
RETURN PHIMH6A.170
END PHIMH6A.171
*ENDIF PHIMH6A.172