New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
vertical_movement_fabm.F90 in branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils305/NEMOGCM/NEMO/TOP_SRC/FABM – NEMO

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils305/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90 @ 11979

Last change on this file since 11979 was 11979, checked in by dford, 4 years ago

Bug fixes: update indexing, and initialise arrays after allocation.
File size: 9.8 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 fabm
18   USE dom_oce
19#if defined key_trdtrc && defined key_iomput
20   USE iom
21   USE trdtrc_oce
22#endif
23
24   IMPLICIT NONE
25
26#  include "domzgr_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   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:)   :: w_if
35   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:)   :: zwgt_if
36   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:)   :: flux_if
37#if defined key_trdtrc && defined key_iomput
38   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:,:) :: tr_vmv
39#endif
40
41   CONTAINS
42
43   SUBROUTINE compute_vertical_movement( kt )
44      !!----------------------------------------------------------------------
45      !!                     ***  compute_vertical_movement  ***
46      !!
47      !! ** Purpose :   compute vertical movement of FABM tracers
48      !!
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.
54      !!----------------------------------------------------------------------
55      !
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
62#if defined key_trdtrc
63      CHARACTER (len=20) :: cltra
64#endif
65
66#if defined key_trdtrc && defined key_iomput
67      IF( lk_trdtrc ) tr_vmv = 0.0_wp
68#endif
69
70      IF( neuler == 0 .AND. kt == nittrc000 ) THEN
71          z2dt = rdt                  ! set time step size (Euler)
72      ELSE
73          z2dt = 2._wp * rdt          ! set time step size (Leapfrog)
74      ENDIF
75      ! Compute interior vertical velocities and include them in source array.
76      DO jj=1,jpj ! j-loop
77         ! 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,:))
80         END DO
81
82         DO ji=1,jpi ! i-loop
83            ! Only process this horizontal point (ji,jj) if number of layers exceeds 1
84            IF (mbkt(ji,jj)>1) THEN ! Level check
85               k_floor=mbkt(ji,jj)
86               ! Linearly interpolate to velocities at the interfaces between layers
87               ! Note:
88               !    - interface k sits between cell centre k and k-1,
89               !    - k [1,jpk] increases downwards
90               !    - 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
100
101               ! Advect:
102               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
205                  END IF
206#endif
207               ENDDO ! State loop
208            END IF ! Level check
209         END DO ! i-loop
210      END DO ! j-loop
211#if defined key_trdtrc && defined key_iomput
212      DO jn=1,jp_fabm ! State loop
213        IF( lk_trdtrc .AND. ln_trdtrc(jp_fabm_m1+jn) ) THEN
214          cltra = 'VMV_'//TRIM(ctrcnm(jp_fabm_m1+jn))
215          CALL iom_put( cltra,  tr_vmv(:,:,:,jn) )
216        END IF
217      ENDDO
218#endif
219
220   END SUBROUTINE compute_vertical_movement
221
222#endif
223END MODULE
Note: See TracBrowser for help on using the repository browser.