source: branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90

Last change on this file was 13576, checked in by dford, 2 weeks ago

Update NEMO-FABM coupler for FABM v1, and introduce two-way NEMO-ERSEM coupling options. See https://code.metoffice.gov.uk/trac/utils/ticket/366.

File size: 10.1 KB
Line 
1MODULE vertical_movement_fabm
2   !!======================================================================
3   !!                         ***  MODULE vertical_movement_fabm  ***
4   !! TOP :   Module for the vertical movement of the FABM tracers
5   !!======================================================================
6
7#if defined key_fabm
8   !!----------------------------------------------------------------------
9   !!   'key_fabm'                                               FABM tracers
10   !!----------------------------------------------------------------------
11   !! compute_vertical_movement : compute vertical movement of FABM fields
12   !!----------------------------------------------------------------------
13   USE par_trc
14   USE oce_trc
15   USE trc
16   USE par_fabm
17   USE dom_oce
18#if defined key_trdtrc && defined key_iomput
19   USE iom
20   USE trdtrc_oce
21#endif
22
23   IMPLICIT NONE
24
25#  include "domzgr_substitute.h90"
26#  include "vectopt_loop_substitute.h90"
27
28   PRIVATE
29
30   PUBLIC compute_vertical_movement
31
32   ! Work arrays for vertical advection (residual movement/sinking/floating)
33   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:) :: w_ct
34#if defined key_trdtrc && defined key_iomput
35   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:,:) :: tr_vmv
36#endif
37
38   CONTAINS
39
40   SUBROUTINE compute_vertical_movement( kt, method )
41      !!----------------------------------------------------------------------
42      !!                     ***  compute_vertical_movement  ***
43      !!
44      !! ** Purpose : compute vertical movement of FABM tracers through the water
45      !!              (sinking/floating/active movement)
46      !!
47      !! ** Method  : Retrieves additional vertical velocity field and applies
48      !!              advection scheme.
49      !!----------------------------------------------------------------------
50      !
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)
56#if defined key_trdtrc
57      CHARACTER (len=20) :: cltra
58#endif
59
60#if defined key_trdtrc && defined key_iomput
61      IF( lk_trdtrc ) tr_vmv = 0.0_wp
62#endif
63
64      IF( neuler == 0 .AND. kt == nittrc000 ) THEN
65         z2dt = rdt                  ! set time step size (Euler)
66      ELSE
67         z2dt = 2._wp * rdt          ! set time step size (Leapfrog)
68      ENDIF
69
70      ! Compute interior vertical velocities and include them in source array.
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,:))
75         END DO
76         DO ji=fs_2,fs_jpim1 ! i-loop
77            ! Only process this horizontal point (ji,jj) if number of layers exceeds 1
78            k_floor = mbkt(ji,jj)
79            IF (k_floor > 1) THEN
80               ! Linearly interpolate to velocities at the interfaces between layers
81               ! Note:
82               !    - interface k sits between cell centre k and k+1 (k=0 for surface)
83               !    - k [1,jpkm1] increases downwards
84               !    - upward velocity is positive, downward velocity is negative
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))
87
88               ! Advect:
89               DO jn=1,jp_fabm ! State loop
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))
100                  END IF
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
110            END IF ! Level check
111         END DO ! i-loop
112      END DO ! j-loop
113
114#if defined key_trdtrc && defined key_iomput
115      DO jn=1,jp_fabm ! State loop
116        IF( lk_trdtrc .AND. ln_trdtrc(jp_fabm_m1+jn) ) THEN
117          cltra = 'VMV_'//TRIM(ctrcnm(jp_fabm_m1+jn))
118          CALL iom_put( cltra,  tr_vmv(:,:,:,jn) )
119        END IF
120      ENDDO
121#endif
122
123   END SUBROUTINE compute_vertical_movement
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
279#endif
280END MODULE
Note: See TracBrowser for help on using the repository browser.