*IF DEF,A03_3A ASJ4F401.11
C ******************************COPYRIGHT****************************** GTS2F400.8839
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.8840
C GTS2F400.8841
C Use, duplication or disclosure of this code is subject to the GTS2F400.8842
C restrictions as set forth in the contract. GTS2F400.8843
C GTS2F400.8844
C Meteorological Office GTS2F400.8845
C London Road GTS2F400.8846
C BRACKNELL GTS2F400.8847
C Berkshire UK GTS2F400.8848
C RG12 2SZ GTS2F400.8849
C GTS2F400.8850
C If no contract has been raised with this copy of the code, the use, GTS2F400.8851
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.8852
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.8853
C Modelling at the above address. GTS2F400.8854
C ******************************COPYRIGHT****************************** GTS2F400.8855
C GTS2F400.8856
C*LL SUBROUTINE SF_EXCH------------------------------------------------ SFEXCH3A.3
CLL SFEXCH3A.4
CLL Purpose: Calculate coefficients of turbulent exchange between SFEXCH3A.5
CLL the surface and the lowest atmospheric layer, and SFEXCH3A.6
CLL "explicit" fluxes between the surface and this layer. SFEXCH3A.7
CLL SFEXCH3A.8
CLL Suitable for Single Column use if *IF definition IBM is selected. SFEXCH3A.9
CLL SFEXCH3A.10
CLL Canopy evaporation made implicit SFEXCH3A.11
CLL with respect to canopy water content (requiring TIMESTEP to be SFEXCH3A.12
CLL passed in). SFEXCH3A.13
CLL SFEXCH3A.14
CLL CW, FH <- programmer of some or all of previous code or changes SFEXCH3A.15
CLL SFEXCH3A.16
CLL Model Modification history: SFEXCH3A.17
CLL version Date SFEXCH3A.18
CLL 3.4 18/10/94 *DECK inserted into UM version 3.4. S Jackson SFEXCH3A.19
CLL SFEXCH3A.20
CLL 4.0 05/01/95 rhostar*cD*vshr before interpolation output ASJ1F400.22
CLL via argument list. R.N.B.Smith ASJ1F400.23
CLL ASJ1F400.24
CLL 4.0 30/12/94 Changes to formulation of cH and 10m wind ARN1F400.7
CLL diagnostic; z0h set to z0(veg). ARN1F400.8
CLL R.N.B.Smith ARN1F400.9
CLL 4.1 08/05/96 decks A03_2C and A03_3B removed ASJ4F401.12
CLL S D Jackson ASJ4F401.13
CLL ARN1F400.10
CLL 4.1 08/05/96 Logical switch for rapidly mixing boundary layer ASJ3F401.10
CLL S D Jackson ASJ3F401.11
CLL ASJ3F401.12
!LL 4.1 12/7/96 Added MPP code P.Burton APBGF401.94
CLL 4.2 Oct. 96 T3E migration - *DEF CRAY removed, WHENIMD removed GSS1F402.70
CLL S J Swarbrick GSS1F402.71
!LL 4.3 14/01/97 MPP code : Corrected setting of polar rows GPB1F403.5
!LL P.Burton GPB1F403.6
CLL 4.3 09/06/97 Add swapbounds for RHOKM_1 & CDR10M. ASJ1F403.6
CLL D.Sexton/RTHBarnes ASJ1F403.7
CLL 4.4 08/09/97 L_BL_LSPICE specifies mixed phase precipitation ADM3F404.219
CLL scheme. D. Wilson ADM3F404.220
CLL 4.5 Jul. 98 Kill the IBM specific lines (JCThil) AJC1F405.110
CLL Programming standard: Unified Model Documentation Paper No 4, SFEXCH3A.21
CLL Version 2, dated 18/1/90. SFEXCH3A.22
CLL SFEXCH3A.23
CLL System component covered: Part of P243. SFEXCH3A.24
CLL SFEXCH3A.25
CLL Project task: SFEXCH3A.26
CLL SFEXCH3A.27
CLL Documentation: UM Documentation Paper No 24, section P243. SFEXCH3A.28
CLL See especially sub-section (ix). SFEXCH3A.29
CLL SFEXCH3A.30
CLLEND------------------------------------------------------------------ SFEXCH3A.31
C* SFEXCH3A.32
C*L Arguments --------------------------------------------------------- SFEXCH3A.33
SUBROUTINE SF_EXCH ( 4,99SFEXCH3A.34
+ P_POINTS,LAND_PTS,U_POINTS,ROW_LENGTH,P_ROWS,U_ROWS SFEXCH3A.35
+,LAND_INDEX,P1,GATHER SFEXCH3A.37
+,AK_1,BK_1,TIMESTEP SFEXCH3A.39
+,CANOPY,CATCH,CF_1,ICE_FRACT,LAND_MASK,PSTAR,Q_1 SFEXCH3A.40
+,QCF_1,QCL_1,RESIST,ROOTD,SMC,SMVCCL,SMVCWT,LYING_SNOW,T_1,TSTAR SFEXCH3A.41
+,U_1,V_1,U_1_P,V_1_P,U_0,V_0,Z0V,SIL_OROG,Z1,Z0MSEA,HO2R2_OROG SFEXCH3A.42
&,BQ_1,BT_1,BF_1,CD,CH,EA,ES,ESL,FQW_1,FQW_LEAD,FTL_1,FTL_LEAD ADM3F404.221
+,TAUX_1,TAUY_1,QW_1,RHOKE,RHOKEA,RHOKES,RHOKESL,RHOKH_1,RHOKM_1 SFEXCH3A.44
+,RIB,TL_1,TSTAR_NL,VSHR,Z0H,Z0M,Z0M_EFF,H_BLEND ARN1F400.11
&,T1_SD,Q1_SD SFEXCH3A.46
&,RHO_CD_MODV1 ASJ1F400.25
+,CDR10M,CHR1P5M,CER1P5M,FME SFEXCH3A.47
+,SU10,SV10,SQ1P5,ST1P5,SFME SFEXCH3A.48
+,NRML SFEXCH3A.49
&,L_Z0_OROG,L_RMBL,L_BL_LSPICE,ERROR,LTIMER ADM3F404.222
*IF DEF,SCMA AJC0F405.77
& ,FACTOR_RHOKH,OBS AJC0F405.78
*ENDIF AJC0F405.79
+) SFEXCH3A.51
IMPLICIT NONE SFEXCH3A.52
C SFEXCH3A.53
C Input variables. All fields are on P grid except where noted. SFEXCH3A.54
C Fxxx in a comment indicates the file from which the data are taken. SFEXCH3A.55
C SFEXCH3A.56
C SFEXCH3A.58
C GENERAL NOTES ABOUT GRID-DEFINITION INPUT VARIABLES. SFEXCH3A.59
C ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SFEXCH3A.60
C For global data :- SFEXCH3A.61
C SFEXCH3A.62
C An Arakawa B-grid is assumed in which each pole is represented by a SFEXCH3A.63
C row of P-grid points. These polar rows are omitted in the input and SFEXCH3A.64
C output of the present subroutine, so that the argument P_ROWS is two SFEXCH3A.65
C less than the total number of P-rows in the grid. Land specific SFEXCH3A.66
C variables that are required as INput by the higher level routine, SFEXCH3A.67
C BDY_LAYR, are stored on P_grid land points only and land pts on polar SFEXCH3A.68
C rows are not input or output by this routine ; diagnostic variables SFEXCH3A.69
C must be defined on land and sea points for post processing. SFEXCH3A.70
C If defined variable IBM is selected then land point calculations are SFEXCH3A.71
C performed using the array LAND_INDEX to select land points. But note SFEXCH3A.72
C that elements of LAND_INDEX define land points on the full field SFEXCH3A.73
C (ie including polar rows). SFEXCH3A.74
C SFEXCH3A.75
C Entire fields of UV-grid values are taken as input, but the two SFEXCH3A.76
C polemost rows are (a) not updated, in the case of INOUT fields, or SFEXCH3A.77
C (b) set to zero, in the case of OUT fields. SFEXCH3A.78
C SFEXCH3A.79
C For limited-area data :- SFEXCH3A.80
C SFEXCH3A.81
C The above applies, but for "polar rows", etc., read "rows at the SFEXCH3A.82
C north and south boundaries of the area", etc. E.g. if you want to SFEXCH3A.83
C do calculations in UV-rows n to m inclusive, the input data will be SFEXCH3A.84
C on P-rows n to m+1, and UV-rows n-1 to m+1. P-rows n to m will SFEXCH3A.85
C then be updated. Land specific variables are processed as for global SFEXCH3A.86
C data. SFEXCH3A.87
C SFEXCH3A.88
C For both cases, the following equalities apply amongst the input SFEXCH3A.89
C grid-definition variables :- SFEXCH3A.90
C SFEXCH3A.91
C P_POINTS = P_ROWS * ROW_LENGTH SFEXCH3A.92
C U_POINTS = U_ROWS * ROW_LENGTH SFEXCH3A.93
C U_ROWS = P_ROWS + 1 SFEXCH3A.94
C LAND_PTS <= P_POINTS SFEXCH3A.95
C SFEXCH3A.96
C An error condition is returned if the input variables don't satisfy SFEXCH3A.97
C these equalities. (There is of course redundancy here; a compromise SFEXCH3A.98
C between economy, clarity and easy dimensioning is intended.) SFEXCH3A.99
C SFEXCH3A.100
C NB: All this has severe implications for batching/macrotasking; SFEXCH3A.101
C effectively it can't be done on a shared-memory machine without SFEXCH3A.102
C either rewriting this routine or using expensive synchronizations SFEXCH3A.103
C (or other messy and/or undesirable subterfuges). SFEXCH3A.104
C SFEXCH3A.105
C SFEXCH3A.120
LOGICAL LTIMER SFEXCH3A.121
*IF DEF,SCMA AJC0F405.80
LOGICAL OBS ! Switch for OBS forcing AJC0F405.81
*ENDIF AJC0F405.82
C SFEXCH3A.122
INTEGER ! Variables defining grid. SFEXCH3A.123
+ P_POINTS ! IN Number of P-grid points to be processed. SFEXCH3A.124
+,LAND_PTS ! IN Number of land points to be processed. SFEXCH3A.125
+,U_POINTS ! IN Number of UV-grid points. SFEXCH3A.126
+,ROW_LENGTH ! IN No. of points in latitude row (inclusive SFEXCH3A.130
C ! of endpoints for ltd. area model). SFEXCH3A.131
+,P_ROWS ! IN Number of rows of data on P-grid. SFEXCH3A.135
+,U_ROWS ! IN Number of rows of data on UV-grid. SFEXCH3A.139
+,LAND_INDEX(LAND_PTS)! IN Index for compressed land point array; SFEXCH3A.143
C ! ith element holds position in the FULL SFEXCH3A.144
C ! field of the ith land pt to be processed SFEXCH3A.145
+,P1 ! IN First P-point to be processed. SFEXCH3A.146
LOGICAL SFEXCH3A.147
+ GATHER ! IN If true then leads variables are comp- SFEXCH3A.148
C ! ressed for sea-ice calculations. This SFEXCH3A.149
C ! saves duplicating calculations if there SFEXCH3A.150
C ! are a relatively few of sea-ice points. SFEXCH3A.151
C ! Set to false for a limited area run SFEXCH3A.152
C ! with a high proportion of sea-ice. SFEXCH3A.153
REAL SFEXCH3A.155
+ AK_1 ! IN Hybrid "A" for lowest model layer. SFEXCH3A.156
+,BK_1 ! IN Hybrid "B" for lowest model layer. SFEXCH3A.157
+,TIMESTEP ! IN Timestep in seconds for EPDT calc. SFEXCH3A.158
+,CANOPY(LAND_PTS) ! IN Surface water (kg per sq metre). F642. SFEXCH3A.159
+,CATCH(LAND_PTS) ! IN Surface capacity (max. surface water) SFEXCH3A.160
C ! (kg per sq metre). F6416. SFEXCH3A.161
+,CF_1(P_POINTS) ! IN Cloud fraction for lowest atmospheric SFEXCH3A.162
C ! layer (decimal fraction). SFEXCH3A.163
*IF DEF,SCMA AJC0F405.83
& ,FACTOR_RHOKH(P_FIELD) ! IN Factor for modifying surface AJC0F405.84
! fluxes if OBS forcing used AJC0F405.85
AJC0F405.86
C AJC0F405.87
C If OBS run use RHOKH_1 and FACTOR_RHOKH AJC0F405.88
C (from FLUX_H and FLUX_E input by namelist) AJC0F405.89
C AJC0F405.90
If (.NOT.OBS) AJC0F405.91
& RHOKH_1(I) = RHOSTAR * CH(I) * VSHR(I) AJC0F405.92
*ENDIF AJC0F405.93
+,ICE_FRACT(P_POINTS) ! IN Fraction of gridbox which is sea-ice. SFEXCH3A.164
+,LYING_SNOW(P_POINTS)! IN Lying snow amount (kg per sq metre). SFEXCH3A.165
+,PSTAR(P_POINTS) ! IN Surface pressure (Pascals). SFEXCH3A.166
+,Q_1(P_POINTS) ! IN Specific humidity for lowest atmospheric SFEXCH3A.167
C ! layer (kg water per kg air). SFEXCH3A.168
+,QCF_1(P_POINTS) ! IN Cloud ice for lowest atmospheric layer SFEXCH3A.169
C ! (kg water per kg air). SFEXCH3A.170
+,QCL_1(P_POINTS) ! IN Cloud liquid water for lowest atm layer SFEXCH3A.171
C ! (kg water per kg air). SFEXCH3A.172
+,RESIST(LAND_PTS) ! IN "Stomatal" resistance to evaporation SFEXCH3A.173
C ! (seconds per metre). F6415. SFEXCH3A.174
+,ROOTD(LAND_PTS) ! IN "Root depth" (metres). F6412. SFEXCH3A.175
+,SMC(LAND_PTS) ! IN Soil moisture content (kg per sq m). SFEXCH3A.176
C ! F621. SFEXCH3A.177
+,SMVCCL(LAND_PTS) ! IN Critical volumetric SMC (cubic metres SFEXCH3A.178
C ! per cubic metre of soil). F6232. SFEXCH3A.179
+,SMVCWT(LAND_PTS) ! IN Volumetric wilting point (cubic m of SFEXCH3A.180
C ! water per cubic m of soil). F6231. SFEXCH3A.181
C SFEXCH3A.182
C Note: (SMVCCL - SMVCWT) is the critical volumetric available soil SFEXCH3A.183
C moisture content. ~~~~~~~~~ SFEXCH3A.184
C SFEXCH3A.185
REAL ! (Split to avoid > 19 continuations.) SFEXCH3A.186
+ T_1(P_POINTS) ! IN Temperature for lowest atmospheric layer SFEXCH3A.187
C ! (Kelvin). SFEXCH3A.188
+,TSTAR(P_POINTS) ! IN Mean gridsquare surface temperature (K). SFEXCH3A.189
+,U_1(U_POINTS) ! IN West-to-east wind component for lowest SFEXCH3A.190
C ! atmospheric layer (m/s). On UV grid. SFEXCH3A.191
+,V_1(U_POINTS) ! IN South-to-north wind component for lowest SFEXCH3A.192
C ! atmospheric layer (m/s). On UV grid. SFEXCH3A.193
+,U_1_P(P_POINTS) ! IN West-to-east wind component for lowest SFEXCH3A.194
C ! atmospheric layer (m/s). On P grid. SFEXCH3A.195
+,V_1_P(P_POINTS) ! IN South-to-north wind component for lowest SFEXCH3A.199
C ! atmospheric layer (m/s). On P grid. SFEXCH3A.200
+,U_0(U_POINTS) ! IN West-to-east component of ocean surface SFEXCH3A.204
C ! current (m/s; ASSUMED zero over land). SFEXCH3A.205
C ! UV grid. F615. SFEXCH3A.206
+,V_0(U_POINTS) ! IN South-to-north component of ocean surface SFEXCH3A.207
C ! current (m/s; ASSUMED zero over land). SFEXCH3A.208
C ! UV grid. F616. SFEXCH3A.209
+,Z0V(P_POINTS) ! IN Vegetative roughness length (m). F6418. SFEXCH3A.210
+,SIL_OROG(LAND_PTS) ! IN Silhouette area of unresolved orography SFEXCH3A.211
C ! per unit horizontal area SFEXCH3A.212
+,Z1(P_POINTS) ! IN Height of lowest atmospheric level (m). SFEXCH3A.213
&,HO2R2_OROG(LAND_PTS) ! IN Peak to trough height of unresolved SFEXCH3A.214
C ! orography devided by 2SQRT(2) (m). SFEXCH3A.215
LOGICAL SFEXCH3A.216
+ LAND_MASK(P_POINTS) ! IN .TRUE. for land; .FALSE. elsewhere. F60. SFEXCH3A.217
+,SU10 ! IN STASH flag for 10-metre W wind. SFEXCH3A.218
+,SV10 ! IN STASH flag for 10-metre S wind. SFEXCH3A.219
+,SQ1P5 ! IN STASH flag for 1.5-metre sp humidity. SFEXCH3A.220
+,ST1P5 ! IN STASH flag for 1.5-metre temperature. SFEXCH3A.221
+,SFME ! IN STASH flag for wind mixing energy flux. SFEXCH3A.222
+,L_Z0_OROG ! IN .TRUE. to use orographic roughness. SFEXCH3A.223
+,L_RMBL ! IN .TRUE. to use rapidly mixing boundary ASJ3F401.13
C ! scheme ASJ3F401.14
&,L_BL_LSPICE ! IN ADM3F404.223
! TRUE Use scientific treatment of mixed ADM3F404.224
! phase precip scheme. ADM3F404.225
! FALSE Do not use mixed phase precip ADM3F404.226
! considerations ADM3F404.227
C SFEXCH3A.224
C Modified (INOUT) variables. SFEXCH3A.225
C SFEXCH3A.226
REAL SFEXCH3A.227
+ Z0MSEA(P_POINTS) ! INOUT Sea-surface roughness length for SFEXCH3A.228
C ! momentum (m). F617. SFEXCH3A.229
C SFEXCH3A.230
C Output variables. SFEXCH3A.231
C SFEXCH3A.232
REAL SFEXCH3A.233
+ BQ_1(P_POINTS) ! OUT A buoyancy parameter for lowest atm level SFEXCH3A.234
C ! ("beta-q twiddle"). SFEXCH3A.235
+,BT_1(P_POINTS) ! OUT A buoyancy parameter for lowest atm level. SFEXCH3A.236
C ! ("beta-T twiddle"). SFEXCH3A.237
&,BF_1(P_POINTS) ADM3F404.228
! OUT A buoyancy parameter for lowest atm level. ADM3F404.229
! ("beta-F twiddle"). ADM3F404.230
+,CD(P_POINTS) ! OUT Bulk transfer coefficient for momentum. SFEXCH3A.238
+,CH(P_POINTS) ! OUT Bulk transfer coefficient for heat and/or SFEXCH3A.239
C ! moisture. SFEXCH3A.240
+,CDR10M(U_POINTS) ! OUT Reqd for calculation of 10m wind (u & v). SFEXCH3A.241
C ! NBB: This is output on the UV-grid, but SFEXCH3A.242
C ! with the first and last rows set to a SFEXCH3A.243
C ! "missing data indicator". SFEXCH3A.244
C ! Sea-ice leads ignored. See 3.D.7 below. SFEXCH3A.245
+,CHR1P5M(P_POINTS) ! OUT Reqd for calculation of 1.5m temperature. SFEXCH3A.246
C ! Sea-ice leads ignored. See 3.D.7 below. SFEXCH3A.247
+,CER1P5M(P_POINTS) ! OUT Reqd for calculation of 1.5m sp humidity. SFEXCH3A.248
C ! Sea-ice leads ignored. See 3.D.7 below. SFEXCH3A.249
&,RHO_CD_MODV1(P_POINTS) ASJ1F400.26
C ! OUT rhostar*cD*vshr before horizontal ASJ1F400.27
C ! interpolation output as a diagnostic. ASJ1F400.28
+,EA(P_POINTS) ! OUT "Explicit" evaporation with only aero- SFEXCH3A.250
C ! dynamic resistance (+ve), or condensation SFEXCH3A.251
C ! (-ve), averaged over grid box (kg/m2/s). SFEXCH3A.252
+,ES(P_POINTS) ! OUT "Explicit" evapotranspiration through a SFEXCH3A.253
C ! non-aerodynamic resistance. Always non- SFEXCH3A.254
C ! negative (kg/m2/s). SFEXCH3A.255
+,ESL(P_POINTS) ! OUT ES without fractional weighting factor SFEXCH3A.256
C ! "L" is for "local". SFEXCH3A.257
REAL ! (Split to avoid > 19 continuations.) SFEXCH3A.258
+ FQW_1(P_POINTS) ! OUT "Explicit" surface flux of QW (i.e. SFEXCH3A.259
C ! evaporation), on P-grid (kg/m2/s). SFEXCH3A.260
+,FQW_LEAD(P_POINTS) ! OUT FQW for leads fraction of gridsquare. SFEXCH3A.261
C ! Missing data indicator at non sea-ice. SFEXCH3A.262
+,FTL_1(P_POINTS) ! OUT "Explicit" surface flux of TL = H/CP. SFEXCH3A.263
C ! (sensible heat / CP). SFEXCH3A.264
+,FTL_LEAD(P_POINTS) ! OUT FTL for leads fraction of gridsquare, SFEXCH3A.265
C ! Missing data indicator at non sea-ice. SFEXCH3A.266
+,TAUX_1(U_POINTS) ! OUT "Explicit" x-component of surface SFEXCH3A.267
C ! turbulent stress; on UV-grid; first and SFEXCH3A.268
C ! last rows set to a "missing data SFEXCH3A.269
C ! indicator". (Newtons per square metre) SFEXCH3A.270
+,TAUY_1(U_POINTS) ! OUT "Explicit" y-component of surface SFEXCH3A.271
C ! turbulent stress; on UV-grid; first and SFEXCH3A.272
C ! last rows set to a "missing data SFEXCH3A.273
C ! indicator". (Newtons per square metre) SFEXCH3A.274
+,QW_1(P_POINTS) ! OUT Total water content of lowest SFEXCH3A.275
C ! atmospheric layer (kg per kg air). SFEXCH3A.276
C SFEXCH3A.277
REAL ! Surface exchange coefficients;passed to subroutine IMPL_CAL SFEXCH3A.278
+ RHOKE(P_POINTS) ! OUT For FQW, then *GAMMA(1) for implicit calcs SFEXCH3A.279
+,RHOKEA(P_POINTS) ! OUT For EA, then *GAMMA(1) for implicit calcs. SFEXCH3A.280
+,RHOKES(P_POINTS) ! OUT For ES, then *GAMMA(1) for implicit calcs SFEXCH3A.281
+,RHOKESL(P_POINTS) ! OUT For ESL,then *GAMMA(1) for implicit calcs SFEXCH3A.282
+,RHOKH_1(P_POINTS) ! OUT For FTL,then *GAMMA(1) for implicit calcs SFEXCH3A.283
+,RHOKM_1(U_POINTS) ! OUT For momentum, then *GAMMA(1) for implicit SFEXCH3A.284
C ! calculations. NBB: This is output on the SFEXCH3A.285
C ! UV-grid, but with the first and last SFEXCH3A.286
C ! rows set to a "missing data indicator". SFEXCH3A.287
REAL ! End of surface exchange coefficients. SFEXCH3A.288
+ RIB(P_POINTS) ! OUT Bulk Richardson number for lowest layer. SFEXCH3A.289
+,TL_1(P_POINTS) ! OUT Liquid/frozen water temperature for SFEXCH3A.290
C ! lowest atmospheric layer (K). SFEXCH3A.291
+,TSTAR_NL(P_POINTS) ! OUT TSTAR No Leads: surface temperature SFEXCH3A.292
C ! over sea-ice fraction of gridsquare. SFEXCH3A.293
C ! =TSTAR over non sea-ice points. SFEXCH3A.294
+,VSHR(P_POINTS) ! OUT Magnitude of surface-to-lowest-lev. wind SFEXCH3A.295
+,Z0H(P_POINTS) ! OUT Roughness length for heat and moisture m SFEXCH3A.296
+,Z0M(P_POINTS) ! OUT Roughness length for momentum (m). SFEXCH3A.297
+,Z0M_EFF(P_POINTS) ! OUT Effective roughness length for momentum SFEXCH3A.298
+,H_BLEND(P_POINTS) ! OUT Blending height SFEXCH3A.301
&,T1_SD(P_POINTS) ! OUT Standard deviation of turbulent SFEXCH3A.302
C ! fluctuations of surface layer SFEXCH3A.303
C ! temperature (K). SFEXCH3A.304
&,Q1_SD(P_POINTS) ! OUT Standard deviation of turbulent SFEXCH3A.305
C ! fluctuations of surface layer SFEXCH3A.306
C ! specific humidity (kg/kg). SFEXCH3A.307
+,FME(P_POINTS) ! OUT Wind mixing energy flux (Watts/sq m). SFEXCH3A.308
INTEGER SFEXCH3A.309
+ NRML(P_POINTS) ! OUT 1 if surface layer unstable, else 0. SFEXCH3A.310
+,ERROR ! OUT 1 if grid definition faulty; else 0. SFEXCH3A.311
C* SFEXCH3A.312
C*L Symbolic constants ------------------------------------------------ SFEXCH3A.313
C SFEXCH3A.314
C (1) UM-wide common parameters. SFEXCH3A.315
C SFEXCH3A.316
*CALL C_0_DG_C 
SFEXCH3A.317
*CALL C_LHEAT SFEXCH3A.318
*CALL C_G SFEXCH3A.319
*CALL C_R_CP SFEXCH3A.320
*CALL C_EPSLON SFEXCH3A.321
*CALL C_VKMAN SFEXCH3A.322
C SFEXCH3A.323
C (2) Boundary Layer local parameters. SFEXCH3A.324
C SFEXCH3A.325
*CALL C_CHARNK SFEXCH3A.326
*CALL C_DENSTY SFEXCH3A.327
*CALL C_GAMMA SFEXCH3A.328
*CALL C_HT_M SFEXCH3A.329
*CALL C_ROUGH SFEXCH3A.330
*CALL C_SURF SFEXCH3A.331
*CALL C_MDI SFEXCH3A.332
C SFEXCH3A.333
C (3) Derived local parameters. SFEXCH3A.334
C SFEXCH3A.335
REAL ETAR,GRCP,LCRCP,LFRCP,LS,LSRCP,H_BLEND_MIN,H_BLEND_MAX ARN1F400.12
PARAMETER ( SFEXCH3A.337
+ ETAR=1./(1.-EPSILON) ! Used in calc of buoyancy parameter BETAC. SFEXCH3A.338
+,GRCP=G/CP ! Used in calc of dT across surface layer. SFEXCH3A.339
+,LCRCP=LC/CP ! Evaporation-to-dT conversion factor. SFEXCH3A.340
+,LFRCP=LF/CP ! Freezing-to-dT conversion factor. SFEXCH3A.341
+,LS=LF+LC ! Latent heat of sublimation. SFEXCH3A.342
+,LSRCP=LS/CP ! Sublimation-to-dT conversion factor. SFEXCH3A.343
+,H_BLEND_MIN=0.0 ! Minimum blending height. SFEXCH3A.344
&,H_BLEND_MAX=1000.0 ! Maximum blending height (m). ARN1F400.13
+) SFEXCH3A.345
C* SFEXCH3A.346
*IF DEF,MPP GPB1F403.1
! MPP Common block GPB1F403.2
*CALL PARVARS GPB1F403.3
*ENDIF GPB1F403.4
C*L SFEXCH3A.347
C External subprograms called. SFEXCH3A.348
C SFEXCH3A.349
EXTERNAL FCDCH,QSAT,QSAT_WAT,SFL_INT ADM3F404.231
*IF -DEF,SCMA AJC1F405.111
EXTERNAL P_TO_UV,UV_TO_P GSS1F403.52
*ENDIF SFEXCH3A.353
EXTERNAL TIMER SFEXCH3A.354
C* SFEXCH3A.355
C SFEXCH3A.356
C Define local storage. SFEXCH3A.357
C SFEXCH3A.358
C (a) Workspace. SFEXCH3A.359
C SFEXCH3A.360
C*L Workspace --------------------------------------------------------- SFEXCH3A.361
C 25 blocks of real workspace are required, as follows. SFEXCH3A.363
REAL SFEXCH3A.364
+ CD_LEAD(P_POINTS) ! Bulk transfer coefficient for momentum SFEXCH3A.365
C ! over sea-ice leads.Missing data over non SFEXCH3A.366
C ! sea-ice points.(Temporary store for Z0MIZ) SFEXCH3A.367
+,CD_MIZ(P_POINTS) ! Bulk transfer coefficient for momentum SFEXCH3A.368
C ! over the sea-ice Marginal Ice Zone. SFEXCH3A.369
C ! Missing data indicator over non sea-ice. SFEXCH3A.370
+,CH_LEAD(P_POINTS) ! Bulk transfer coefficient for heat and SFEXCH3A.371
C ! or moisture over sea ice leads. SFEXCH3A.372
C ! Missing data indicator over non sea-ice. SFEXCH3A.373
+,CH_MIZ(P_POINTS) ! Bulk transfer coefficient for heat and SFEXCH3A.374
C ! or moisture over the Marginal Ice Zone. SFEXCH3A.375
C ! Missing data indicator over non sea-ice. SFEXCH3A.376
+,CD_STD(P_POINTS) ! Local drag coefficient for SFEXCH3A.377
C ! calculation of interpolation coefficients SFEXCH3A.378
+,DQ(P_POINTS) ! Sp humidity difference between surface SFEXCH3A.385
C ! and lowest atmospheric level (Q1 - Q*). SFEXCH3A.386
C ! Holds value over sea-ice where ICE_FRACT SFEXCH3A.387
C ! >0 i.e. Leads contribution not included. SFEXCH3A.388
&,DQI(P_POINTS) ADM3F404.232
! Ice water difference between surface ADM3F404.233
! and lowest atmospheric level (Q1 - Q*). ADM3F404.234
! Holds value over sea-ice where ICE_FRACT ADM3F404.235
! >0 i.e. Leads contribution not included. ADM3F404.236
+,DQ_LEAD(P_POINTS) ! DQ for leads fraction of gridsquare. SFEXCH3A.389
C ! Missing data indicator over non sea-ice. SFEXCH3A.390
&,DQI_LEAD(P_POINTS) ADM3F404.237
! DQI for leads fraction of gridsquare. ADM3F404.238
! Missing data indicator over non sea-ice. ADM3F404.239
+,DTEMP(P_POINTS) ! Liquid/ice static energy difference SFEXCH3A.391
C ! between surface and lowest atmospheric SFEXCH3A.392
C ! level, divided by CP (a modified SFEXCH3A.393
C ! temperature difference). SFEXCH3A.394
C ! Holds value over sea-ice where ICE_FRACT SFEXCH3A.395
C ! >0 i.e. Leads contribution not included. SFEXCH3A.396
+,DTEMP_LEAD(P_POINTS) ! DTEMP for leads fraction of gridsquare. SFEXCH3A.397
C ! Missing data indicator over non sea-ice. SFEXCH3A.398
+,FRACA(P_POINTS) ! Fraction of surface moisture flux with SFEXCH3A.399
C ! only aerodynamic resistance. SFEXCH3A.400
+,PSIS(P_POINTS) ! Soil moisture availability factor. SFEXCH3A.401
+,QSL(P_POINTS) ! Saturated sp humidity at liquid/ice SFEXCH3A.402
C ! temperature and pressure of lowest SFEXCH3A.403
C ! atmospheric level. SFEXCH3A.404
+,QSTAR(P_POINTS) ! Surface saturated sp humidity. Holds SFEXCH3A.405
C ! value over sea-ice where ICE_FRACT > 0. SFEXCH3A.406
C ! i.e. Leads contribution not included. SFEXCH3A.407
+,QSTAR_LEAD(P_POINTS) ! QSTAR for sea-ice leads. SFEXCH3A.408
C ! Missing data indicator over non sea-ice. SFEXCH3A.409
+,RESFS(P_POINTS) ! Combined soil, stomatal and aerodynamic SFEXCH3A.410
C ! resistance factor = PSIS/(1+RS/RA) for SFEXCH3A.411
C ! fraction 1 - FRACA of the gridbox SFEXCH3A.412
+,RESFT(P_POINTS) ! Total resistance factor for moisture SFEXCH3A.413
C ! transfer to/from the surface SFEXCH3A.414
C ! = FRACA + (1-FRACA)*RESFS so that SFEXCH3A.415
C ! 1/RT = RESFT/RA SFEXCH3A.416
+,RIB_LEAD(P_POINTS) ! Bulk Richardson no. for sea-ice leads at SFEXCH3A.417
C ! lowest layer. At non sea-ice points holds SFEXCH3A.418
C ! RIB for FCDCH calculation, then set to SFEXCH3A.419
C ! to missing data indicator. SFEXCH3A.420
+,RS(P_POINTS) ! Stomatal resistance to evaporation. SFEXCH3A.421
+,U_0_P(P_POINTS) ! West-to-east component of ocean surface SFEXCH3A.422
C ! current (m/s; zero over land if U_0 OK). SFEXCH3A.423
C ! P grid. F615. SFEXCH3A.424
+,V_0_P(P_POINTS) ! South-to-north component of ocean surface SFEXCH3A.425
C ! current (m/s; zero over land if V_0 OK). SFEXCH3A.426
C ! P grid. F616. SFEXCH3A.427
+,Z0F(P_POINTS) ! Roughness length for free-convective heat SFEXCH3A.428
C ! and moisture transport. SFEXCH3A.429
+,Z0FS(P_POINTS) ! Roughness length for free-convective heat SFEXCH3A.430
C ! and moisture transport over sea. SFEXCH3A.431
+,Z0HS(P_POINTS) ! Roughness length for heat and moisture SFEXCH3A.432
C ! transport over sea. SFEXCH3A.433
&,WIND_PROFILE_FACTOR(P_POINTS) ARN1F400.14
C ! For transforming effective surface transfer ARN1F400.15
C ! coefficients to those excluding form drag. ARN1F400.16
C SFEXCH3A.434
C Workspace (reqd for compression). SFEXCH3A.435
INTEGER SFEXCH3A.436
+ SICE_INDEX(P_POINTS) ! Index vector for gather to sea-ice points SFEXCH3A.437
LOGICAL ITEST(P_POINTS) ! Used as 'logical' for compression. SFEXCH3A.438
C* SFEXCH3A.503
C SFEXCH3A.504
C (b) Scalars. SFEXCH3A.505
C SFEXCH3A.506
INTEGER SFEXCH3A.507
+ I ! Loop counter (horizontal field index). SFEXCH3A.508
+,J ! Offset counter within I-loop. SFEXCH3A.509
+,L ! Loop counter (land point field index). SFEXCH3A.510
+,NSICE ! Number of sea-ice points. SFEXCH3A.511
+,SI ! Loop counter (sea-ice field index). SFEXCH3A.515
REAL SFEXCH3A.516
+ AL ! Temporary in calculation of buoyancy parameters. SFEXCH3A.517
+,ALPHAL ! Temporary in calculation of buoyancy parameters. SFEXCH3A.518
+,BETAC ! Temporary in calculation of buoyancy parameters. SFEXCH3A.519
+,BETAQ ! Temporary in calculation of buoyancy parameters. SFEXCH3A.520
+,BETAT ! Temporary in calculation of buoyancy parameters. SFEXCH3A.521
+,CHN ! Neutral-stability value of CH, used as a first SFEXCH3A.522
C ! approximation to the "true" CH. SFEXCH3A.523
REAL ! (Split to avoid having > 19 continuation lines.) SFEXCH3A.524
+ RHOSTAR ! Surface air density in kg per cubic metre. SFEXCH3A.525
+,TAU ! Magnitude of surface wind stress over sea. SFEXCH3A.526
+,SMCRIT ! "Critical" available SMC in kg per sq m. SFEXCH3A.527
+,USHEAR ! U-component of surface-to-lowest-level wind shear. SFEXCH3A.528
+,VSHEAR ! V-component of surface-to-lowest-level wind shear. SFEXCH3A.529
+,VSHR2 ! Square of magnitude of surface-to-lowest-level SFEXCH3A.530
C ! wind shear. SFEXCH3A.531
+,ZETAM ! Temporary in calculation of CHN. SFEXCH3A.532
+,ZETAH ! Temporary in calculation of CHN. SFEXCH3A.533
+,ZETA1 ! Temporary in calculation of land roughness SFEXCH3A.534
+,ZETA2 ! Temporary in calculation of land roughness SFEXCH3A.535
+,ZETA3 ! Temporary in calculation of land roughness SFEXCH3A.536
+,Z0 ! Roughness length over land. SFEXCH3A.537
+,EPDT ! "Potential" Evaporation * Timestep SFEXCH3A.538
&,VS ! Surface layer friction velocity SFEXCH3A.539
&,VSF1_CUBED ! Cube of surface layer free convective scaling SFEXCH3A.540
C ! velocity SFEXCH3A.541
&,WS1 ! Turbulent velocity scale for surface layer SFEXCH3A.542
C SFEXCH3A.543
C----------------------------------------------------------------------- SFEXCH3A.544
CL 0. Check that the scalars input to define the grid are consistent. SFEXCH3A.545
C----------------------------------------------------------------------- SFEXCH3A.546
C SFEXCH3A.547
IF (LTIMER) THEN SFEXCH3A.548
CALL TIMER('SFEXCH ',3) SFEXCH3A.549
ENDIF SFEXCH3A.550
SFEXCH3A.551
ERROR=0 SFEXCH3A.552
*IF DEF,MPP AJC1F405.112
IF ( AJC1F405.113
*ELSEIF DEF,SCMA AJC1F405.114
IF ( U_ROWS .NE. (P_ROWS) .OR. AJC1F405.115
*ELSE AJC1F405.116
IF ( U_ROWS .NE. (P_ROWS+1) .OR. AJC1F405.117
*ENDIF AJC1F405.118
+ U_POINTS .NE. (U_ROWS*ROW_LENGTH) .OR. SFEXCH3A.571
+ P_POINTS .NE. (P_ROWS*ROW_LENGTH) .OR. SFEXCH3A.572
+ LAND_PTS .GT. P_POINTS ) THEN SFEXCH3A.573
ERROR=1 SFEXCH3A.574
GOTO6 SFEXCH3A.575
ENDIF SFEXCH3A.576
C SFEXCH3A.577
C----------------------------------------------------------------------- SFEXCH3A.578
CL 1. Construct SICE_INDEX for compression onto sea points in SFEXCH3A.579
CL sea-ice leads calculations. SFEXCH3A.580
C----------------------------------------------------------------------- SFEXCH3A.581
C SFEXCH3A.582
DO I = 1,P_POINTS SFEXCH3A.583
ITEST(I) = .FALSE. SFEXCH3A.584
IF (ICE_FRACT(I).GT.0.0 .AND. .NOT.LAND_MASK(I)) SFEXCH3A.585
& ITEST(I) = .TRUE. SFEXCH3A.586
ENDDO SFEXCH3A.587
C SFEXCH3A.588
C Routine whenimd is functionally equivalent to WHENILE, so ITEST is SFEXCH3A.589
C 1 for "False", 0 for "True". SFEXCH3A.590
C SFEXCH3A.591
C GSS1F402.72
NSICE = 0 GSS1F402.73
DO I=1,P_POINTS GSS1F402.74
IF(ITEST(I))THEN GSS1F402.75
NSICE = NSICE + 1 GSS1F402.76
SICE_INDEX(NSICE) = I GSS1F402.77
END IF GSS1F402.78
END DO GSS1F402.79
C SFEXCH3A.594
C----------------------------------------------------------------------- SFEXCH3A.595
CL 2. Calculate QSAT values required later and components of ocean SFEXCH3A.596
CL current. SFEXCH3A.597
C Done here to avoid loop splitting. SFEXCH3A.598
C QSTAR 'borrowed' to store P at level 1 (just this once). SFEXCH3A.599
C PSIS 'borrowed' to store leads and non sea-ice surface temp. SFEXCH3A.600
C----------------------------------------------------------------------- SFEXCH3A.601
C SFEXCH3A.602
C----------------------------------------------------------------------- SFEXCH3A.604
CL 2.1 IF (GATHER) THEN SFEXCH3A.605
CL Calculate temperatures and pressures for QSAT calculations. SFEXCH3A.606
CL Calculate QSAT values. For sea-ice points, separate values SFEXCH3A.607
CL are required for the leads (QSTAR_LEAD) and sea-ice (QSTAR) SFEXCH3A.608
CL fractions respectively. QSTAR_LEAD = missing data, elsewhere. SFEXCH3A.609
CL Use RS to store compressed PSTAR for this section only. SFEXCH3A.610
CL NB Unlike QSTAR, TSTAR values at sea-ice points are gridsq. SFEXCH3A.611
CL means and so include the leads contribution. SFEXCH3A.612
CL ELSE SFEXCH3A.613
CL As above with QSTAR_LEAD done on full field. SFEXCH3A.614
CL ENDIF SFEXCH3A.615
C----------------------------------------------------------------------- SFEXCH3A.616
IF (GATHER) THEN SFEXCH3A.617
DO I = 1,P_POINTS SFEXCH3A.618
IF (L_BL_LSPICE) THEN ADM3F404.240
TL_1(I) = T_1(I) - LCRCP*QCL_1(I) ! P243.9 ADM3F404.241
ELSE ADM3F404.242
TL_1(I) = T_1(I) - LCRCP*QCL_1(I) - LSRCP*QCF_1(I) !P243.9 ADM3F404.243
ENDIF ADM3F404.244
TSTAR_NL(I) = TSTAR(I) SFEXCH3A.620
QSTAR_LEAD(I) = 1.0E30 ! Missing data indicator SFEXCH3A.621
QSTAR(I) = AK_1 + BK_1*PSTAR(I) SFEXCH3A.622
ENDDO SFEXCH3A.623
IF (NSICE.GT.0) THEN SFEXCH3A.624
CDIR$ IVDEP SFEXCH3A.625
! Fujitsu vectorization directive GRB0F405.465
!OCL NOVREC GRB0F405.466
DO SI = 1,NSICE SFEXCH3A.626
I = SICE_INDEX(SI) SFEXCH3A.627
TSTAR_NL(I) = (TSTAR(I)-(1.0-ICE_FRACT(I)) *TFS) SFEXCH3A.628
+ / ICE_FRACT(I) ! P2430.1 SFEXCH3A.629
PSIS(SI) = TFS SFEXCH3A.630
RS(SI) = PSTAR(I) SFEXCH3A.631
ENDDO SFEXCH3A.632
ENDIF SFEXCH3A.633
IF (L_BL_LSPICE) THEN ADM3F404.245
CALL QSAT_WAT(QSL,TL_1,QSTAR,P_POINTS) ADM3F404.246
ELSE ADM3F404.247
CALL QSAT(QSL,TL_1,QSTAR,P_POINTS) ADM3F404.248
ENDIF ADM3F404.249
SFEXCH3A.635
CALL QSAT(QSTAR,TSTAR_NL,PSTAR,P_POINTS) SFEXCH3A.636
C ...values at sea-ice points contain ice contribution only SFEXCH3A.637
IF (NSICE.GT.0) CALL QSAT(QSTAR_LEAD,PSIS,RS,NSICE) SFEXCH3A.638
C ...values at sea-ice points only SFEXCH3A.639
ELSE SFEXCH3A.640
C----------------------------------------------------------------------- SFEXCH3A.642
CL 2.1 Single Column Model selected. SFEXCH3A.643
CL Calculate temperatures and pressures for QSAT calculations. SFEXCH3A.644
CL If there is sea-ice, separate values of surface saturated SFEXCH3A.645
CL specific humidity are required for the leads (QSTAR_LEAD) SFEXCH3A.646
CL and sea-ice (QSTAR) fractions respectively. SFEXCH3A.647
CL NB Unlike QSTAR, TSTAR values at sea-ice points are gridsq. SFEXCH3A.648
CL means and so include the leads contribution. SFEXCH3A.649
C----------------------------------------------------------------------- SFEXCH3A.650
DO I = 1,P_POINTS SFEXCH3A.652
IF (L_BL_LSPICE) THEN ADM3F404.250
TL_1(I) = T_1(I) - LCRCP*QCL_1(I) ! P243.9 ADM3F404.251
ELSE ADM3F404.252
TL_1(I) = T_1(I) - LCRCP*QCL_1(I) - LSRCP*QCF_1(I) !P243.9 ADM3F404.253
ENDIF ADM3F404.254
TSTAR_NL(I) = TSTAR(I) SFEXCH3A.654
C Set to missing data at non sea-ice points after QSAT. SFEXCH3A.655
PSIS(I) = TSTAR(I) SFEXCH3A.656
QSTAR(I) = AK_1 + BK_1*PSTAR(I) SFEXCH3A.657
IF ( ICE_FRACT(I).GT.0.0 .AND. .NOT.LAND_MASK(I) ) THEN SFEXCH3A.658
TSTAR_NL(I) = (TSTAR(I)-(1.0-ICE_FRACT(I)) *TFS) SFEXCH3A.659
+ / ICE_FRACT(I) ! P2430.1 SFEXCH3A.660
PSIS(I) = TFS SFEXCH3A.661
ENDIF SFEXCH3A.662
ENDDO SFEXCH3A.663
IF (L_BL_LSPICE) THEN ADM3F404.255
CALL QSAT_WAT(QSL,TL_1,QSTAR,P_POINTS) ADM3F404.256
ELSE ADM3F404.257
CALL QSAT(QSL,TL_1,QSTAR,P_POINTS) ADM3F404.258
ENDIF ADM3F404.259
SFEXCH3A.665
CALL QSAT(QSTAR,TSTAR_NL,PSTAR,P_POINTS) SFEXCH3A.666
C ...values at sea-ice points contain ice contribution only SFEXCH3A.667
IF (NSICE.GT.0) CALL QSAT(QSTAR_LEAD,PSIS,PSTAR,P_POINTS) SFEXCH3A.668
C ...values at sea-ice points contain leads contribution only SFEXCH3A.669
DO I=1,P_POINTS SFEXCH3A.670
IF ( .NOT.(ICE_FRACT(I).GT.0.0 .AND. .NOT.LAND_MASK(I)) ) SFEXCH3A.671
+ QSTAR_LEAD(I) = 1.0E30 SFEXCH3A.672
ENDDO SFEXCH3A.673
ENDIF ! End of IF (GATHER) THEN... ELSE. GSS1F402.80
C----------------------------------------------------------------------- SFEXCH3A.677
CL 2.2 Set components of ocean surface current. SFEXCH3A.678
C----------------------------------------------------------------------- SFEXCH3A.679
*IF DEF,SCMA AJC1F405.119
DO I = 1, U_POINTS AJC1F405.120
U_0_P(I) = U_0(I) AJC1F405.121
V_0_P(I) = V_0(I) AJC1F405.122
ENDDO AJC1F405.123
*ELSE SFEXCH3A.683
CALL UV_TO_P(U_0,U_0_P,U_POINTS,P_POINTS,ROW_LENGTH,U_ROWS) SFEXCH3A.684
CALL UV_TO_P(V_0,V_0_P,U_POINTS,P_POINTS,ROW_LENGTH,U_ROWS) SFEXCH3A.685
*ENDIF SFEXCH3A.686
C SFEXCH3A.687
C----------------------------------------------------------------------- SFEXCH3A.688
CL 3. Loop round gridpoints to be processed, performing calculations SFEXCH3A.689
CL BEFORE call to FCDCH which necessitates splitting of loop. SFEXCH3A.690
C SFEXCH3A.691
C Note that there are a number of redundant operations in this SFEXCH3A.692
C loop, mainly setting "default" values before entering an IF- SFEXCH3A.693
C block. These are necessary to prevent the flow control SFEXCH3A.694
C becoming too complicated for the compiler to vectorize, and SFEXCH3A.695
C should therefore be left in. SFEXCH3A.696
C SFEXCH3A.697
C Note too that this section has been split into several small SFEXCH3A.698
C loops instead of being kept as one larger loop. This is also SFEXCH3A.699
C done to enable the compiler to vectorize. SFEXCH3A.700
C----------------------------------------------------------------------- SFEXCH3A.701
C SFEXCH3A.702
C----------------------------------------------------------------------- SFEXCH3A.703
CL 3.1 Fix roughness lengths for the various surface types. SFEXCH3A.704
C----------------------------------------------------------------------- SFEXCH3A.705
DO I = 1,P_POINTS SFEXCH3A.706
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SFEXCH3A.707
CL 3.1.1 Liquid sea. Overwrite sea-ice and land in 3.1.2, 3.1.3. SFEXCH3A.708
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SFEXCH3A.709
Z0M(I) = Z0MSEA(I) ! P243.B5 SFEXCH3A.710
Z0H(I) = Z0HSEA ! " SFEXCH3A.711
Z0M_EFF(I) = Z0MSEA(I) SFEXCH3A.712
Z0F(I) = Z0FSEA ! " SFEXCH3A.714
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SFEXCH3A.715
CL 3.1.2 Sea ice: Z0MIZ set on all points for input to FCDCH routine SFEXCH3A.716
CL in CD_MIZ,CH_MIZ calculations. Similarily Z0HSEA,Z0FSEA for SFEXCH3A.717
CL CD_LEAD,CH_LEAD calculations. Z0SICE for CD,CH over sea-ice. SFEXCH3A.718
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SFEXCH3A.719
CD_LEAD(I) = Z0MIZ SFEXCH3A.720
Z0HS(I) = Z0HSEA SFEXCH3A.721
Z0FS(I) = Z0FSEA SFEXCH3A.722
IF (ICE_FRACT(I).GT.0.0 .AND. .NOT.LAND_MASK(I)) THEN SFEXCH3A.723
Z0M(I) = Z0SICE ! P243.B4 SFEXCH3A.724
Z0H(I) = Z0SICE ! " SFEXCH3A.725
Z0M_EFF(I) = Z0SICE SFEXCH3A.726
Z0F(I) = Z0SICE ! " SFEXCH3A.728
ENDIF SFEXCH3A.729
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SFEXCH3A.730
CL 3.1.2a Specify blending height for all points. Set to minimum value SFEXCH3A.731
CL so that LAMBDA_EFF = LAMBDA over the sea in KMKH. SFEXCH3A.732
CL This avoids having to pass land-sea mask into KMKH. SFEXCH3A.733
CL Also set the wind profile factor to the default value of 1.0 ARN1F400.20
CL for non-land points. ARN1F400.21
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SFEXCH3A.734
SFEXCH3A.735
H_BLEND(I) = H_BLEND_MIN SFEXCH3A.736
WIND_PROFILE_FACTOR(I) = 1.0 ARN1F400.22
C ARN1F400.23
ENDDO SFEXCH3A.738
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SFEXCH3A.740
CL 3.1.3 Land. Reduce roughness if there is any snow lying. SFEXCH3A.741
CL Eqns P243.B1, B2. SFEXCH3A.742
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SFEXCH3A.743
CDIR$ IVDEP SFEXCH3A.748
! Fujitsu vectorization directive GRB0F405.467
!OCL NOVREC GRB0F405.468
DO L = 1,LAND_PTS SFEXCH3A.749
I = LAND_INDEX(L) - (P1-1) SFEXCH3A.750
IF (LYING_SNOW(I) .LT. 5.0E3) THEN ! Only reduce non-orographic SFEXCH3A.752
C ! roughness for land points SFEXCH3A.753
C ! without permanent snow. SFEXCH3A.754
C SFEXCH3A.755
Z0 = Z0V(I) - 4.0E-4 * LYING_SNOW(I) SFEXCH3A.756
ZETA1 = MIN( 5.0E-4 , Z0V(I) ) SFEXCH3A.757
Z0M(I) = MAX( ZETA1 , Z0 ) SFEXCH3A.758
Z0H(I) = MIN( Z0V(I) , Z0M(I) ) ARN1F400.24
Z0F(I) = Z0H(I) SFEXCH3A.760
ELSE ! for permanent land-ice Z0V is appropriate SFEXCH3A.761
Z0M(I) = Z0V(I) ! P243.B1, based on P243.B2 (2nd case) SFEXCH3A.762
Z0H(I) = Z0V(I) ! " , " " " ( " " ) SFEXCH3A.763
Z0F(I) = Z0V(I) ! " , " " " ( " " ) SFEXCH3A.764
ENDIF SFEXCH3A.765
C SFEXCH3A.766
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SFEXCH3A.767
CL 3.1.4 Orographic roughness. Calculate Z0M_EFF in neutral case. ARN1F400.25
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SFEXCH3A.770
IF (L_Z0_OROG) THEN SFEXCH3A.771
C ARN1F400.26
C ! Set blending height, effective roughness length and ARN1F400.27
C ! wind profile factor at land points. ARN1F400.28
C SFEXCH3A.775
SFEXCH3A.776
IF (SIL_OROG(L) .EQ. RMDI) THEN SFEXCH3A.779
ZETA1 = 0.0 ARN1F400.29
ELSE SFEXCH3A.781
ZETA1 = 0.5 * OROG_DRAG_PARAM * SIL_OROG(L) SFEXCH3A.782
ENDIF SFEXCH3A.783
ZETA2 = LOG ( 1.0 + Z1(I)/Z0M(I) ) ARN1F400.30
ZETA3 = 1.0 / SQRT ( ZETA1/(VKMAN*VKMAN) + 1.0/(ZETA2*ZETA2) ) ARN1F400.31
ZETA2 = 1.0 / EXP(ZETA3) ARN1F400.32
H_BLEND(I) = MAX ( Z1(I) / (1.0 - ZETA2) , ARN1F400.33
& HO2R2_OROG(L) * SQRT(2.0) ) ARN1F400.34
H_BLEND(I) = MIN ( H_BLEND_MAX, H_BLEND(I) ) ARN1F400.35
SFEXCH3A.784
Z0M_EFF(I) = H_BLEND(I) / EXP ( VKMAN / SQRT ( ZETA1 + SFEXCH3A.785
& (VKMAN / LOG ( H_BLEND(I) / Z0M(I) ) ) * ARN1F400.36
& (VKMAN / LOG ( H_BLEND(I) / Z0M(I) ) ) ) ) ARN1F400.37
SFEXCH3A.787
WIND_PROFILE_FACTOR(I) = LOG ( H_BLEND(I) / Z0M_EFF(I) ) / ARN1F400.38
& LOG ( H_BLEND(I) / Z0M(I) ) ARN1F400.39
SFEXCH3A.788
ELSE ! Orographic roughness not represented so ARN1F400.40
C ! leave blending height and wind profile factor at their ARN1F400.41
C ! default values and set effective roughness length to its ARN1F400.42
C ! value based on land type. ARN1F400.43
C ARN1F400.44
Z0M_EFF(I) = Z0M(I) ARN1F400.45
ENDIF ARN1F400.46
SFEXCH3A.798
ENDDO SFEXCH3A.803
C----------------------------------------------------------------------- SFEXCH3A.805
CL 3.2 Calculate buoyancy parameters for the lowest model level. SFEXCH3A.806
C----------------------------------------------------------------------- SFEXCH3A.807
DO I=1,P_POINTS SFEXCH3A.808
IF (L_BL_LSPICE) THEN ADM3F404.260
QW_1(I) = Q_1(I) + QCL_1(I) ! P243.10 ADM3F404.261
ELSE ADM3F404.262
QW_1(I) = Q_1(I) + QCL_1(I) + QCF_1(I) ! P243.10 ADM3F404.263
ENDIF ADM3F404.264
BETAT = 1.0 / T_1(I) ! P243.19 (1st eqn) SFEXCH3A.810
BETAQ = C_VIRTUAL / SFEXCH3A.811
+ ( 1.0 + C_VIRTUAL*Q_1(I) - QCL_1(I) - QCF_1(I) ) SFEXCH3A.812
C ... P243.19 (2nd eqn) SFEXCH3A.813
C SFEXCH3A.814
IF (TL_1(I).GT.TM.OR.L_BL_LSPICE) THEN ADM3F404.265
ALPHAL = (EPSILON * LC * QSL(I)) / (R * TL_1(I) * TL_1(I)) SFEXCH3A.816
C ... P243.20 (Clausius-Clapeyron) SFEXCH3A.817
C SFEXCH3A.818
AL = 1.0 / (1.0 + LCRCP*ALPHAL) ! P243.21 SFEXCH3A.819
BETAC = CF_1(I) * AL * (LCRCP*BETAT - ETAR*BETAQ) SFEXCH3A.820
C ... P243.19 (3rd eqn) SFEXCH3A.821
C SFEXCH3A.822
ELSE SFEXCH3A.823
ALPHAL = (EPSILON * LS * QSL(I)) / (R * TL_1(I) * TL_1(I)) SFEXCH3A.824
C ... P243.20 (Clausius-Clapeyron) SFEXCH3A.825
C SFEXCH3A.826
AL = 1.0 / (1.0 + LSRCP*ALPHAL) ! P243.21 SFEXCH3A.827
BETAC = CF_1(I) * AL * (LSRCP*BETAT - ETAR*BETAQ) SFEXCH3A.828
C ... P243.19 (3rd eqn) SFEXCH3A.829
C SFEXCH3A.830
ENDIF SFEXCH3A.831
BT_1(I) = BETAT - ALPHAL*BETAC ! P243.18 (1st eqn) SFEXCH3A.832
BQ_1(I) = BETAQ + BETAC ! P243.18 (2nd eqn) SFEXCH3A.833
BF_1(I) = BETAQ *EPSILON*ETAR ! P243.18 (2nd eqn) ADM3F404.266
ENDDO SFEXCH3A.834
C----------------------------------------------------------------------- SFEXCH3A.835
CL 3.3 Calculate temperature (strictly, liquid/ice static energy) and SFEXCH3A.836
CL humidity jumps, and wind shear, across the surface layer. SFEXCH3A.837
CL Separate values are required for the leads and ice fractions SFEXCH3A.838
CL of sea-ice grid-squares. SFEXCH3A.839
C----------------------------------------------------------------------- SFEXCH3A.840
IF (GATHER) THEN SFEXCH3A.842
DO I=1,P_POINTS SFEXCH3A.843
DTEMP(I) = TL_1(I) - TSTAR_NL(I) SFEXCH3A.844
& + GRCP * ( Z1(I) + Z0M_EFF(I) - Z0H(I) ) ARN1F400.47
C ! P243.118 SFEXCH3A.846
DQ(I) = QW_1(I) - QSTAR(I) ! P243.119 SFEXCH3A.847
DTEMP_LEAD(I) = 1.0E30 ! Missing data indicator SFEXCH3A.848
DQ_LEAD(I) = 1.0E30 ! Missing data indicator SFEXCH3A.849
USHEAR = U_1_P(I) - U_0_P(I) SFEXCH3A.850
VSHEAR = V_1_P(I) - V_0_P(I) SFEXCH3A.851
VSHR2 = MAX (1.0E-6 , USHEAR*USHEAR + VSHEAR*VSHEAR) SFEXCH3A.852
VSHR(I) = SQRT(VSHR2) SFEXCH3A.853
C ... P243.120 (previous 4 lines of code) SFEXCH3A.854
ENDDO SFEXCH3A.855
C----------------------------------------------------------------------- SFEXCH3A.856
CL 3.3.1 Calculate leads values by looping round sea-ice points only. SFEXCH3A.857
C Avoids an if test in the above loop, so code can run faster. SFEXCH3A.858
C----------------------------------------------------------------------- SFEXCH3A.859
CDIR$ IVDEP SFEXCH3A.860
! Fujitsu vectorization directive GRB0F405.469
!OCL NOVREC GRB0F405.470
DO SI=1,NSICE SFEXCH3A.861
I = SICE_INDEX(SI) SFEXCH3A.862
DTEMP_LEAD(I) = TL_1(I)-TFS + GRCP*(Z1(I)+Z0MSEA(I)-Z0HS(I)) SFEXCH3A.863
DQ_LEAD(I) = QW_1(I) - QSTAR_LEAD(SI) SFEXCH3A.864
ENDDO SFEXCH3A.865
ELSE SFEXCH3A.866
DO I=1,P_POINTS SFEXCH3A.868
USHEAR = U_1_P(I) - U_0_P(I) SFEXCH3A.873
VSHEAR = V_1_P(I) - V_0_P(I) SFEXCH3A.874
VSHR2 = MAX (1.0E-6 , USHEAR*USHEAR + VSHEAR*VSHEAR) SFEXCH3A.876
VSHR(I) = SQRT(VSHR2) SFEXCH3A.877
C ... P243.120 (previous 4 lines of code) SFEXCH3A.878
DTEMP(I) = TL_1(I) - TSTAR_NL(I) SFEXCH3A.879
& + GRCP * ( Z1(I) + Z0M_EFF(I) - Z0H(I) ) ARN1F400.48
C ! P243.118 SFEXCH3A.881
DQ(I) = QW_1(I) - QSTAR(I) ! P243.119 SFEXCH3A.882
IF ( ICE_FRACT(I).GT.0.0 .AND. .NOT.LAND_MASK(I) ) THEN SFEXCH3A.883
DTEMP_LEAD(I) = TL_1(I)-TFS + GRCP*(Z1(I)+Z0MSEA(I)-Z0HS(I)) SFEXCH3A.884
DQ_LEAD(I) = QW_1(I) - QSTAR_LEAD(I) SFEXCH3A.885
ELSE SFEXCH3A.886
DTEMP_LEAD(I) = 1.0E30 ! Missing data indicator SFEXCH3A.887
DQ_LEAD(I) = 1.0E30 ! Missing data indicator SFEXCH3A.888
ENDIF SFEXCH3A.889
C SFEXCH3A.890
ENDDO SFEXCH3A.891
ENDIF ! End of IF (GATHER) THEN... ELSE... SFEXCH3A.893
C----------------------------------------------------------------------- SFEXCH3A.895
CL 3.4 Evaporation over land surfaces without snow is limited by SFEXCH3A.896
CL soil moisture availability and stomatal resistance. SFEXCH3A.897
C Set FRACA (= fA in the documentation) according to P243.68, SFEXCH3A.898
C PSIS according to P243.65, and RESFS (= fS) according to P243.75 SFEXCH3A.899
C and P243.61, using neutral-stability value of CH (as explained SFEXCH3A.900
C in section (v) of the P243 documentation). SFEXCH3A.901
C----------------------------------------------------------------------- SFEXCH3A.902
DO I=1,P_POINTS SFEXCH3A.903
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SFEXCH3A.904
CL 3.4.1 Set parameters (workspace) to values relevant for non-land SFEXCH3A.905
CL points. Only aerodynamic resistance applies. SFEXCH3A.906
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SFEXCH3A.907
FRACA(I) = 1.0 SFEXCH3A.908
PSIS(I) = 0.0 SFEXCH3A.909
RESFT(I) = 1.0 SFEXCH3A.910
RS(I) = 0.0 SFEXCH3A.911
RESFS(I) = 0.0 SFEXCH3A.912
ENDDO SFEXCH3A.914
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SFEXCH3A.916
CL 3.4.2 Over-write workspace for other points. SFEXCH3A.917
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SFEXCH3A.918
CDIR$ IVDEP SFEXCH3A.923
! Fujitsu vectorization directive GRB0F405.471
!OCL NOVREC GRB0F405.472
DO L=1,LAND_PTS SFEXCH3A.924
I = LAND_INDEX(L) - (P1-1) SFEXCH3A.925
C SFEXCH3A.927
C Calculate the soil moisture availability factor, PSIS. SFEXCH3A.928
C SFEXCH3A.929
SMCRIT = RHO_WATER * ROOTD(L) * (SMVCCL(L)-SMVCWT(L)) SFEXCH3A.930
C ... P243.66 SFEXCH3A.931
C SFEXCH3A.932
PSIS(I) = 0.0 SFEXCH3A.933
IF (SMC(L).GE.SMCRIT .AND. SMCRIT.GT.0.0) SFEXCH3A.934
+ PSIS(I) = 1.0 SFEXCH3A.935
IF (SMC(L).LT.SMCRIT .AND. SMC(L).GT.0.0) SFEXCH3A.936
+ PSIS(I) = SMC(L)/SMCRIT SFEXCH3A.937
C SFEXCH3A.938
C For snow-covered land or land points with negative moisture flux SFEXCH3A.939
C set the fraction of the flux with only aerodynamic resistance to 1 SFEXCH3A.940
C (no surface/stomatal resistance to evap from snow or condensation). SFEXCH3A.941
C SFEXCH3A.942
FRACA(I) = 1.0 SFEXCH3A.943
C SFEXCH3A.944
C When there is positive moisture flux from snow-free land, calculate SFEXCH3A.945
C the fraction of the flux from the "canopy". SFEXCH3A.946
C SFEXCH3A.947
IF (DQ(I).LT.0.0 .AND. LYING_SNOW(I).LE.0.0) FRACA(I) = 0.0 SFEXCH3A.948
IF (DQ(I).LT.0.0.AND.LYING_SNOW(I).LE.0.0.AND.CATCH(L).GT.0.0) SFEXCH3A.949
+ FRACA(I) = CANOPY(L)/CATCH(L) SFEXCH3A.950
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SFEXCH3A.951
CL 3.4.2.1 Calculate neutral stability value of CH (CHN), as an SFEXCH3A.952
CL approximation to CH. SFEXCH3A.953
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SFEXCH3A.954
ZETAM = LOG ( (Z1(I) + Z0M_EFF(I))/Z0M_EFF(I) ) SFEXCH3A.955
ZETAH = LOG ( (Z1(I) + Z0M_EFF(I))/Z0H(I) ) ARN1F400.49
CHN = (VKMAN/ZETAH) * (VKMAN/ZETAM) * WIND_PROFILE_FACTOR(I) ARN1F400.50
C ... P243.41 (previous 3 lines of code) SFEXCH3A.958
C SFEXCH3A.959
RS(I) = RESIST(L) SFEXCH3A.960
RESFS(I) = PSIS(I) / ( 1.0 + CHN*VSHR(I)*RS(I) ) SFEXCH3A.961
RESFT(I) = FRACA(I) + (1.0 - FRACA(I)) * RESFS(I) SFEXCH3A.962
ENDDO ! Evaporation over land points only, section 3.4.2 SFEXCH3A.967
C----------------------------------------------------------------------- SFEXCH3A.969
CL 3.5 Calculate bulk Richardson numbers for the surface layer. SFEXCH3A.970
CL At sea-ice points RIB contains value for ice only (not leads). SFEXCH3A.971
CL Initialise RIB_LEAD to RIB so that it contains sensible SFEXCH3A.972
CL values at non sea ice points for the FCDCH calculation below. SFEXCH3A.973
C----------------------------------------------------------------------- SFEXCH3A.974
IF (GATHER) THEN SFEXCH3A.976
DO I=1,P_POINTS SFEXCH3A.977
RIB(I) = G * Z1(I) * SFEXCH3A.978
+ ( BT_1(I)*DTEMP(I) + BQ_1(I)*RESFT(I)*DQ(I) ) / SFEXCH3A.979
+ ( VSHR(I)*VSHR(I) ) SFEXCH3A.980
RIB_LEAD(I) = RIB(I) SFEXCH3A.981
ENDDO SFEXCH3A.982
C----------------------------------------------------------------------- SFEXCH3A.983
CL 3.5.1 Calculate bulk Richardson no. for leads at sea-ice points SFEXCH3A.984
CL only. SFEXCH3A.985
C----------------------------------------------------------------------- SFEXCH3A.986
CDIR$ IVDEP SFEXCH3A.987
! Fujitsu vectorization directive GRB0F405.473
!OCL NOVREC GRB0F405.474
DO SI = 1,NSICE SFEXCH3A.988
I = SICE_INDEX(SI) SFEXCH3A.989
RIB_LEAD(I) = G * Z1(I) * SFEXCH3A.990
+ ( BT_1(I) * DTEMP_LEAD(I) + SFEXCH3A.991
+ BQ_1(I) * RESFT(I) * DQ_LEAD(I) ) / SFEXCH3A.992
+ ( VSHR(I) * VSHR(I) ) SFEXCH3A.993
C ... P2430.2, for sea-ice leads. SFEXCH3A.994
ENDDO SFEXCH3A.995
ELSE SFEXCH3A.996
DO I=1,P_POINTS SFEXCH3A.998
RIB(I) = G * Z1(I) * SFEXCH3A.999
+ ( BT_1(I)*DTEMP(I) + BQ_1(I)*RESFT(I)*DQ(I) ) / SFEXCH3A.1000
+ ( VSHR(I)*VSHR(I) ) SFEXCH3A.1001
C ... P243.43 (G times middle line is surface SFEXCH3A.1002
C layer buoyancy difference, P243.25) SFEXCH3A.1003
C SFEXCH3A.1004
RIB_LEAD(I) = RIB(I) SFEXCH3A.1005
IF ( ICE_FRACT(I).GT.0.0 .AND. .NOT.LAND_MASK(I) ) THEN SFEXCH3A.1006
RIB_LEAD(I) = G * Z1(I) * SFEXCH3A.1007
+ ( BT_1(I) * DTEMP_LEAD(I) + SFEXCH3A.1008
+ BQ_1(I) * RESFT(I) * DQ_LEAD(I) ) / SFEXCH3A.1009
+ ( VSHR(I) * VSHR(I) ) SFEXCH3A.1010
C ... P2430.2, for sea-ice leads. SFEXCH3A.1011
ENDIF SFEXCH3A.1012
ENDDO SFEXCH3A.1013
ENDIF ! End of IF (GATHER) THEN... ELSE... SFEXCH3A.1015
C SFEXCH3A.1017
C----------------------------------------------------------------------- SFEXCH3A.1018
CL 3.6 Calculate stability corrected effective roughness length. ARN1F400.51
CL Simple linear interpolation when RIB between 0 and RIB_CRIT (>0) for ARN1F400.52
CL form drag term. ARN1F400.53
C----------------------------------------------------------------------- SFEXCH3A.1022
SFEXCH3A.1023
IF (L_Z0_OROG) THEN SFEXCH3A.1024
CDIR$ IVDEP SFEXCH3A.1029
! Fujitsu vectorization directive GRB0F405.475
!OCL NOVREC GRB0F405.476
DO L = 1,LAND_PTS SFEXCH3A.1030
I = LAND_INDEX(L) - (P1-1) SFEXCH3A.1031
SFEXCH3A.1033
IF (RIB(I) .GE. RI_CRIT) THEN ARN1F400.54
Z0M_EFF(I) = Z0M(I) SFEXCH3A.1035
ELSEIF (RIB(I) .GT. 0.0 ) THEN ARN1F400.55
IF (SIL_OROG(L) .EQ. RMDI) THEN ARN1F400.56
ZETA1 = 0.0 ARN1F400.57
ELSE ARN1F400.58
ZETA1 = 0.5 * OROG_DRAG_PARAM * SIL_OROG(L) * ARN1F400.59
& ( 1.0 - RIB(I) / RI_CRIT ) ARN1F400.60
ENDIF ARN1F400.61
Z0M_EFF(I) = H_BLEND(I) / EXP ( VKMAN / SQRT ( ZETA1 + ARN1F400.62
& (VKMAN / LOG ( H_BLEND(I) / Z0M(I) ) ) * ARN1F400.63
& (VKMAN / LOG ( H_BLEND(I) / Z0M(I) ) ) ) ) ARN1F400.64
ARN1F400.65
ENDIF SFEXCH3A.1039
SFEXCH3A.1040
WIND_PROFILE_FACTOR(I) = LOG ( H_BLEND(I) / Z0M_EFF(I) ) / ARN1F400.66
& LOG ( H_BLEND(I) / Z0M(I) ) ARN1F400.67
SFEXCH3A.1041
ARN1F400.68
ENDDO SFEXCH3A.1050
ENDIF SFEXCH3A.1052
SFEXCH3A.1053
C SFEXCH3A.1054
C----------------------------------------------------------------------- SFEXCH3A.1055
CL 3.7 Calculate CD, CH via routine FCDCH. SFEXCH3A.1056
CL Calculate CD_MIZ,CH_MIZ,CD_LEAD,CH_LEAD on full field then set SFEXCH3A.1057
CL non sea-ice points to missing data (contain nonsense after FCDCH) SFEXCH3A.1058
C Unlike the QSAT calculations above, arrays are not compressed to SFEXCH3A.1060
C sea-ice points for FCDCH. This is because it would require extra SFEXCH3A.1061
C work space and initial tests showed that with with the extra SFEXCH3A.1062
C compression calculations required no time was saved. SFEXCH3A.1063
C NB CD_LEAD stores Z0MIZ for calculation of CD_MIZ,CH_MIZ. SFEXCH3A.1065
C----------------------------------------------------------------------- SFEXCH3A.1066
C SFEXCH3A.1067
CALL FCDCH(RIB,CD_LEAD,CD_LEAD,CD_LEAD,Z1,WIND_PROFILE_FACTOR, ARN1F400.69
& P_POINTS,CD_MIZ,CH_MIZ,CD_STD,LTIMER) ARN1F400.70
C ! Marginal Ice Zone.P2430.9 ARN1F400.71
CALL FCDCH(RIB_LEAD,Z0MSEA,Z0HS,Z0FS,Z1,WIND_PROFILE_FACTOR, ARN1F400.72
& P_POINTS,CD_LEAD,CH_LEAD,CD_STD,LTIMER) ARN1F400.73
C ! Sea-ice leads.P2430.8 ARN1F400.74
CALL FCDCH(RIB,Z0M_EFF,Z0H,Z0F,Z1,WIND_PROFILE_FACTOR, ARN1F400.75
& P_POINTS,CD,CH,CD_STD,LTIMER) ARN1F400.76
DO I=1,P_POINTS SFEXCH3A.1073
C IF ( an ordinary sea points (no sea-ice) or a land point) SFEXCH3A.1074
IF (.NOT.(ICE_FRACT(I).GT.0.0 .AND. .NOT.LAND_MASK(I)) ) THEN SFEXCH3A.1075
CD_MIZ(I) = 1.0E30 SFEXCH3A.1076
CH_MIZ(I) = 1.0E30 SFEXCH3A.1077
CD_LEAD(I) = 1.0E30 SFEXCH3A.1078
CH_LEAD(I) = 1.0E30 SFEXCH3A.1079
RIB_LEAD(I) = 1.0E30 SFEXCH3A.1080
ENDIF SFEXCH3A.1081
ENDDO SFEXCH3A.1082
C SFEXCH3A.1083
C----------------------------------------------------------------------- SFEXCH3A.1084
CL 4. Loop round gridpoints to be processed, performing calculations SFEXCH3A.1085
CL AFTER call to FCDCH which necessitates splitting of loop. SFEXCH3A.1086
C----------------------------------------------------------------------- SFEXCH3A.1087
CL 4.1 Recalculate RESFS using "true" CH. SFEXCH3A.1088
C----------------------------------------------------------------------- SFEXCH3A.1089
CDIR$ IVDEP SFEXCH3A.1095
! Fujitsu vectorization directive GRB0F405.477
!OCL NOVREC GRB0F405.478
DO L = 1,LAND_PTS SFEXCH3A.1096
I = LAND_INDEX(L) - (P1-1) SFEXCH3A.1097
RESFS(I) = PSIS(I) / ( 1.0 + CH(I)*VSHR(I)*RS(I) ) ! P243.75 SFEXCH3A.1099
EPDT = -PSTAR(I)/(R*TSTAR(I))*CH(I)*VSHR(I)*DQ(I)*TIMESTEP SFEXCH3A.1100
IF (EPDT .GT. 0.0 .AND. LYING_SNOW(I) .LE. 0.0 SFEXCH3A.1101
& .AND. CATCH(L) .GT. 0.0) THEN SFEXCH3A.1102
FRACA(I) = CANOPY(L)/(CATCH(L)+EPDT) SFEXCH3A.1103
RESFT(I) = FRACA(I) + (1.0 - FRACA(I)) * RESFS(I) SFEXCH3A.1104
ENDIF ! Positive evaporation from snow-free land SFEXCH3A.1105
ENDDO ! Loop over land-points SFEXCH3A.1110
C----------------------------------------------------------------------- SFEXCH3A.1112
CL 4.D Call SFL_INT to calculate CDR10M, CHR1P5M and CER1P5M - SFEXCH3A.1113
CL interpolation coefficients used in SF_EVAP and IMPL_CAL to SFEXCH3A.1114
CL calculate screen temperature, specific humidity and 10m winds. SFEXCH3A.1115
C----------------------------------------------------------------------- SFEXCH3A.1118
C SFEXCH3A.1119
IF (SU10 .OR. SV10 .OR. SQ1P5 .OR. ST1P5) THEN SFEXCH3A.1120
SFEXCH3A.1129
CALL SFL_INT ( SFEXCH3A.1130
& P_POINTS,U_POINTS,RIB,Z1,Z0M,Z0M_EFF,Z0H,Z0F,CD_STD,CD,CH, ARN1F400.77
& RESFT,WIND_PROFILE_FACTOR, ARN1F400.78
& CDR10M,CHR1P5M,CER1P5M, ARN1F400.79
+ SU10,SV10,ST1P5,SQ1P5,LTIMER SFEXCH3A.1133
+ ) SFEXCH3A.1134
ENDIF SFEXCH3A.1135
C----------------------------------------------------------------------- SFEXCH3A.1136
CL 4.2 Now that diagnostic calculations are over, update CD and CH SFEXCH3A.1137
CL to their correct values (i.e. gridsquare means). SFEXCH3A.1138
C----------------------------------------------------------------------- SFEXCH3A.1139
DO I = 1,P_POINTS SFEXCH3A.1140
IF ( ICE_FRACT(I).GT.0.0 .AND. .NOT.LAND_MASK(I) ) THEN SFEXCH3A.1141
IF ( ICE_FRACT(I).LT. 0.7 ) THEN SFEXCH3A.1142
CD(I) = ( ICE_FRACT(I)*CD_MIZ(I) + SFEXCH3A.1143
+ (0.7-ICE_FRACT(I))*CD_LEAD(I) ) / 0.7 ! P2430.5 SFEXCH3A.1144
CH(I) = ( ICE_FRACT(I)*CH_MIZ(I) + SFEXCH3A.1145
+ (0.7-ICE_FRACT(I))*CH_LEAD(I) ) / 0.7 ! P2430.4 SFEXCH3A.1146
ELSE SFEXCH3A.1147
CD(I) = ( (1.0-ICE_FRACT(I))*CD_MIZ(I) + SFEXCH3A.1148
+ (ICE_FRACT(I)-0.7)*CD(I) ) / 0.3 ! P2430.7 SFEXCH3A.1149
CH(I) = ( (1.0-ICE_FRACT(I))*CH_MIZ(I) + SFEXCH3A.1150
+ (ICE_FRACT(I)-0.7)*CH(I) ) / 0.3 ! P2430.7 SFEXCH3A.1151
ENDIF SFEXCH3A.1152
ENDIF SFEXCH3A.1153
C----------------------------------------------------------------------- SFEXCH3A.1154
CL 4.3 Calculate the surface exchange coefficients RHOK(*). SFEXCH3A.1155
C----------------------------------------------------------------------- SFEXCH3A.1156
RHOSTAR = PSTAR(I) / ( R*TSTAR(I) ) SFEXCH3A.1157
C ... surface air density from ideal gas equation SFEXCH3A.1158
C SFEXCH3A.1159
RHOKM_1(I) = RHOSTAR * CD(I) * VSHR(I) ! P243.124 SFEXCH3A.1160
RHOKH_1(I) = RHOSTAR * CH(I) * VSHR(I) ! P243.125 SFEXCH3A.1161
*IF DEF,SCMA AJC0F405.94
If (OBS) then AJC0F405.95
RHOKEA(I) = FACTOR_RHOKH(I) * FRACA(I) AJC0F405.96
else AJC0F405.97
*ENDIF AJC0F405.98
RHOKEA(I) = RHOKH_1(I) * FRACA(I) ! P243.126 SFEXCH3A.1162
*IF DEF,SCMA AJC0F405.99
endif AJC0F405.100
*ENDIF AJC0F405.101
RHOKESL(I) = RHOKH_1(I) * RESFS(I) ! P243.127 SFEXCH3A.1163
RHOKES(I) = RHOKESL(I) * (1.0 - FRACA(I)) ! P243.128 SFEXCH3A.1164
RHOKE(I) = RHOKEA(I) + RHOKES(I) ! P243.129 SFEXCH3A.1165
C ASJ1F400.29
C RHOSTAR * CD * VSHR stored for diagnostic output before ASJ1F400.30
C horizontal interpolation. ASJ1F400.31
C ASJ1F400.32
RHO_CD_MODV1(I) = RHOKM_1(I) ASJ1F400.33
C ASJ1F400.34
C----------------------------------------------------------------------- SFEXCH3A.1166
CL 4.4 Calculate the "explicit" fluxes of sensible heat and moisture. SFEXCH3A.1167
C----------------------------------------------------------------------- SFEXCH3A.1168
C SFEXCH3A.1169
C Firstly for land points, sea points with no sea-ice (ICE_FRACT=0.0) SFEXCH3A.1170
C and the sea-ice parts of the sea-ice points. SFEXCH3A.1171
FTL_1(I) = -RHOKH_1(I) * DTEMP(I) ! P243.134 SFEXCH3A.1172
FQW_1(I) = -RHOKE(I) * DQ(I) ! P243.135 SFEXCH3A.1173
EA(I) = -RHOKEA(I) * DQ(I) ! P243.136 SFEXCH3A.1174
ES(I) = -RHOKES(I) * DQ(I) ! P243.137 SFEXCH3A.1175
ESL(I) = -RHOKESL(I) * DQ(I) ! P243.138 SFEXCH3A.1176
FTL_LEAD(I) = 1.0E30 ! Missing data indicator SFEXCH3A.1177
FQW_LEAD(I) = 1.0E30 ! Missing data indicator SFEXCH3A.1178
C SFEXCH3A.1179
C Secondly for sea ice points. SFEXCH3A.1180
IF ( (ICE_FRACT(I).GT.0.0) .AND. .NOT.LAND_MASK(I) ) THEN SFEXCH3A.1181
FTL_LEAD(I) = -RHOKH_1(I)*DTEMP_LEAD(I) * (1.0-ICE_FRACT(I)) SFEXCH3A.1182
C ...P2430.11 SFEXCH3A.1183
FTL_1(I) = FTL_LEAD(I) + FTL_1(I)*ICE_FRACT(I) ! P2430.13/15 SFEXCH3A.1184
FQW_LEAD(I) = -RHOKE(I)*DQ_LEAD(I) * (1.0-ICE_FRACT(I)) SFEXCH3A.1185
C ...P2430.12 SFEXCH3A.1186
FQW_1(I) = FQW_LEAD(I) + FQW_1(I)*ICE_FRACT(I) ! P2430.14/16 SFEXCH3A.1187
ENDIF SFEXCH3A.1188
C----------------------------------------------------------------------- SFEXCH3A.1189
CL 4.4.1 Set indicator for unstable suface layer (buoyancy flux +ve.) ASJ3F401.15
CL if required by logical L_RMBL. ASJ3F401.16
C----------------------------------------------------------------------- SFEXCH3A.1191
IF ( L_RMBL.AND.BT_1(I)*FTL_1(I)+BQ_1(I)*FQW_1(I).GT.0.0) THEN ASJ3F401.17
NRML(I) = 1 SFEXCH3A.1193
ELSE SFEXCH3A.1194
NRML(I) = 0 SFEXCH3A.1195
ENDIF SFEXCH3A.1196
C----------------------------------------------------------------------- SFEXCH3A.1197
CL 4.5 Multiply surface exchange coefficients that are on the P-grid SFEXCH3A.1198
CL by GAMMA(1).Needed for implicit calculations in P244(IMPL_CAL). SFEXCH3A.1199
CL RHOKM_1 dealt with in section 4.1 below. SFEXCH3A.1200
C----------------------------------------------------------------------- SFEXCH3A.1201
RHOKE(I) = RHOKE(I) * GAMMA(1) SFEXCH3A.1202
RHOKEA(I) = RHOKEA(I) * GAMMA(1) SFEXCH3A.1203
RHOKES(I) = RHOKES(I) * GAMMA(1) SFEXCH3A.1204
RHOKESL(I) = RHOKESL(I) * GAMMA(1) SFEXCH3A.1205
RHOKH_1(I) = RHOKH_1(I) * GAMMA(1) SFEXCH3A.1206
C----------------------------------------------------------------------- SFEXCH3A.1207
CL 4.5.1 Calculate the standard deviations of layer 1 turbulent SFEXCH3A.1208
CL fluctuations of temperature and humidity using approximate SFEXCH3A.1209
CL formulae from first order closure. SFEXCH3A.1210
C----------------------------------------------------------------------- SFEXCH3A.1211
VS = SQRT ( RHOKM_1(I)/RHOSTAR * VSHR(I) ) SFEXCH3A.1212
VSF1_CUBED = 1.25 * ( Z1(I) + Z0M(I) ) * G * SFEXCH3A.1213
& ( BT_1(I)*FTL_1(I) + BQ_1(I)*FQW_1(I) ) / RHOSTAR SFEXCH3A.1214
C !--------------------------------------------------------------- SFEXCH3A.1215
C ! Only calculate standard deviations for unstable surface layers SFEXCH3A.1216
C !--------------------------------------------------------------- SFEXCH3A.1217
IF (VSF1_CUBED .GT. 0.0) THEN SFEXCH3A.1218
WS1 = ( VSF1_CUBED + VS * VS * VS ) ** (1.0/3.0) SFEXCH3A.1219
T1_SD(I) = MAX ( 0.0 , 1.93 * FTL_1(I) / (RHOSTAR * WS1) ) SFEXCH3A.1220
Q1_SD(I) = MAX ( 0.0 , 1.93 * FQW_1(I) / (RHOSTAR * WS1) ) SFEXCH3A.1221
ELSE SFEXCH3A.1222
T1_SD(I) = 0.0 SFEXCH3A.1223
Q1_SD(I) = 0.0 SFEXCH3A.1224
ENDIF SFEXCH3A.1225
C----------------------------------------------------------------------- SFEXCH3A.1226
CL 4.6 For sea points, calculate the wind mixing energy flux and the SFEXCH3A.1227
CL sea-surface roughness length on the P-grid, using time-level n SFEXCH3A.1228
CL quantities. Use CD_LEAD at sea-ice points. SFEXCH3A.1229
C----------------------------------------------------------------------- SFEXCH3A.1230
IF (.NOT.LAND_MASK(I)) THEN SFEXCH3A.1231
TAU = RHOKM_1(I) * VSHR(I) ! P243.130 SFEXCH3A.1232
IF (ICE_FRACT(I) .GT. 0.0) SFEXCH3A.1233
& TAU = RHOSTAR * CD_LEAD(I) * VSHR(I) * VSHR(I) SFEXCH3A.1234
IF (SFME) FME(I) = (1.0-ICE_FRACT(I)) * TAU * SQRT(TAU/RHOSEA) SFEXCH3A.1235
C ! P243.96 SFEXCH3A.1236
Z0MSEA(I) = MAX ( Z0HSEA , SFEXCH3A.1237
& (CHARNOCK/G) * (TAU / RHOSTAR) ) SFEXCH3A.1238
C ... P243.B6 (Charnock formula) SFEXCH3A.1239
C TAU/RHOSTAR is "mod VS squared", see eqn P243.131 SFEXCH3A.1240
C SFEXCH3A.1241
ENDIF ! of IF (.NOT. LAND_MASK), land-points done in next loop. SFEXCH3A.1245
ENDDO ! Loop over points for sections 4.2 - 4.6 SFEXCH3A.1246
DO L=1,LAND_PTS SFEXCH3A.1247
I = LAND_INDEX(L) - (P1-1) SFEXCH3A.1248
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SFEXCH3A.1250
CL 4.7 Set Z0MSEA to Z0V, FME to zero for land points. SFEXCH3A.1251
C (Former because UM uses same storage for Z0V SFEXCH3A.1252
C and Z0MSEA.) SFEXCH3A.1253
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SFEXCH3A.1254
Z0MSEA(I) = Z0V(I) SFEXCH3A.1255
IF (SFME) FME(I) = 0.0 SFEXCH3A.1256
ENDDO ! Loop over points for section 4.7 SFEXCH3A.1262
*IF -DEF,SCMA AJC1F405.124
C SFEXCH3A.1263
C----------------------------------------------------------------------- SFEXCH3A.1264
CL 5. Calculate "explicit" surface fluxes of momentum (on UV-grid). SFEXCH3A.1265
C----------------------------------------------------------------------- SFEXCH3A.1266
CL 5.1 Interpolate exchange coefficient to UV-grid, then mutiply SFEXCH3A.1267
CL by GAMMA(1) to be passed to subroutine IMPL_CAL (P244) which SFEXCH3A.1268
CL only uses RHOKM_1 when mulitplied by GAMMA(1). SFEXCH3A.1269
C----------------------------------------------------------------------- SFEXCH3A.1270
C SFEXCH3A.1271
C PSIS used purely as spare workspace here. SFEXCH3A.1272
C SFEXCH3A.1273
*IF DEF,MPP ASJ1F403.8
! RHOKM_1 contains incorrect data in halos. The P_TO_UV can interpolate ASJ1F403.9
! this into the real data, so first we must update east/west halos. ASJ1F403.10
CALL SWAPBOUNDS(RHOKM_1,ROW_LENGTH,U_POINTS/ROW_LENGTH,1,0,1) ASJ1F403.11
ASJ1F403.12
*ENDIF ASJ1F403.13
CALL P_TO_UV(RHOKM_1,PSIS,P_POINTS,U_POINTS,ROW_LENGTH,P_ROWS) SFEXCH3A.1274
DO I=1,U_POINTS-2*ROW_LENGTH SFEXCH3A.1275
J = I+ROW_LENGTH SFEXCH3A.1276
RHOKM_1(J) = PSIS(I) SFEXCH3A.1277
TAUX_1(J) = RHOKM_1(J) * ( U_1(J) - U_0(J) ) ! P243.132 SFEXCH3A.1278
TAUY_1(J) = RHOKM_1(J) * ( V_1(J) - V_0(J) ) ! P243.133 SFEXCH3A.1279
RHOKM_1(J) = GAMMA(1) * RHOKM_1(J) SFEXCH3A.1280
ENDDO SFEXCH3A.1281
C----------------------------------------------------------------------- SFEXCH3A.1282
CL 5.2 Set first and last rows to "missing data indicator". SFEXCH3A.1283
C----------------------------------------------------------------------- SFEXCH3A.1284
*IF DEF,MPP GPB1F403.7
IF (attop) THEN GPB1F403.8
*ENDIF GPB1F403.9
DO I=1,ROW_LENGTH GPB1F403.10
RHOKM_1(I) = 1.0E30 GPB1F403.11
TAUX_1(I) = 1.0E30 GPB1F403.12
TAUY_1(I) = 1.0E30 GPB1F403.13
ENDDO GPB1F403.14
*IF DEF,MPP GPB1F403.15
ENDIF GPB1F403.16
GPB1F403.17
IF (atbase) THEN GPB1F403.18
*ENDIF GPB1F403.19
DO I= (U_ROWS-1)*ROW_LENGTH + 1 , U_ROWS*ROW_LENGTH GPB1F403.20
RHOKM_1(I) = 1.0E30 GPB1F403.21
TAUX_1(I) = 1.0E30 GPB1F403.22
TAUY_1(I) = 1.0E30 GPB1F403.23
ENDDO GPB1F403.24
*IF DEF,MPP GPB1F403.25
ENDIF GPB1F403.26
*ENDIF GPB1F403.27
C----------------------------------------------------------------------- SFEXCH3A.1294
CL 5.D Interpolate CDR10M to UV-grid. SFEXCH3A.1295
C----------------------------------------------------------------------- SFEXCH3A.1296
*IF DEF,MPP ASJ1F403.14
! CDR10M contains incorrect data in halos. The P_TO_UV can interpolate ASJ1F403.15
! this into the real data, so first we must update east/west halos. ASJ1F403.16
CALL SWAPBOUNDS(CDR10M,ROW_LENGTH,U_POINTS/ROW_LENGTH,1,0,1) ASJ1F403.17
ASJ1F403.18
*ENDIF ASJ1F403.19
IF (SU10 .OR. SV10) THEN SFEXCH3A.1297
CALL P_TO_UV(CDR10M,PSIS,P_POINTS,U_POINTS,ROW_LENGTH,P_ROWS) SFEXCH3A.1298
DO I=1,U_POINTS-2*ROW_LENGTH SFEXCH3A.1299
J = I + ROW_LENGTH ARN1F400.80
CDR10M(J) = PSIS(I) ARN1F400.81
ENDDO SFEXCH3A.1301
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SFEXCH3A.1302
CL 5.D.1 Set first and last rows to "missing data indicator". SFEXCH3A.1303
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SFEXCH3A.1304
*IF DEF,MPP GPB1F403.28
IF (attop) THEN GPB1F403.29
*ENDIF GPB1F403.30
DO I=1,ROW_LENGTH GPB1F403.31
CDR10M(I) = 1.0E30 GPB1F403.32
ENDDO GPB1F403.33
*IF DEF,MPP GPB1F403.34
ENDIF GPB1F403.35
GPB1F403.36
IF (atbase) THEN GPB1F403.37
*ENDIF GPB1F403.38
DO I= (U_ROWS-1)*ROW_LENGTH + 1 , U_ROWS*ROW_LENGTH GPB1F403.39
CDR10M(I) = 1.0E30 GPB1F403.40
ENDDO GPB1F403.41
*IF DEF,MPP GPB1F403.42
ENDIF GPB1F403.43
*ENDIF GPB1F403.44
ENDIF SFEXCH3A.1310
*ELSE SFEXCH3A.1311
C SFEXCH3A.1312
C----------------------------------------------------------------------- SFEXCH3A.1313
CL 5. Calculate "explicit" surface fluxes of momentum, then overwrite SFEXCH3A.1314
CL coefficient with GAMMA(1)*RHOKM_1 to be passed out for implicit SFEXCH3A.1315
CL calculations in P244 (subroutine IMPL_CAL). This routine only SFEXCH3A.1316
CL uses RHOKM_1when multiplied by GAMMA(1). SFEXCH3A.1317
C----------------------------------------------------------------------- SFEXCH3A.1318
C SFEXCH3A.1319
DO I=1,U_POINTS SFEXCH3A.1320
TAUX_1(I) = RHOKM_1(I) * ( U_1(I) - U_0(I) ) ! P243.132 SFEXCH3A.1321
TAUY_1(I) = RHOKM_1(I) * ( V_1(I) - V_0(I) ) ! P243.133 SFEXCH3A.1322
RHOKM_1(I) = GAMMA(1) * RHOKM_1(I) SFEXCH3A.1323
ENDDO SFEXCH3A.1324
*ENDIF SFEXCH3A.1325
6 CONTINUE ! Branch for error exit. SFEXCH3A.1326
IF (LTIMER) THEN SFEXCH3A.1327
CALL TIMER('SFEXCH ',4) SFEXCH3A.1328
ENDIF SFEXCH3A.1329
RETURN SFEXCH3A.1330
END SFEXCH3A.1331
*ENDIF SFEXCH3A.1332