*IF DEF,RECON RCINITW1.2
C *****************************COPYRIGHT****************************** RCINITW1.3
C (c) CROWN COPYRIGHT 1996, METEOROLOGICAL OFFICE, All Rights Reserved. RCINITW1.4
C RCINITW1.5
C Use, duplication or disclosure of this code is subject to the RCINITW1.6
C restrictions as set forth in the contract. RCINITW1.7
C RCINITW1.8
C Meteorological Office RCINITW1.9
C London Road RCINITW1.10
C BRACKNELL RCINITW1.11
C Berkshire UK RCINITW1.12
C RG12 2SZ RCINITW1.13
C RCINITW1.14
C If no contract has been raised with this copy of the code, the use, RCINITW1.15
C duplication or disclosure of it is strictly prohibited. Permission RCINITW1.16
C to do so must first be obtained in writing from the Head of Numerical RCINITW1.17
C Modelling at the above address. RCINITW1.18
C ******************************COPYRIGHT****************************** RCINITW1.19
!+ Initialise vertical velocity in reconfiguration for Charney-Phil grid RCINITW1.20
! RCINITW1.21
! Subroutine Interface: RCINITW1.22
SUBROUTINE RC_INIT_W(rowlength,rows,levels, 1,2RCINITW1.23
& firstlat,firstlong,dlat,dlong, RCINITW1.24
& Rho_th,Rho_p, RCINITW1.25
& U,V,W, RCINITW1.26
& R_W_Levels,R_P_Levels, RCINITW1.27
& Global, RCINITW1.28
& Errorstatus) RCINITW1.29
RCINITW1.30
IMPLICIT NONE RCINITW1.31
! RCINITW1.32
! Description: RCINITW1.33
! Generate a 3D field of vertical velocity w from 3D fields of RCINITW1.34
! u,v wind components and density on theta and pressure levels. RCINITW1.35
! This provides an initialisation of w as a prognostic variable RCINITW1.36
! in the new dynamics scheme (Charney-Phillips vertical staggering RCINITW1.37
! and Arakawa C grid). An equatorial lat-long grid is assumed. RCINITW1.38
! RCINITW1.39
! Method: RCINITW1.40
! Calculate vertical velocity from an equation derived from the RCINITW1.41
! continuity equation with w=0 at the upper boundary. RCINITW1.42
! Horizontal spacings dx and dy are obtained from the grid descriptors RCINITW1.43
! and vertical spacings from the radial vertical co-ordinate. w is RCINITW1.44
! found by integration (summation) downwards from the top level where RCINITW1.45
! w=0 provides a boundary condition. RCINITW1.46
! Note that v is only defined for (rows-1) rows. du/dx and dv/dy are RCINITW1.47
! assumed to be zero for polar rows. RCINITW1.48
! RCINITW1.49
! <--- dx(2,1) ----> RCINITW1.50
! p(1,1) u(1,1) p(2,1) u(2,1) RCINITW1.51
! v(1,1) v(2,1)^ ^ RCINITW1.52
! p(1,2) u(1,2) p(2,2)^ u(2,2) ^ = dy(2,2) RCINITW1.53
! v(1,2) v(2,2)^ ^ RCINITW1.54
! RCINITW1.55
! Current Code Owner: I Edmond RCINITW1.56
! RCINITW1.57
! History: RCINITW1.58
! Version Date Comment RCINITW1.59
! ------- ---- ------- RCINITW1.60
! 4.1 15/6/96 Original code. Ian Edmond RCINITW1.61
! RCINITW1.62
! Code Description: RCINITW1.63
! Language: FORTRAN 77 + common extensions. RCINITW1.64
! This code is written to UMDP3 v6 programming standards. RCINITW1.65
! RCINITW1.66
! System component covered: <appropriate code> RCINITW1.67
! System Task: <appropriate code> RCINITW1.68
! RCINITW1.69
! Declarations: RCINITW1.70
! These are of the form:- RCINITW1.71
! INTEGER ExampleVariable !Description of variable RCINITW1.72
! RCINITW1.73
! 1.0 Global variables (*Called COMDECKs etc...): RCINITW1.74
*CALL C_PI
RCINITW1.75
RCINITW1.76
! 2.0 Subroutine arguments RCINITW1.77
! 2.1 Scalar arguments with intent(in): RCINITW1.78
INTEGER RCINITW1.79
& RowLength ! Number of points along a horizontal row. RCINITW1.80
&,Rows ! Number of rows RCINITW1.81
&,Levels ! Number of vertical levels. RCINITW1.82
RCINITW1.83
REAL RCINITW1.84
& Firstlat ! Latitude (degrees) of first point in field. RCINITW1.85
&,Firstlong ! Longitude (degrees) of first point in field. RCINITW1.86
&,DLat ! Latitude (degrees) grid spacing. RCINITW1.87
&,DLong ! Longitude (degrees) grid spacing. RCINITW1.88
RCINITW1.89
LOGICAL RCINITW1.90
& Global ! =T global grid, =F limited area grid RCINITW1.91
RCINITW1.92
! 2.2 Array arguments with intent(in): RCINITW1.93
REAL RCINITW1.94
& U(RowLength,Rows,Levels) ! u wind component on C grid RCINITW1.95
&,Rho_th(RowLength,Rows,Levels) ! dry density on C grid theta level RCINITW1.96
&,Rho_p(RowLength,Rows,Levels) ! dry density on C grid pressure le RCINITW1.97
&,Ru(RowLength,Rows) ! hgt of u wind component on C grid RCINITW1.98
&,V(RowLength,Rows-1,Levels) ! v wind component on C grid RCINITW1.99
&,Rv(RowLength,Rows-1) ! hgt of v wind component on C g RCINITW1.100
&,R_W_Levels(RowLength,Rows,Levels)! radius on w-levels RCINITW1.101
&,R_P_Levels(RowLength,Rows,Levels)! radius on p-levels RCINITW1.102
RCINITW1.103
! 2.3 Scalar arguments with intent(InOut): RCINITW1.104
RCINITW1.105
! 2.5 Scalar arguments with intent(out): RCINITW1.106
RCINITW1.107
! 2.6 Array arguments with intent(out): RCINITW1.108
REAL RCINITW1.109
* W(RowLength,Rows,Levels) ! w wind component RCINITW1.110
RCINITW1.111
! 2.7 ErrorStatus <Delete If ErrorStatus not used> RCINITW1.112
INTEGER ErrorStatus ! Error flag (0 = OK) RCINITW1.113
RCINITW1.114
! 3.0 Local parameters: RCINITW1.115
RCINITW1.116
! 4.0 Local scalars: RCINITW1.117
INTEGER i,j,k RCINITW1.118
REAL RCINITW1.119
& Latu ! Latitude (radians) of u points. RCINITW1.120
&,Latv1 ! Latitude (radians) of v points. RCINITW1.121
&,Latv2 ! Latitude (radians) of v+1 points. RCINITW1.122
&,pole1 ! Increment of W on the top row. RCINITW1.123
&,pole2 ! Increment of W on the bottom row. RCINITW1.124
RCINITW1.125
! 5.0 l dynamic arrays: RCINITW1.126
REAL RCINITW1.127
& Dx(Rows) ! dx horizontal increment RCINITW1.128
&,Dy(RowLength) ! dy horizontal increment RCINITW1.129
&,DrhoxuByDx(RowLength,Rows) ! du/dx u wind gradient RCINITW1.130
&,DrhoxvByDy(RowLength,Rows) ! dv/dy v wind gradient RCINITW1.131
RCINITW1.132
! Function & Subroutine calls: RCINITW1.133
RCINITW1.134
!- End of header RCINITW1.135
RCINITW1.136
! 1. Calculate x and y increments RCINITW1.137
RCINITW1.138
Do k=Levels,2,-1 RCINITW1.139
! Find the height at each u point on the Charney Phillips grid. RCINITW1.140
Call cu_hgt
(R_P_Levels, !(IN) Heights of pressure surfaces RCINITW1.141
& Ru, !(OUT))Height of u points on CP gr RCINITW1.142
& RowLength*Rows, !(IN) No. of points per level. RCINITW1.143
& Levels, !(IN) No. of levels. RCINITW1.144
& k, RCINITW1.145
& RowLength, !(IN) No. of columns. RCINITW1.146
& (.NOT.Global)) !(IN) Global or LAM grid. RCINITW1.147
RCINITW1.148
! Find the height at each v point on the Charney Phillips grid. RCINITW1.149
Call cv_hgt
(R_P_Levels, !(IN) Heights of pressure surfaces RCINITW1.150
& Rv, !(OUT)Height of v points on CP gri RCINITW1.151
& RowLength*Rows, !(IN) No. of points per level. RCINITW1.152
& Levels, !(IN) No. of levels. RCINITW1.153
& k, RCINITW1.154
& RowLength, !(IN) No. of columns. RCINITW1.155
& (.NOT.Global)) !(IN) Global or LAM grid. RCINITW1.156
RCINITW1.157
! 2. Calculate w from integration of horizontal wind gradients. RCINITW1.158
RCINITW1.159
! 2.1 Find d(Rho_p * u)/dx RCINITW1.160
RCINITW1.161
! General case d(Rho_p*u)/dx RCINITW1.162
Do j=2,Rows-1 RCINITW1.163
RCINITW1.164
Do i=2,RowLength-1 RCINITW1.165
RCINITW1.166
Latu = ( (j-1) * DLat + FirstLat ) * pi_over_180 RCINITW1.167
RCINITW1.168
DrhoxuByDx(i,j) = ( RCINITW1.169
& Ru(i,j) * U(i,j,k) RCINITW1.170
& * ( Rho_p(i+1,j,k) + Rho_p(i,j,k) ) RCINITW1.171
& - Ru(i-1,j) * U(i-1,j,k) RCINITW1.172
& * ( Rho_p(i,j,k) + Rho_p(i-1,j,k) ) RCINITW1.173
& ) RCINITW1.174
& / ( RCINITW1.175
& 2 * R_P_Levels(i,j,k)**2 RCINITW1.176
& * cos(Latu) * DLong * pi_over_180 RCINITW1.177
& ) RCINITW1.178
RCINITW1.179
RCINITW1.180
End do ! i 2,RowLength RCINITW1.181
RCINITW1.182
End do ! j 2,Rows-1 RCINITW1.183
RCINITW1.184
! First E-W boundary column d(Rho_p*u)/dx RCINITW1.185
RCINITW1.186
Do j=2,Rows-1 RCINITW1.187
RCINITW1.188
Latu = ( (j-1) * DLat + FirstLat ) * pi_over_180 RCINITW1.189
RCINITW1.190
DrhoxuByDx(1,j) = ( RCINITW1.191
& Ru(1,j) * U(1,j,k) RCINITW1.192
& * ( Rho_p(2,j,k) + Rho_p(1,j,k) ) RCINITW1.193
& - Ru(RowLength,j) * U(RowLength,j,k) RCINITW1.194
& * ( Rho_p(1,j,k) + Rho_p(RowLength,j,k) ) RCINITW1.195
& ) RCINITW1.196
& / ( RCINITW1.197
& 2 * R_P_Levels(1,j,k)**2 RCINITW1.198
& * cos(Latu) * DLong * pi_over_180 RCINITW1.199
& ) RCINITW1.200
RCINITW1.201
DrhoxuByDx(RowLength,j) = RCINITW1.202
& ( RCINITW1.203
& Ru(RowLength,j) * U(RowLength,j,k) RCINITW1.204
& * ( Rho_p(1,j,k) + Rho_p(RowLength,j,k) ) RCINITW1.205
& - Ru(RowLength-1,j) * U(RowLength-1,j,k) RCINITW1.206
& * ( Rho_p(RowLength,j,k) + Rho_p(RowLength-1,j,k) ) RCINITW1.207
& ) RCINITW1.208
& / ( RCINITW1.209
& 2 * R_P_Levels(RowLength,j,k)**2 RCINITW1.210
& * cos(Latu) * DLong * pi_over_180 RCINITW1.211
& ) RCINITW1.212
RCINITW1.213
End do ! j 2,Rows-1 RCINITW1.214
RCINITW1.215
Do i=2,Rowlength-1 RCINITW1.216
RCINITW1.217
DrhoxuByDx(i,1) = 0.0 RCINITW1.218
DrhoxuByDx(i,rows) = 0.0 RCINITW1.219
RCINITW1.220
End do ! i 2,Rowlength-1 RCINITW1.221
RCINITW1.222
! First and last N-S rows d(Rho_p*u)/dx RCINITW1.223
RCINITW1.224
Do i=1,RowLength RCINITW1.225
RCINITW1.226
DrhoxuByDx(i,1) = 0.0 RCINITW1.227
DrhoxuByDx(i,Rows) = 0.0 RCINITW1.228
RCINITW1.229
End do ! i 1,RowLength RCINITW1.230
RCINITW1.231
RCINITW1.232
! 2.2 Find d(Rho_p * v)/dy RCINITW1.233
RCINITW1.234
! General case d(Rho_p*v)/dy RCINITW1.235
Do j=2,Rows-1 RCINITW1.236
RCINITW1.237
Do i=1,RowLength RCINITW1.238
RCINITW1.239
Latu = ( (j-1) * DLat + FirstLat ) * pi_over_180 RCINITW1.240
Latv1 = ( (j-1) * DLat + FirstLat - (DLat/2.0) ) RCINITW1.241
& * pi_over_180 RCINITW1.242
Latv2 = ( (j-1) * DLat + FirstLat + (DLat/2.0) ) RCINITW1.243
& * pi_over_180 RCINITW1.244
RCINITW1.245
DrhoxvByDy(i,j) = RCINITW1.246
& ( RCINITW1.247
& Rv(i,j-1) * V(i,j-1,k) * cos(Latv1) RCINITW1.248
& * ( Rho_p(i,j,k) + Rho_p(i,j-1,k) ) RCINITW1.249
& - Rv(i,j) * V(i,j,k) * cos(Latv2) RCINITW1.250
& * ( Rho_p(i,j,k) + Rho_p(i,j+1,k) ) RCINITW1.251
& ) RCINITW1.252
& / ( RCINITW1.253
& 2 * R_P_Levels(i,j,k)**2 RCINITW1.254
& * cos(Latu) * DLat * pi_over_180 RCINITW1.255
& ) RCINITW1.256
RCINITW1.257
End do ! i 1,RowLength RCINITW1.258
RCINITW1.259
End do ! j 2,Rows-1 RCINITW1.260
RCINITW1.261
! First and last N-S boundary rows d(Rho_p*v)/dy RCINITW1.262
Do i=1,RowLength RCINITW1.263
RCINITW1.264
Latu = 0.125 * ( DLat + FirstLat ) * pi_over_180 RCINITW1.265
Latv1 = ( FirstLat + (DLat/2.0) ) * pi_over_180 RCINITW1.266
Latv2 = ( (Rows-1) * DLat + FirstLat- (DLat/2.0) ) RCINITW1.267
& * pi_over_180 RCINITW1.268
RCINITW1.269
DrhoxvByDy(i,1) = RCINITW1.270
& ( RCINITW1.271
& -Rv(i,1) * V(i,1,k) * cos(Latv1) RCINITW1.272
& * ( Rho_p(i,2,k) + Rho_p(i,1,k) ) RCINITW1.273
& ) RCINITW1.274
& / ( RCINITW1.275
& 2 * R_P_Levels(i,1,k)**2 RCINITW1.276
& * cos(Latu) * DLat * pi_over_180 RCINITW1.277
& ) RCINITW1.278
RCINITW1.279
Latu = 0.125 * ( (Rows-2) * DLat + FirstLat) RCINITW1.280
& * pi_over_180 RCINITW1.281
RCINITW1.282
DrhoxvByDy(i,Rows) = RCINITW1.283
& ( RCINITW1.284
& Rv(i,Rows-1) * V(i,Rows-1,k) * cos(Latv2) RCINITW1.285
& * ( Rho_p(i,Rows,k) + Rho_p(i,Rows-1,k) ) RCINITW1.286
& ) RCINITW1.287
& / ( RCINITW1.288
& 2 * R_P_Levels(i,Rows,k)**2 RCINITW1.289
& * cos(Latu) * DLat * pi_over_180 RCINITW1.290
& ) RCINITW1.291
RCINITW1.292
End do ! i 1,RowLength RCINITW1.293
RCINITW1.294
! 2.3 W found using an equation obtained from VAR doc. UMDP101 part 3. RCINITW1.295
RCINITW1.296
pole1 = 0.0 RCINITW1.297
pole2 = 0.0 RCINITW1.298
RCINITW1.299
! Top level W is 0. RCINITW1.300
If (k .EQ. Levels) then RCINITW1.301
RCINITW1.302
Do j=1,Rows RCINITW1.303
RCINITW1.304
Do i=1,RowLength RCINITW1.305
RCINITW1.306
W(i,j,k) = 0.0 RCINITW1.307
RCINITW1.308
End do ! i 1,RowLength RCINITW1.309
RCINITW1.310
End do ! j 1,Rows RCINITW1.311
RCINITW1.312
End if RCINITW1.313
RCINITW1.314
RCINITW1.315
! W at level k-1. RCINITW1.316
Do j=2,Rows-1 RCINITW1.317
RCINITW1.318
Do i=1,RowLength RCINITW1.319
RCINITW1.320
W(i,j,k-1) = ( Rho_th(i,j,k) / Rho_th(i,j,k-1) ) * W(i,j,k) RCINITW1.321
& + ( ( R_W_Levels(i,j,k) - R_W_Levels(i,j,k-1) ) RCINITW1.322
& / Rho_th(i,j,k-1) RCINITW1.323
& ) * ( DrhoxuByDx(i,j) + DrhoxvByDy(i,j) ) RCINITW1.324
& RCINITW1.325
RCINITW1.326
End do ! i 1,RowLength RCINITW1.327
RCINITW1.328
End do ! j 1,Rows RCINITW1.329
RCINITW1.330
RCINITW1.331
! First and last N-S rows W at level k-1. RCINITW1.332
If (Global) then RCINITW1.333
RCINITW1.334
Do i=1,RowLength RCINITW1.335
RCINITW1.336
pole1 = pole1 RCINITW1.337
& + ( RCINITW1.338
& ( Rho_th(i,1,k) / Rho_th(i,1,k-1) ) * W(i,1,k) RCINITW1.339
& + ( ( R_W_Levels(i,1,k) UIE2F404.1384
& - R_W_Levels(i,1,k-1) UIE2F404.1385
& ) / Rho_th(i,1,k-1) UIE2F404.1386
& ) * ( DrhoxuByDx(i,1) + DrhoxvByDy(i,1) ) RCINITW1.342
& ) RCINITW1.343
RCINITW1.344
pole2 = pole2 RCINITW1.345
& + ( RCINITW1.346
& ( Rho_th(i,Rows,k) / Rho_th(i,Rows,k-1) ) RCINITW1.347
& * W(i,Rows,k) RCINITW1.348
& + ( ( R_W_Levels(i,Rows,k) RCINITW1.349
& - R_W_Levels(i,Rows,k-1) RCINITW1.350
& ) / Rho_th(i,Rows,k-1) RCINITW1.351
& ) * ( DrhoxuByDx(i,Rows) + DrhoxvByDy(i,Rows) ) RCINITW1.352
& ) RCINITW1.353
RCINITW1.354
End do ! i 1,RowLength RCINITW1.355
RCINITW1.356
! If global data then W at poles = average of W values at each RCINITW1.357
Do i=1,RowLength RCINITW1.358
RCINITW1.359
W(i,1,k-1) = pole1 / RowLength RCINITW1.360
W(i,Rows,k-1) = pole2 / RowLength RCINITW1.361
RCINITW1.362
End do ! i 1,RowLength RCINITW1.363
RCINITW1.364
RCINITW1.365
else ! LAM RCINITW1.366
RCINITW1.367
! If LAM data W at poles = 0.0 RCINITW1.368
Do i=1,RowLength RCINITW1.369
RCINITW1.370
W(i,1,k-1) = 0.0 RCINITW1.371
W(i,Rows,k-1) = 0.0 RCINITW1.372
RCINITW1.373
End do ! i RCINITW1.374
RCINITW1.375
Do j=1,Rows RCINITW1.376
RCINITW1.377
W(1,j,k-1) = 0.0 RCINITW1.378
W(Rowlength,j,k-1) = 0.0 RCINITW1.379
RCINITW1.380
End do RCINITW1.381
RCINITW1.382
RCINITW1.383
End if RCINITW1.384
RCINITW1.385
End do ! k RCINITW1.386
RCINITW1.387
RETURN RCINITW1.388
END RCINITW1.389
*ENDIF RCINITW1.390