*IF DEF,A17_1A GRVSET1A.2
C *****************************COPYRIGHT****************************** GRVSET1A.3
C (c) CROWN COPYRIGHT 1997, METEOROLOGICAL OFFICE, All Rights Reserved. GRVSET1A.4
C GRVSET1A.5
C Use, duplication or disclosure of this code is subject to the GRVSET1A.6
C restrictions as set forth in the contract. GRVSET1A.7
C GRVSET1A.8
C Meteorological Office GRVSET1A.9
C London Road GRVSET1A.10
C BRACKNELL GRVSET1A.11
C Berkshire UK GRVSET1A.12
C RG12 2SZ GRVSET1A.13
C GRVSET1A.14
C If no contract has been raised with this copy of the code, the use, GRVSET1A.15
C duplication or disclosure of it is strictly prohibited. Permission GRVSET1A.16
C to do so must first be obtained in writing from the Head of Numerical GRVSET1A.17
C Modelling at the above address. GRVSET1A.18
C ******************************COPYRIGHT****************************** GRVSET1A.19
! Subroutine GRAVSETT ---------------------------------------------- GRVSET1A.20
! GRVSET1A.21
! Purpose: To perform gravitational settlement of tracer particles GRVSET1A.22
! down to the lowest layer of the model. GRVSET1A.23
! This version allows tracers to fall through 1 or 2 layers. GRVSET1A.24
! GRVSET1A.25
! Current owners of code: S Woodward, M Woodage GRVSET1A.26
! GRVSET1A.27
! History: GRVSET1A.28
! Version Date Comment GRVSET1A.29
! ------- ---- ------- GRVSET1A.30
! 4.4 03/10/97 Original code S Woodward, M Woodage GRVSET1A.31
! GRVSET1A.32
! Code description: GRVSET1A.33
! Language: FORTRAN77 + extensions GRVSET1A.34
! Programming standard: UMDP 3 Vn 6 GRVSET1A.35
! GRVSET1A.36
! System components covered: GRVSET1A.37
! GRVSET1A.38
! System task: GRVSET1A.39
! GRVSET1A.40
!Documentation: Not yet available GRVSET1A.41
! GRVSET1A.42
!----------------------------------------------------------------------- GRVSET1A.43
! GRVSET1A.44
SUBROUTINE GRAVSETT( 2,1GRVSET1A.45
& PFIELD,NLEVS,TRACFLD,DIAM,RHOP,PSTAR,AK,BK,DELTA_AK,DELTA_BK,T, GRVSET1A.46
& FIRST_POINT,LAST_POINT,DT,DRYDEP) GRVSET1A.47
! GRVSET1A.48
! GRVSET1A.49
IMPLICIT NONE GRVSET1A.50
! GRVSET1A.51
INTEGER NLEVS !IN number of model levels GRVSET1A.52
INTEGER PFIELD !IN number of grid points GRVSET1A.53
INTEGER FIRST_POINT !IN first point for calcns to be done GRVSET1A.54
INTEGER LAST_POINT !IN last point for calcns to be done GRVSET1A.55
! GRVSET1A.56
REAL DIAM !IN tracer particle diameter GRVSET1A.57
REAL RHOP !IN tracer particle density GRVSET1A.58
REAL PSTAR(PFIELD) !IN Pstar GRVSET1A.59
REAL AK(NLEVS) !IN A vals on levs GRVSET1A.60
REAL BK(NLEVS) !IN B vals on levs GRVSET1A.61
REAL DELTA_AK(NLEVS) !IN A(lev+1/2) - A(lev-1/2) GRVSET1A.62
REAL DELTA_BK(NLEVS) !IN B(lev+1/2) - B(lev-1/2) GRVSET1A.63
REAL T(PFIELD,NLEVS) !IN temperature GRVSET1A.64
REAL DT !IN timestep s GRVSET1A.65
! GRVSET1A.66
REAL TRACFLD(PFIELD,NLEVS) !IN/OUT tracer field GRVSET1A.67
! GRVSET1A.68
REAL DRYDEP(PFIELD) !OUT deposition flux from layer2 (kg m-2 s-1) GRVSET1A.69
! GRVSET1A.70
GRVSET1A.71
! Include COMDECKS GRVSET1A.72
! GRVSET1A.73
*IF DEF,MPP GRVSET1A.74
! Parameters and Common blocks GRVSET1A.75
*CALL PARVARS 
GRVSET1A.76
*ENDIF GRVSET1A.77
! GRVSET1A.78
*CALL C_R_CP GRVSET1A.79
*CALL C_G GRVSET1A.80
*CALL C_SULCHM GRVSET1A.81
! GRVSET1A.82
! External subroutines called GRVSET1A.83
EXTERNAL VGRAV GRVSET1A.84
! GRVSET1A.85
! Local variables GRVSET1A.86
! GRVSET1A.87
INTEGER K !LOC loop counter for levels GRVSET1A.88
INTEGER J !LOC loop counter for points GRVSET1A.89
! GRVSET1A.90
REAL P(PFIELD) !LOC Pressure GRVSET1A.91
REAL VRHOCDT(PFIELD) !LOC v*rho*tracer*deltat @lev GRVSET1A.92
REAL RHOK2(PFIELD) !LOC rho(lev+2) GRVSET1A.93
REAL RHOK1(PFIELD) !LOC rho(lev+1) GRVSET1A.94
REAL RHOK(PFIELD) !LOC rho(lev) GRVSET1A.95
REAL DZK(PFIELD) !LOC thickness of layer lev GRVSET1A.96
REAL DZK1(PFIELD) !LOC thickness of layer lev+1 GRVSET1A.97
REAL DZK2(PFIELD) !LOC thickness of layer lev+2 GRVSET1A.98
REAL V(PFIELD,NLEVS) !LOC deposition velocity (vstokes corrected) GRVSET1A.99
REAL MASSOUT2K2(PFIELD) !LOC flux falling 2 levs from lev k+2 GRVSET1A.100
REAL MASSOUT1K2(PFIELD) !LOC flux falling 1 levs from lev k+2 GRVSET1A.101
REAL MASSOUT2K1(PFIELD) !LOC flux falling 2 levs from lev k+1 GRVSET1A.102
REAL MASSOUT1K1(PFIELD) !LOC flux falling 1 levs from lev k+1 GRVSET1A.103
REAL MASSOUT2K(PFIELD) !LOC flux falling 2 levs from lev k GRVSET1A.104
REAL MASSOUT1K(PFIELD) !LOC flux falling 1 levs from lev k GRVSET1A.105
REAL DUMMY1(PFIELD,NLEVS) !LOC GRVSET1A.106
REAL DUMMY2(PFIELD,NLEVS) !LOC GRVSET1A.107
! GRVSET1A.108
! GRVSET1A.109
! Calculate settlement velocity GRVSET1A.110
! GRVSET1A.111
CALL VGRAV(PFIELD,NLEVS,DIAM,RHOP,PSTAR,AK,BK,T,V,DUMMY1,DUMMY2, GRVSET1A.112
& FIRST_POINT,LAST_POINT) GRVSET1A.113
! GRVSET1A.114
! Calculate new tracer mixing ratios GRVSET1A.115
! GRVSET1A.116
! Initialise deposition flux to zero GRVSET1A.117
DO J = 1,PFIELD GRVSET1A.118
DRYDEP(J)=0. GRVSET1A.119
ENDDO GRVSET1A.120
! GRVSET1A.121
! Level 1 (K at start of loop) GRVSET1A.122
! GRVSET1A.123
DO J = FIRST_POINT,LAST_POINT GRVSET1A.124
P(J)=AK(1)+BK(1)*PSTAR(J) GRVSET1A.125
RHOK(J)=P(J)/(R*T(J,1)) GRVSET1A.126
DZK(J)=-(DELTA_AK(1)+DELTA_BK(1)*PSTAR(J))/(RHOK(J)*G) GRVSET1A.127
MASSOUT2K(J)=0. GRVSET1A.128
MASSOUT1K(J)=0. GRVSET1A.129
ENDDO !J GRVSET1A.130
! GRVSET1A.131
! Level 2 (K+1 at start of loop) GRVSET1A.132
! NB deposit tracer direct to ground from lev 2 if V high enough GRVSET1A.133
! GRVSET1A.134
DO J = FIRST_POINT,LAST_POINT GRVSET1A.135
! GRVSET1A.136
P(J)=AK(2)+BK(2)*PSTAR(J) GRVSET1A.137
RHOK1(J)=P(J)/(R*T(J,2)) GRVSET1A.138
DZK1(J)=-(DELTA_AK(2)+DELTA_BK(2)*PSTAR(J))/(RHOK1(J)*G) GRVSET1A.139
! GRVSET1A.140
! check for deposition : GRVSET1A.141
IF (V(J,2)*DT .GT. DZK(J)) THEN GRVSET1A.142
! some tracer deposited onto ground GRVSET1A.143
! GRVSET1A.144
IF ( V(J,2)*DT .GT. DZK1(J)+DZK(J) ) THEN GRVSET1A.145
! all deposited to ground GRVSET1A.146
MASSOUT2K1(J)=RHOK1(J)*TRACFLD(J,2)*DZK1(J) GRVSET1A.147
MASSOUT1K1(J)=0. GRVSET1A.148
ELSE IF ( V(J,2)*DT .GT. DZK1(J) ) THEN GRVSET1A.149
! some deposited to ground, some to layer 1 GRVSET1A.150
MASSOUT2K1(J)=RHOK1(J)*TRACFLD(J,2)*(V(J,2)*DT-DZK(J)) GRVSET1A.151
MASSOUT1K1(J)=RHOK1(J)*TRACFLD(J,2)*DZK1(J)-MASSOUT2K1(J) GRVSET1A.152
ELSE GRVSET1A.153
! some deposited to ground, some to layer1, some left in layer2 GRVSET1A.154
MASSOUT2K1(J)=RHOK1(J)*TRACFLD(J,2)*(V(J,2)*DT-DZK(J)) GRVSET1A.155
MASSOUT1K1(J)=RHOK1(J)*TRACFLD(J,2)*DZK(J) GRVSET1A.156
ENDIF GRVSET1A.157
! GRVSET1A.158
DRYDEP(J)=MASSOUT2K1(J)/DT GRVSET1A.159
! GRVSET1A.160
ELSE GRVSET1A.161
! only falls into layer 1 GRVSET1A.162
MASSOUT2K1(J)=0. GRVSET1A.163
IF ( V(J,2)*DT .GT. DZK1(J)) THEN GRVSET1A.164
! all falls into layer 1 GRVSET1A.165
MASSOUT1K1(J)=RHOK1(J)*TRACFLD(J,2)*DZK1(J) GRVSET1A.166
ELSE GRVSET1A.167
! some to layer 1 , some left in layer2 GRVSET1A.168
MASSOUT1K1(J)=RHOK1(J)*TRACFLD(J,2)*V(J,2)*DT GRVSET1A.169
ENDIF GRVSET1A.170
! GRVSET1A.171
ENDIF GRVSET1A.172
! GRVSET1A.173
ENDDO !END J LOOP GRVSET1A.174
! GRVSET1A.175
! Main loop through levels, from bottom up GRVSET1A.176
! GRVSET1A.177
DO K = 1,NLEVS-2 GRVSET1A.178
! GRVSET1A.179
DO J = FIRST_POINT,LAST_POINT GRVSET1A.180
! GRVSET1A.181
P(J)=AK(K+2)+BK(K+2)*PSTAR(J) GRVSET1A.182
RHOK2(J)=P(J)/(R*T(J,K+2)) GRVSET1A.183
DZK2(J)=-(DELTA_AK(K+2)+DELTA_BK(K+2)*PSTAR(J))/ GRVSET1A.184
& (RHOK2(J)*G) GRVSET1A.185
! GRVSET1A.186
! Calculate mass of tracer falling between levels GRVSET1A.187
! GRVSET1A.188
! limit fall to 2 levs GRVSET1A.189
IF (V(J,K+2)*DT.GT.(DZK1(J)+DZK(J))) GRVSET1A.190
& V(J,K+2)=(DZK1(J)+DZK(J))/DT GRVSET1A.191
! GRVSET1A.192
! check how far tracer falls: GRVSET1A.193
IF ( V(J,K+2)*DT .GT. DZK1(J) ) THEN GRVSET1A.194
! it falls through more than 1 layer GRVSET1A.195
IF ( V(J,K+2)*DT .GT. (DZK2(J)+DZK1(J)) ) THEN GRVSET1A.196
! all into layer k GRVSET1A.197
MASSOUT2K2(J)=RHOK2(J)*TRACFLD(J,K+2)*DZK2(J) GRVSET1A.198
MASSOUT1K2(J)=0. GRVSET1A.199
ELSE IF ( V(J,K+2)*DT .GT. DZK2(J) ) THEN GRVSET1A.200
! some into k+1, some into k GRVSET1A.201
MASSOUT2K2(J)= GRVSET1A.202
& RHOK2(J)*TRACFLD(J,K+2)*(V(J,K+2)*DT-DZK1(J)) GRVSET1A.203
MASSOUT1K2(J)= GRVSET1A.204
& RHOK2(J)*TRACFLD(J,K+2)*DZK2(J)-MASSOUT2K2(J) GRVSET1A.205
ELSE GRVSET1A.206
! some left in k+2, some into k+1, some into k GRVSET1A.207
MASSOUT2K2(J)= GRVSET1A.208
& RHOK2(J)*TRACFLD(J,K+2)*(V(J,K+2)*DT-DZK1(J)) GRVSET1A.209
MASSOUT1K2(J)=RHOK2(J)*TRACFLD(J,K+2)*DZK1(J) GRVSET1A.210
ENDIF GRVSET1A.211
! GRVSET1A.212
ELSE GRVSET1A.213
! falls no more than 1 layer GRVSET1A.214
MASSOUT2K2(J)=0. GRVSET1A.215
IF (V(J,K+2)*DT .GT. DZK2(J)) THEN GRVSET1A.216
! all falls into layer k+1 GRVSET1A.217
MASSOUT1K2(J)=RHOK2(J)*TRACFLD(J,K+2)*DZK2(J) GRVSET1A.218
ELSE GRVSET1A.219
! some falls into k+1, some left in k+2 GRVSET1A.220
MASSOUT1K2(J)=RHOK2(J)*TRACFLD(J,K+2)*V(J,K+2)*DT GRVSET1A.221
ENDIF GRVSET1A.222
! GRVSET1A.223
ENDIF GRVSET1A.224
! GRVSET1A.225
! Update tracer field GRVSET1A.226
! GRVSET1A.227
TRACFLD(J,K)=TRACFLD(J,K)+(MASSOUT2K2(J)+MASSOUT1K1(J)- GRVSET1A.228
& MASSOUT2K(J)-MASSOUT1K(J))/(RHOK(J)*DZK(J)) GRVSET1A.229
! GRVSET1A.230
! Put k+2 vals in k+1's & k+1's in k's GRVSET1A.231
MASSOUT1K(J)=MASSOUT1K1(J) GRVSET1A.232
MASSOUT1K1(J)=MASSOUT1K2(J) GRVSET1A.233
MASSOUT2K(J)=MASSOUT2K1(J) GRVSET1A.234
MASSOUT2K1(J)=MASSOUT2K2(J) GRVSET1A.235
DZK(J)=DZK1(J) GRVSET1A.236
DZK1(J)=DZK2(J) GRVSET1A.237
RHOK(J)=RHOK1(J) GRVSET1A.238
RHOK1(J)=RHOK2(J) GRVSET1A.239
! GRVSET1A.240
ENDDO !END J LOOP GRVSET1A.241
! GRVSET1A.242
ENDDO !END K LOOP GRVSET1A.243
! GRVSET1A.244
! Top 2 levels GRVSET1A.245
! GRVSET1A.246
DO J=FIRST_POINT,LAST_POINT GRVSET1A.247
! GRVSET1A.248
TRACFLD(J,NLEVS-1)=TRACFLD(J,NLEVS-1)+ GRVSET1A.249
& (MASSOUT1K1(J)-MASSOUT2K(J)-MASSOUT1K(J))/(RHOK(J)*DZK(J)) GRVSET1A.250
TRACFLD(J,NLEVS)=TRACFLD(J,NLEVS)- GRVSET1A.251
& (MASSOUT2K1(J)+MASSOUT1K1(J))/ GRVSET1A.252
& (RHOK1(J)*DZK1(J)) GRVSET1A.253
! GRVSET1A.254
ENDDO !J GRVSET1A.255
! GRVSET1A.256
RETURN GRVSET1A.257
END GRVSET1A.258
*ENDIF GRVSET1A.259