Changeset 13489 for branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1_v2/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90
- Timestamp:
- 2020-09-17T19:13:57+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1_v2/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90
r12352 r13489 15 15 USE trc 16 16 USE par_fabm 17 USE fabm18 17 USE dom_oce 19 18 #if defined key_trdtrc && defined key_iomput … … 25 24 26 25 # include "domzgr_substitute.h90" 26 # include "vectopt_loop_substitute.h90" 27 27 28 28 PRIVATE … … 32 32 ! Work arrays for vertical advection (residual movement/sinking/floating) 33 33 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:) :: w_ct 34 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: w_if35 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: zwgt_if36 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: flux_if37 34 #if defined key_trdtrc && defined key_iomput 38 35 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:,:) :: tr_vmv … … 41 38 CONTAINS 42 39 43 SUBROUTINE compute_vertical_movement( kt )40 SUBROUTINE compute_vertical_movement( kt, method ) 44 41 !!---------------------------------------------------------------------- 45 42 !! *** compute_vertical_movement *** 46 43 !! 47 !! ** Purpose : compute vertical movement of FABM tracers 44 !! ** Purpose : compute vertical movement of FABM tracers through the water 45 !! (sinking/floating/active movement) 48 46 !! 49 !! ** Method : Sets additional vertical velocity field and computes 50 !! resulting advection using a conservative 3rd upwind 51 !! scheme with QUICKEST TVD limiter, based on GOTM 52 !! module adv_center.F90 (www.gotm.net). Currently assuming 53 !! zero flux at sea surface and sea floor. 47 !! ** Method : Retrieves additional vertical velocity field and applies 48 !! advection scheme. 54 49 !!---------------------------------------------------------------------- 55 50 ! 56 INTEGER, INTENT(in) :: kt ! ocean time-step index 57 INTEGER :: ji,jj,jk,jn,k_floor,n_iter,n_count 58 INTEGER,PARAMETER :: n_itermax=100 59 REAL(wp) :: cmax_no,z2dt 60 REAL(wp),DIMENSION(jpk) :: tr_it,tr_u,tr_d,tr_c,tr_slope,c_no,flux_lim 61 REAL(wp),DIMENSION(jpk) :: phi_lim,x_fac 51 INTEGER, INTENT(in) :: kt ! ocean time-step index 52 INTEGER, INTENT(in) :: method ! advection method (1: 1st order upstream, 3: 3rd order TVD with QUICKEST limiter) 53 54 INTEGER :: ji,jj,jk,jn,k_floor 55 REAL(wp) :: zwgt_if(1:jpkm1-1), dc(1:jpkm1), w_if(1:jpkm1-1), z2dt, h(1:jpkm1) 62 56 #if defined key_trdtrc 63 57 CHARACTER (len=20) :: cltra … … 69 63 70 64 IF( neuler == 0 .AND. kt == nittrc000 ) THEN 71 65 z2dt = rdt ! set time step size (Euler) 72 66 ELSE 73 67 z2dt = 2._wp * rdt ! set time step size (Leapfrog) 74 68 ENDIF 69 75 70 ! Compute interior vertical velocities and include them in source array. 76 DO jj= 1,jpj! j-loop77 ! Get vertical velocities at layer centres (entire 1:jpi,1:jpk slice).78 DO jk=1,jpk 79 CALL fabm_get_vertical_movement(model,1,jpi,jj,jk,w_ct(:,jk,:))71 DO jj=2,jpjm1 ! j-loop 72 ! Get vertical velocities at layer centres (entire i-k slice). 73 DO jk=1,jpkm1 74 CALL model%get_vertical_movement(fs_2,fs_jpim1,jj,jk,w_ct(:,jk,:)) 80 75 END DO 81 82 DO ji=1,jpi ! i-loop 76 DO ji=fs_2,fs_jpim1 ! i-loop 83 77 ! Only process this horizontal point (ji,jj) if number of layers exceeds 1 84 IF (mbkt(ji,jj)>1) THEN ! Level check85 k_floor=mbkt(ji,jj)78 k_floor = mbkt(ji,jj) 79 IF (k_floor > 1) THEN 86 80 ! Linearly interpolate to velocities at the interfaces between layers 87 81 ! Note: 88 ! - interface k sits between cell centre k and k -1,89 ! - k [1,jpk ] increases downwards82 ! - interface k sits between cell centre k and k+1 (k=0 for surface) 83 ! - k [1,jpkm1] increases downwards 90 84 ! - upward velocity is positive, downward velocity is negative 91 zwgt_if(1,:)=0._wp ! surface 92 w_if(1,:)=0._wp ! surface 93 zwgt_if(2:k_floor,:)=spread(& 94 fse3t(ji,jj,2:k_floor)/ (fse3t(ji,jj,1:k_floor-1)+fse3t(ji,jj,2:k_floor))& 95 ,2,jp_fabm) 96 w_if(2:k_floor,:) = zwgt_if(2:k_floor,:)*w_ct(ji,1:k_floor-1,:)& 97 +(1._wp-zwgt_if(1:k_floor-1,:))*w_ct(ji,2:k_floor,:) 98 zwgt_if(k_floor+1:,:)=0._wp ! sea floor and below 99 w_if(k_floor+1:,:)=0._wp ! sea floor and below 85 h(1:k_floor) = fse3t(ji,jj,1:k_floor) 86 zwgt_if(1:k_floor-1) = h(2:k_floor) / (h(1:k_floor-1) + h(2:k_floor)) 100 87 101 88 ! Advect: 102 89 DO jn=1,jp_fabm ! State loop 103 ! get maximum Courant number: 104 c_no(2:k_floor)=abs(w_if(2:k_floor,jn))*z2dt/ & 105 ( 0.5_wp*(fse3t(ji,jj,2:k_floor) + & 106 fse3t(ji,jj,1:k_floor-1)) ) 107 cmax_no=MAXVAL(c_no(2:k_floor)) 108 109 ! number of iterations: 110 n_iter=min(n_itermax,int(cmax_no)+1) 111 IF (ln_ctl.AND.(n_iter .gt. 1)) THEN 112 WRITE(numout,*) 'vertical_movement_fabm():' 113 WRITE(numout,*) ' Maximum Courant number is ',cmax_no,'.' 114 WRITE(numout,*) ' ',n_iter,' iterations used for vertical advection.' 115 ENDIF 116 117 ! effective Courant number: 118 c_no=c_no/n_iter 119 120 tr_it(1:k_floor)=trb(ji,jj,1:k_floor,jp_fabm_m1+jn) 121 DO n_count=1,n_iter ! Iterative loop 122 !Compute slope ratio 123 IF (k_floor.gt.2) THEN !More than 2 vertical wet points 124 IF (k_floor.gt.3) THEN 125 WHERE (w_if(3:k_floor-1,jn).ge.0._wp) !upward movement 126 tr_u(3:k_floor-1)=tr_it(4:k_floor) 127 tr_c(3:k_floor-1)=tr_it(3:k_floor-1) 128 tr_d(3:k_floor-1)=tr_it(2:k_floor-2) 129 ELSEWHERE !downward movement 130 tr_u(3:k_floor-1)=tr_it(1:k_floor-3) 131 tr_c(3:k_floor-1)=tr_it(2:k_floor-2) 132 tr_d(3:k_floor-1)=tr_it(3:k_floor-1) 133 ENDWHERE 134 ENDIF 135 IF (w_if(2,jn).ge.0._wp) THEN 136 tr_u(2)=tr_it(3) 137 tr_c(2)=tr_it(2) 138 tr_d(2)=tr_it(1) 139 ELSE 140 tr_u(2)=tr_it(1) 141 tr_c(2)=tr_it(1) 142 tr_d(2)=tr_it(2) 143 ENDIF 144 IF (w_if(k_floor,jn).ge.0._wp) THEN 145 tr_u(k_floor)=tr_it(k_floor) 146 tr_c(k_floor)=tr_it(k_floor) 147 tr_d(k_floor)=tr_it(k_floor-1) 148 ELSE 149 tr_u(k_floor)=tr_it(k_floor-2) 150 tr_c(k_floor)=tr_it(k_floor-1) 151 tr_d(k_floor)=tr_it(k_floor) 152 ENDIF 153 ELSE !only 2 vertical wet points, i.e. only 1 interface 154 IF (w_if(k_floor,jn).ge.0._wp) THEN 155 tr_u(2)=tr_it(2) 156 tr_c(2)=tr_it(2) 157 tr_d(2)=tr_it(1) 158 ELSE 159 tr_u(2)=tr_it(1) 160 tr_c(2)=tr_it(1) 161 tr_d(2)=tr_it(2) 162 ENDIF 163 ENDIF 164 WHERE (abs(tr_d(2:k_floor)-tr_c(2:k_floor)).gt.1.e-10_wp) 165 tr_slope(2:k_floor)= & 166 (tr_c(2:k_floor)-tr_u(2:k_floor))/ & 167 (tr_d(2:k_floor)-tr_c(2:k_floor)) 168 ELSEWHERE 169 tr_slope(2:k_floor)=SIGN(1._wp,w_if(2:k_floor,jn))* & 170 (tr_c(2:k_floor)-tr_u(2:k_floor))*1.e10_wp 171 ENDWHERE 172 173 !QUICKEST flux limiter: 174 x_fac(2:k_floor)=(1._wp-2._wp*c_no(2:k_floor))/6._wp 175 phi_lim(2:k_floor)=(0.5_wp+x_fac(2:k_floor)) + & 176 (0.5_wp-x_Fac(2:k_floor))*tr_slope(2:k_floor) 177 flux_lim(2:k_floor)=max( 0._wp, & 178 min( phi_lim(2:k_floor),2._wp/(1._wp-c_no(2:k_floor)), & 179 2._wp*tr_slope(2:k_floor)/(c_no(2:k_floor)+1.e-10_wp)) ) 180 181 ! Compute limited flux: 182 flux_if(2:k_floor,jn) = w_if(2:k_floor,jn)* & 183 ( tr_c(2:k_floor) + & 184 0.5_wp*flux_lim(2:k_floor)*(1._wp-c_no(2:k_floor))* & 185 (tr_d(2:k_floor)-tr_c(2:k_floor)) ) 186 187 ! Compute pseudo update for trend aggregation: 188 tr_it(1:k_floor-1) = tr_it(1:k_floor-1) + & 189 z2dt/float(n_iter)/fse3t(ji,jj,1:k_floor-1)* & 190 flux_if(2:k_floor,jn) 191 tr_it(2:k_floor) = tr_it(2:k_floor) - & 192 z2dt/float(n_iter)/fse3t(ji,jj,2:k_floor)* & 193 flux_if(2:k_floor,jn) 194 195 ENDDO ! Iterative loop 196 197 ! Estimate rate of change from pseudo state updates (source 198 ! splitting): 199 tra(ji,jj,1:k_floor,jp_fabm_m1+jn) = & 200 tra(ji,jj,1:k_floor,jp_fabm_m1+jn) + & 201 (tr_it(1:k_floor) - trb(ji,jj,1:k_floor,jp_fabm_m1+jn))/z2dt 202 #if defined key_trdtrc && defined key_iomput 203 IF( lk_trdtrc .AND. ln_trdtrc( jp_fabm_m1+jn ) ) THEN 204 tr_vmv(ji,jj,1:k_floor,jn)=(tr_it(1:k_floor) - trb(ji,jj,1:k_floor,jn))/z2dt 90 IF (ALL(w_ct(ji,1:k_floor,jn) == 0._wp)) CYCLE 91 92 ! Compute velocities at interfaces 93 w_if(1:k_floor-1) = zwgt_if(1:k_floor-1) * w_ct(ji,1:k_floor-1,jn) + (1._wp - zwgt_if(1:k_floor-1)) * w_ct(ji,2:k_floor,jn) 94 95 ! Compute change (per volume) due to vertical movement per layer 96 IF (method == 1) THEN 97 CALL advect_1(k_floor, trn(ji,jj,1:k_floor,jp_fabm_m1+jn), w_if(1:k_floor-1), h(1:k_floor), z2dt, dc(1:k_floor)) 98 ELSE 99 CALL advect_3(k_floor, trb(ji,jj,1:k_floor,jp_fabm_m1+jn), w_if(1:k_floor-1), h(1:k_floor), z2dt, dc(1:k_floor)) 205 100 END IF 206 #endif 207 ENDDO ! State loop 101 102 ! Incorporate change due to vertical movement in sources-sinks 103 tra(ji,jj,1:k_floor,jp_fabm_m1+jn) = tra(ji,jj,1:k_floor,jp_fabm_m1+jn) + dc(1:k_floor) 104 105 #if defined key_trdtrc && defined key_iomput 106 ! Store change due to vertical movement as diagnostic 107 IF( lk_trdtrc .AND. ln_trdtrc( jp_fabm_m1+jn)) tr_vmv(ji,jj,1:k_floor,jn) = dc(1:k_floor) 108 #endif 109 END DO ! State loop 208 110 END IF ! Level check 209 111 END DO ! i-loop 210 112 END DO ! j-loop 113 211 114 #if defined key_trdtrc && defined key_iomput 212 115 DO jn=1,jp_fabm ! State loop … … 220 123 END SUBROUTINE compute_vertical_movement 221 124 125 SUBROUTINE advect_1(nk, c, w, h, dt, trend) 126 INTEGER, INTENT(IN) :: nk 127 REAL(wp), INTENT(IN) :: c(1:nk) 128 REAL(wp), INTENT(IN) :: w(1:nk-1) 129 REAL(wp), INTENT(IN) :: h(1:nk) 130 REAL(wp), INTENT(IN) :: dt 131 REAL(wp), INTENT(OUT) :: trend(1:nk) 132 133 REAL(wp) :: flux(0:nk) 134 INTEGER :: jk 135 ! Compute fluxes (per surface area) over at interfaces (remember: positive for upwards) 136 flux(0) = 0._wp 137 DO jk=1,nk-1 ! k-loop 138 IF (w(jk) > 0) THEN 139 ! Upward movement (source layer is jk+1) 140 flux(jk) = min(w(jk), h(jk+1)/dt) * c(jk+1) 141 ELSE 142 ! Downward movement (source layer is jk) 143 flux(jk) = max(w(jk), -h(jk)/dt) * c(jk) 144 END IF 145 END DO 146 flux(nk) = 0._wp 147 trend = (flux(1:nk) - flux(0:nk-1)) / h 148 END SUBROUTINE 149 150 SUBROUTINE advect_3(nk, c_old, w, h, dt, trend) 151 INTEGER, INTENT(IN) :: nk 152 REAL(wp), INTENT(IN) :: c_old(1:nk) 153 REAL(wp), INTENT(IN) :: w(1:nk-1) 154 REAL(wp), INTENT(IN) :: h(1:nk) 155 REAL(wp), INTENT(IN) :: dt 156 REAL(wp), INTENT(OUT) :: trend(1:nk) 157 158 INTEGER, PARAMETER :: n_itermax=100 159 REAL(wp) :: cmax_no 160 REAL(wp) :: cfl(1:nk-1) 161 INTEGER :: n_iter, n_count, jk 162 REAL(wp) :: c(1:nk) 163 REAL(wp) :: tr_u(1:nk-1) 164 REAL(wp) :: tr_c(1:nk-1) 165 REAL(wp) :: tr_d(1:nk-1) 166 REAL(wp) :: delta_tr_u(1:nk-1) 167 REAL(wp) :: delta_tr(1:nk-1) 168 REAL(wp) :: ratio(1:nk-1) 169 REAL(wp) :: x_fac(1:nk-1) 170 REAL(wp) :: phi_lim(1:nk-1) 171 REAL(wp) :: limiter(1:nk-1) 172 REAL(wp) :: flux_if(1:nk-1) 173 174 c(:) = c_old(:) 175 176 ! get maximum Courant number: 177 cfl = ABS(w) * dt / (0.5_wp * (h(2:nk) + h(1:nk-1))) 178 cmax_no = MAXVAL(cfl) 179 180 ! number of iterations: 181 n_iter = MIN(n_itermax, INT(cmax_no) + 1) 182 IF (ln_ctl.AND.(n_iter .gt. 1)) THEN 183 WRITE(numout,*) 'compute_vertical_movement::advect_3():' 184 WRITE(numout,*) ' Maximum Courant number is ',cmax_no,'.' 185 WRITE(numout,*) ' ',n_iter,' iterations used for vertical advection.' 186 ENDIF 187 188 ! effective Courant number: 189 cfl = cfl/n_iter 190 191 DO n_count=1,n_iter ! Iterative loop 192 ! Determine tracer concentration at 1.5 upstream (tr_u), 0.5 upstream (tr_c), 0.5 downstream (tr_d) from interface 193 IF (nk.gt.2) THEN 194 ! More than 2 vertical wet points 195 IF (nk.gt.3) THEN 196 WHERE (w(2:nk-2).ge.0._wp) 197 !upward movement 198 tr_u(2:nk-2)=c(4:nk) 199 tr_c(2:nk-2)=c(3:nk-1) 200 tr_d(2:nk-2)=c(2:nk-2) 201 ELSEWHERE 202 ! downward movement 203 tr_u(2:nk-2)=c(1:nk-3) 204 tr_c(2:nk-2)=c(2:nk-2) 205 tr_d(2:nk-2)=c(3:nk-1) 206 ENDWHERE 207 ENDIF 208 209 ! Interface between surface layer and the next 210 IF (w(1).ge.0._wp) THEN 211 ! upward movement 212 tr_u(1)=c(3) 213 tr_c(1)=c(2) 214 tr_d(1)=c(1) 215 ELSE 216 ! downward movement 217 tr_u(1)=c(1) 218 tr_c(1)=c(1) 219 tr_d(1)=c(2) 220 ENDIF 221 222 ! Interface between bottom layer and the previous 223 IF (w(nk-1).ge.0._wp) THEN 224 ! upward movement 225 tr_u(nk-1)=c(nk) 226 tr_c(nk-1)=c(nk) 227 tr_d(nk-1)=c(nk-1) 228 ELSE 229 ! downward movement 230 tr_u(nk-1)=c(nk-2) 231 tr_c(nk-1)=c(nk-1) 232 tr_d(nk-1)=c(nk) 233 ENDIF 234 ELSE 235 ! only 2 vertical wet points, i.e. only 1 interface 236 IF (w(1).ge.0._wp) THEN 237 ! upward movement 238 tr_u(1)=c(2) 239 tr_c(1)=c(2) 240 tr_d(1)=c(1) 241 ELSE 242 ! downward movement 243 tr_u(1)=c(1) 244 tr_c(1)=c(1) 245 tr_d(1)=c(2) 246 ENDIF 247 ENDIF 248 249 delta_tr_u = tr_c - tr_u 250 delta_tr = tr_d - tr_c 251 WHERE (delta_tr * delta_tr_u > 0._wp) 252 ! Monotonic function over tr_u, tr_c, r_d 253 254 ! Compute slope ratio 255 ratio = delta_tr_u / delta_tr 256 257 ! QUICKEST flux limiter 258 x_fac = (1._wp - 2._wp * cfl) / 6._wp 259 phi_lim = (0.5_wp + x_fac) + (0.5_wp - x_fac) * ratio 260 limiter = MIN(phi_lim, 2._wp / (1._wp - cfl), 2._wp * ratio / (cfl + 1.e-10_wp)) 261 262 ! Compute limited flux 263 flux_if = w * (tr_c + 0.5_wp * limiter * (1._wp - cfl) * delta_tr) 264 ELSEWHERE 265 ! Non-monotonic, use 1st order upstream 266 flux_if = w * tr_c 267 ENDWHERE 268 269 ! Compute pseudo update for trend aggregation: 270 c(1:nk-1) = c(1:nk-1) + dt / real(n_iter, wp) / h(1:nk-1) * flux_if 271 c(2:nk) = c(2:nk) - dt / real(n_iter, wp) / h(2:nk) * flux_if 272 273 ENDDO ! Iterative loop 274 275 ! Estimate rate of change from pseudo state updates (source splitting): 276 trend = (c - c_old) / dt 277 END SUBROUTINE 278 222 279 #endif 223 280 END MODULE
Note: See TracChangeset
for help on using the changeset viewer.