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.
Changeset 13489 for branches/UKMO/AMM15_v3_6_STABLE_package_collate_utils366_fabmv1_v2/NEMOGCM/NEMO/TOP_SRC/FABM/vertical_movement_fabm.F90 – NEMO

Ignore:
Timestamp:
2020-09-17T19:13:57+02:00 (4 years ago)
Author:
dford
Message:

Merge in previous branch.

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  
    1515   USE trc 
    1616   USE par_fabm 
    17    USE fabm 
    1817   USE dom_oce 
    1918#if defined key_trdtrc && defined key_iomput 
     
    2524 
    2625#  include "domzgr_substitute.h90" 
     26#  include "vectopt_loop_substitute.h90" 
    2727 
    2828   PRIVATE 
     
    3232   ! Work arrays for vertical advection (residual movement/sinking/floating) 
    3333   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 
    3734#if defined key_trdtrc && defined key_iomput 
    3835   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:,:) :: tr_vmv 
     
    4138   CONTAINS 
    4239 
    43    SUBROUTINE compute_vertical_movement( kt ) 
     40   SUBROUTINE compute_vertical_movement( kt, method ) 
    4441      !!---------------------------------------------------------------------- 
    4542      !!                     ***  compute_vertical_movement  *** 
    4643      !! 
    47       !! ** Purpose :   compute vertical movement of FABM tracers 
     44      !! ** Purpose : compute vertical movement of FABM tracers through the water 
     45      !!              (sinking/floating/active movement) 
    4846      !! 
    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. 
    5449      !!---------------------------------------------------------------------- 
    5550      ! 
    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) 
    6256#if defined key_trdtrc 
    6357      CHARACTER (len=20) :: cltra 
     
    6963 
    7064      IF( neuler == 0 .AND. kt == nittrc000 ) THEN 
    71           z2dt = rdt                  ! set time step size (Euler) 
     65         z2dt = rdt                  ! set time step size (Euler) 
    7266      ELSE 
    73           z2dt = 2._wp * rdt          ! set time step size (Leapfrog) 
     67         z2dt = 2._wp * rdt          ! set time step size (Leapfrog) 
    7468      ENDIF 
     69 
    7570      ! 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,:)) 
     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,:)) 
    8075         END DO 
    81  
    82          DO ji=1,jpi ! i-loop 
     76         DO ji=fs_2,fs_jpim1 ! i-loop 
    8377            ! 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) 
     78            k_floor = mbkt(ji,jj) 
     79            IF (k_floor > 1) THEN 
    8680               ! Linearly interpolate to velocities at the interfaces between layers 
    8781               ! Note: 
    88                !    - interface k sits between cell centre k and k-1, 
    89                !    - k [1,jpk] increases downwards 
     82               !    - interface k sits between cell centre k and k+1 (k=0 for surface) 
     83               !    - k [1,jpkm1] increases downwards 
    9084               !    - 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)) 
    10087 
    10188               ! Advect: 
    10289               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)) 
    205100                  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 
    208110            END IF ! Level check 
    209111         END DO ! i-loop 
    210112      END DO ! j-loop 
     113 
    211114#if defined key_trdtrc && defined key_iomput 
    212115      DO jn=1,jp_fabm ! State loop 
     
    220123   END SUBROUTINE compute_vertical_movement 
    221124 
     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 
    222279#endif 
    223280END MODULE 
Note: See TracChangeset for help on using the changeset viewer.