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.
divhor.F90 in NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN – NEMO

source: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/divhor.F90 @ 15620

Last change on this file since 15620 was 15620, checked in by techene, 2 years ago

#2605 update RK3 loops according to the trunk transformation
  • Property svn:keywords set to Id
File size: 9.0 KB
Line 
1MODULE divhor
2   !!==============================================================================
3   !!                       ***  MODULE  divhor  ***
4   !! Ocean diagnostic variable : now horizontal divergence
5   !!==============================================================================
6   !! History :  1.0  ! 2002-09  (G. Madec, E. Durand)  Free form, F90
7   !!             -   ! 2005-01  (J. Chanut) Unstructured open boundaries
8   !!             -   ! 2003-08  (G. Madec)  merged of cur and div, free form, F90
9   !!             -   ! 2005-01  (J. Chanut, A. Sellar) unstructured open boundaries
10   !!            3.3  ! 2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
11   !!             -   ! 2010-10  (R. Furner, G. Madec) runoff and cla added directly here
12   !!            3.7  ! 2014-01  (G. Madec) suppression of velocity curl from in-core memory
13   !!             -   ! 2014-12  (G. Madec) suppression of cross land advection option
14   !!             -   ! 2015-10  (G. Madec) add velocity and rnf flag in argument of div_hor
15   !!            4.5  ! 2015-10  (S. Techene, G. Madec)  hdiv replaced by e3divUh
16   !!----------------------------------------------------------------------
17
18   !!----------------------------------------------------------------------
19   !!   div_hor    : Compute the horizontal divergence field
20   !!----------------------------------------------------------------------
21   USE oce             ! ocean dynamics and tracers
22   USE dom_oce         ! ocean space and time domain
23   USE sbc_oce, ONLY : ln_rnf      ! river runoff
24   USE sbcrnf , ONLY : sbc_rnf_div ! river runoff
25   USE isf_oce, ONLY : ln_isf      ! ice shelf
26   USE isfhdiv, ONLY : isf_hdiv    ! ice shelf
27#if defined key_asminc   
28   USE asminc          ! Assimilation increment
29#endif
30   !
31   USE in_out_manager  ! I/O manager
32   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
33   USE lib_mpp         ! MPP library
34   USE timing          ! Timing
35
36   IMPLICIT NONE
37   PRIVATE
38
39   !                  !! * Interface
40   INTERFACE div_hor
41      MODULE PROCEDURE div_hor_RK3, div_hor_old
42   END INTERFACE
43
44   PUBLIC   div_hor    ! routine called by ssh_nxt.F90 and istate.F90
45
46   !! * Substitutions
47#  include "do_loop_substitute.h90"
48#  include "domzgr_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
51   !! $Id$
52   !! Software governed by the CeCILL license (see ./LICENSE)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE div_hor_RK3( kt, Kbb, Kmm, puu, pvv, pe3divUh )
57      !!----------------------------------------------------------------------
58      !!                  ***  ROUTINE div_hor_RK3  ***
59      !!                   
60      !! ** Purpose :   compute the horizontal divergence at now time-step
61      !!
62      !! ** Method  :   the now divergence is computed as :
63      !!         hdiv = 1/(e1e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
64      !!      and correct with runoff inflow (div_rnf) and cross land flow (div_cla)
65      !!
66      !! ** Action  : - thickness weighted horizontal divergence of in input velocity (puu,pvv)
67      !!----------------------------------------------------------------------
68      INTEGER                         , INTENT(in   ) ::   kt, Kbb, Kmm   ! ocean time-step & time-level indices
69      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   puu, pvv       ! horizontal velocity
70      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   pe3divUh       ! e3t*div[Uh]
71      !
72      INTEGER  ::   ji, jj, jk    ! dummy loop indices
73      !!----------------------------------------------------------------------
74      !
75      IF( ln_timing )   CALL timing_start('div_hor_RK3')
76      !
77      IF( kt == nit000 ) THEN
78         IF(lwp) WRITE(numout,*)
79         IF(lwp) WRITE(numout,*) 'div_hor_RK3 : thickness weighted horizontal divergence '
80         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
81         hdiv    (:,:,:) = 0._wp    ! initialize hdiv & pe3divUh for the halos and jpk level at the first time step
82      ENDIF
83      !
84      pe3divUh(:,:,:) = 0._wp    !!gm to be applied to the halos only
85      !
86      DO_3D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 1, jpkm1 )
87         hdiv(ji,jj,jk) = (   e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * puu(ji  ,jj,jk)      &
88            &               - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * puu(ji-1,jj,jk)      &
89            &               + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * pvv(ji,jj  ,jk)      &
90            &               - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * pvv(ji,jj-1,jk)  )   &
91            &            * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
92      END_3D
93      !
94      IF( ln_rnf )   CALL sbc_rnf_div( hdiv, Kmm )             !==  + runoffs divergence  ==!
95      !
96#if defined key_asminc 
97      IF( ln_sshinc .AND. ln_asmiau )   &                      !==  + SSH assimilation increment  ==!
98         &           CALL ssh_asm_div( kt, Kbb, Kmm, hdiv )
99#endif
100      !
101      IF( ln_isf )   CALL isf_hdiv( kt, Kmm, hdiv )            !==  + ice-shelf mass exchange ==!
102      !
103      IF( nn_hls==1 )   CALL lbc_lnk( 'divhor', hdiv, 'T', 1._wp )   !   (no sign change)
104      !
105!!gm Patch before suppression of hdiv from all modules that use it
106!      DO_3D( 0, 0, 0, 0, 1, jpkm1 )                            !==  e3t * Horizontal divergence  ==!
107!         pe3divUh(ji,jj,jk) = hdiv(ji,jj,jk) * e3t(ji,jj,jk,Kmm)
108!      END_3D
109!JC: over whole domain, and after lbclnk on hdiv to prevent from reproducibility issues
110      DO_3D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 1, jpkm1 )
111         pe3divUh(ji,jj,jk) = hdiv(ji,jj,jk) * e3t(ji,jj,jk,Kmm)
112      END_3D
113!!gm end
114      !
115      !
116      IF( ln_timing )   CALL timing_stop('div_hor_RK3')
117      !
118   END SUBROUTINE div_hor_RK3
119
120
121   SUBROUTINE div_hor_old( kt, Kbb, Kmm )
122      !!----------------------------------------------------------------------
123      !!                  ***  ROUTINE div_hor  ***
124      !!                   
125      !! ** Purpose :   compute the horizontal divergence at now time-step
126      !!
127      !! ** Method  :   the now divergence is computed as :
128      !!         hdiv = 1/(e1e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
129      !!      and correct with runoff inflow (div_rnf) and cross land flow (div_cla)
130      !!
131      !! ** Action  : - update hdiv, the now horizontal divergence
132      !!----------------------------------------------------------------------
133      INTEGER, INTENT(in) ::   kt        ! ocean time-step index
134      INTEGER, INTENT(in) ::   Kbb, Kmm  ! ocean time level indices
135      !
136      INTEGER  ::   ji, jj, jk    ! dummy loop indices
137      !!----------------------------------------------------------------------
138      !
139      IF( ln_timing )   CALL timing_start('div_hor')
140      !
141      IF( kt == nit000 ) THEN
142         IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
143            IF(lwp) WRITE(numout,*)
144            IF(lwp) WRITE(numout,*) 'div_hor : horizontal velocity divergence '
145            IF(lwp) WRITE(numout,*) '~~~~~~~   '
146         ENDIF
147         DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
148            hdiv(ji,jj,jk) = 0._wp    ! initialize hdiv for the halos at the first time step
149         END_3D
150      ENDIF
151      !
152      DO_3D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 1, jpkm1 )                                          !==  Horizontal divergence  ==!
153         ! round brackets added to fix the order of floating point operations
154         ! needed to ensure halo 1 - halo 2 compatibility
155         hdiv(ji,jj,jk) = (  ( e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * uu(ji  ,jj,jk,Kmm)     &
156            &                - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm)     &
157            &                )                                                             & ! bracket for halo 1 - halo 2 compatibility
158            &              + ( e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * vv(ji,jj  ,jk,Kmm)     &
159            &                - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm)     &
160            &                )                                                             & ! bracket for halo 1 - halo 2 compatibility
161            &             )  * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
162      END_3D
163      !
164      IF( ln_rnf )   CALL sbc_rnf_div( hdiv, Kmm )                               !==  runoffs    ==!   (update hdiv field)
165      !
166#if defined key_asminc 
167      IF( ln_sshinc .AND. ln_asmiau )   CALL ssh_asm_div( kt, Kbb, Kmm, hdiv )   !==  SSH assimilation  ==!   (update hdiv field)
168      !
169#endif
170      IF( ln_isf )                      CALL isf_hdiv( kt, Kmm, hdiv )           !==  ice shelf         ==!   (update hdiv field)
171      !
172      IF( nn_hls==1 )   CALL lbc_lnk( 'divhor', hdiv, 'T', 1.0_wp )   !   (no sign change)
173      !                                                               ! needed for ww in sshwzv
174      IF( ln_timing )   CALL timing_stop('div_hor')
175      !
176   END SUBROUTINE div_hor_old
177   
178   !!======================================================================
179END MODULE divhor
Note: See TracBrowser for help on using the repository browser.