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.
sshwzv.F90 in NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/DYN – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/DYN/sshwzv.F90 @ 11463

Last change on this file since 11463 was 11426, checked in by davestorkey, 5 years ago

branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap : Merge in changes from main IMMERSE branch.
  • Property svn:keywords set to Id
File size: 17.8 KB
RevLine 
[1565]1MODULE sshwzv   
[3]2   !!==============================================================================
[1438]3   !!                       ***  MODULE  sshwzv  ***
4   !! Ocean dynamics : sea surface height and vertical velocity
[3]5   !!==============================================================================
[1438]6   !! History :  3.1  !  2009-02  (G. Madec, M. Leclair)  Original code
[2528]7   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA
8   !!             -   !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
9   !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
[4292]10   !!            3.3  !  2011-10  (M. Leclair) split former ssh_wzv routine and remove all vvl related work
[3]11   !!----------------------------------------------------------------------
[1438]12
[3]13   !!----------------------------------------------------------------------
[6140]14   !!   ssh_nxt       : after ssh
[11057]15   !!   ssh_atf       : filter ans swap the ssh arrays
[6140]16   !!   wzv           : compute now vertical velocity
[1438]17   !!----------------------------------------------------------------------
[6140]18   USE oce            ! ocean dynamics and tracers variables
19   USE dom_oce        ! ocean space and time domain variables
20   USE sbc_oce        ! surface boundary condition: ocean
21   USE domvvl         ! Variable volume
22   USE divhor         ! horizontal divergence
23   USE phycst         ! physical constants
[9019]24   USE bdy_oce , ONLY : ln_bdy, bdytmask   ! Open BounDarY
[6140]25   USE bdydyn2d       ! bdy_ssh routine
[2528]26#if defined key_agrif
[9570]27   USE agrif_oce_interp
[2528]28#endif
[6140]29   !
[10364]30   USE iom 
[6140]31   USE in_out_manager ! I/O manager
32   USE restart        ! only for lrst_oce
33   USE prtctl         ! Print control
34   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
35   USE lib_mpp        ! MPP library
36   USE timing         ! Timing
[9023]37   USE wet_dry        ! Wetting/Drying flux limiting
[592]38
[3]39   IMPLICIT NONE
40   PRIVATE
41
[1438]42   PUBLIC   ssh_nxt    ! called by step.F90
[4292]43   PUBLIC   wzv        ! called by step.F90
[10364]44   PUBLIC   wAimp      ! called by step.F90
[11057]45   PUBLIC   ssh_atf    ! called by step.F90
[3]46
47   !! * Substitutions
[1438]48#  include "vectopt_loop_substitute.h90"
[3]49   !!----------------------------------------------------------------------
[9598]50   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[888]51   !! $Id$
[10068]52   !! Software governed by the CeCILL license (see ./LICENSE)
[592]53   !!----------------------------------------------------------------------
[3]54CONTAINS
55
[10978]56   SUBROUTINE ssh_nxt( kt, Kbb, Kmm, pssh, Kaa )
[3]57      !!----------------------------------------------------------------------
[4292]58      !!                ***  ROUTINE ssh_nxt  ***
[1438]59      !!                   
[10978]60      !! ** Purpose :   compute the after ssh (ssh(Kaa))
[3]61      !!
[4292]62      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment
63      !!      is computed by integrating the horizontal divergence and multiply by
64      !!      by the time step.
[3]65      !!
[10978]66      !! ** action  :   ssh(:,:,Kaa), after sea surface height
[2528]67      !!
68      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
[3]69      !!----------------------------------------------------------------------
[10978]70      INTEGER                         , INTENT(in   ) ::   kt             ! time step
71      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! time level index
72      REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! sea-surface height
[4292]73      !
[7753]74      INTEGER  ::   jk            ! dummy loop indice
[5836]75      REAL(wp) ::   z2dt, zcoef   ! local scalars
[9019]76      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv   ! 2D workspace
[3]77      !!----------------------------------------------------------------------
[3294]78      !
[9019]79      IF( ln_timing )   CALL timing_start('ssh_nxt')
[3294]80      !
[3]81      IF( kt == nit000 ) THEN
82         IF(lwp) WRITE(numout,*)
[4292]83         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
[1438]84         IF(lwp) WRITE(numout,*) '~~~~~~~ '
[3]85      ENDIF
[2528]86      !
[7646]87      z2dt = 2._wp * rdt                          ! set time step size (Euler/Leapfrog)
[2715]88      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
[7646]89      zcoef = 0.5_wp * r1_rau0
[3]90
[1438]91      !                                           !------------------------------!
92      !                                           !   After Sea Surface Height   !
93      !                                           !------------------------------!
[9023]94      IF(ln_wd_il) THEN
[11053]95         CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), z2dt, Kmm, uu, vv )
[7753]96      ENDIF
[7646]97
[10969]98      CALL div_hor( kt, Kbb, Kmm )                     ! Horizontal divergence
[7646]99      !
[7753]100      zhdiv(:,:) = 0._wp
[1438]101      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
[10978]102        zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
[1438]103      END DO
104      !                                                ! Sea surface elevation time stepping
[4338]105      ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
106      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
107      !
[10978]108      pssh(:,:,Kaa) = (  pssh(:,:,Kbb) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
[9023]109      !
110#if defined key_agrif
[11426]111      Kbb_a = Kbb; Kmm_a = Kmm; Krhs_a = Kaa; CALL agrif_ssh( kt )
[9023]112#endif
113      !
[5930]114      IF ( .NOT.ln_dynspg_ts ) THEN
[7646]115         IF( ln_bdy ) THEN
[10978]116            CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1. )    ! Not sure that's necessary
117            CALL bdy_ssh( pssh(:,:,Kaa) )             ! Duplicate sea level across open boundaries
[5930]118         ENDIF
[4292]119      ENDIF
120      !                                           !------------------------------!
121      !                                           !           outputs            !
122      !                                           !------------------------------!
123      !
[10978]124      IF(ln_ctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa)  - : ', mask1=tmask )
[4292]125      !
[9019]126      IF( ln_timing )   CALL timing_stop('ssh_nxt')
[4292]127      !
128   END SUBROUTINE ssh_nxt
129
130   
[10978]131   SUBROUTINE wzv( kt, Kbb, Kmm, pww, Kaa )
[4292]132      !!----------------------------------------------------------------------
133      !!                ***  ROUTINE wzv  ***
134      !!                   
135      !! ** Purpose :   compute the now vertical velocity
136      !!
137      !! ** Method  : - Using the incompressibility hypothesis, the vertical
138      !!      velocity is computed by integrating the horizontal divergence 
139      !!      from the bottom to the surface minus the scale factor evolution.
140      !!        The boundary conditions are w=0 at the bottom (no flux) and.
141      !!
[10978]142      !! ** action  :   pww      : now vertical velocity
[4292]143      !!
144      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
145      !!----------------------------------------------------------------------
[10978]146      INTEGER                         , INTENT(in)    ::   kt             ! time step
147      INTEGER                         , INTENT(in)    ::   Kbb, Kmm, Kaa  ! time level indices
148      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pww            ! now vertical velocity
[4292]149      !
[5836]150      INTEGER  ::   ji, jj, jk   ! dummy loop indices
151      REAL(wp) ::   z1_2dt       ! local scalars
[9019]152      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zhdiv
[4292]153      !!----------------------------------------------------------------------
154      !
[9019]155      IF( ln_timing )   CALL timing_start('wzv')
[5836]156      !
[4292]157      IF( kt == nit000 ) THEN
158         IF(lwp) WRITE(numout,*)
159         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
160         IF(lwp) WRITE(numout,*) '~~~~~ '
161         !
[10978]162         pww(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
[4292]163      ENDIF
164      !                                           !------------------------------!
165      !                                           !     Now Vertical Velocity    !
166      !                                           !------------------------------!
167      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog)
168      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
169      !
170      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
[9019]171         ALLOCATE( zhdiv(jpi,jpj,jpk) ) 
[4292]172         !
173         DO jk = 1, jpkm1
174            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
[4338]175            ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
[4292]176            DO jj = 2, jpjm1
177               DO ji = fs_2, fs_jpim1   ! vector opt.
[5836]178                  zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
[4292]179               END DO
[592]180            END DO
181         END DO
[10425]182         CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
[4292]183         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
[10978]184         !                             ! Same question holds for hdiv. Perhaps just for security
[4292]185         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
186            ! computation of w
[10978]187            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk) + zhdiv(:,:,jk)    &
188               &                         + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) )     ) * tmask(:,:,jk)
[4292]189         END DO
[10978]190         !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
[9019]191         DEALLOCATE( zhdiv ) 
[4292]192      ELSE   ! z_star and linear free surface cases
193         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
[7753]194            ! computation of w
[10978]195            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)                 &
196               &                         + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) )  ) * tmask(:,:,jk)
[4292]197         END DO
[1438]198      ENDIF
[592]199
[7646]200      IF( ln_bdy ) THEN
[4327]201         DO jk = 1, jpkm1
[10978]202            pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:)
[4327]203         END DO
204      ENDIF
[4292]205      !
[9023]206#if defined key_agrif 
207      IF( .NOT. AGRIF_Root() ) THEN
[10978]208         IF ((nbondi ==  1).OR.(nbondi == 2)) pww(nlci-1 , :     ,:) = 0.e0      ! east
209         IF ((nbondi == -1).OR.(nbondi == 2)) pww(2      , :     ,:) = 0.e0      ! west
210         IF ((nbondj ==  1).OR.(nbondj == 2)) pww(:      ,nlcj-1 ,:) = 0.e0      ! north
211         IF ((nbondj == -1).OR.(nbondj == 2)) pww(:      ,2      ,:) = 0.e0      ! south
[9023]212      ENDIF 
213#endif 
[5836]214      !
[9124]215      IF( ln_timing )   CALL timing_stop('wzv')
[9023]216      !
[5836]217   END SUBROUTINE wzv
[592]218
219
[11057]220   SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh )
[1438]221      !!----------------------------------------------------------------------
[11057]222      !!                    ***  ROUTINE ssh_atf  ***
[1438]223      !!
[11057]224      !!      !!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!
225      !!
[1438]226      !! ** Purpose :   achieve the sea surface  height time stepping by
227      !!              applying Asselin time filter and swapping the arrays
[10978]228      !!              ssh(:,:,Kaa)  already computed in ssh_nxt 
[1438]229      !!
[2528]230      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
231      !!              from the filter, see Leclair and Madec 2010) and swap :
[10978]232      !!                ssh(:,:,Kmm) = ssh(:,:,Kaa) + atfp * ( ssh(:,:,Kbb) -2 ssh(:,:,Kmm) + ssh(:,:,Kaa) )
[2528]233      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
[10978]234      !!                ssh(:,:,Kmm) = ssh(:,:,Kaa)
[1438]235      !!
[10978]236      !! ** action  : - ssh(:,:,Kbb), ssh(:,:,Kmm)   : before & now sea surface height
[1438]237      !!                               ready for the next time step
[2528]238      !!
239      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
[1438]240      !!----------------------------------------------------------------------
[11057]241      INTEGER                         , INTENT(in   ) ::   kt             ! ocean time-step index
242      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! ocean time level indices
243      REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! SSH field
[6140]244      !
245      REAL(wp) ::   zcoef   ! local scalar
[1438]246      !!----------------------------------------------------------------------
[3294]247      !
[11057]248      IF( ln_timing )   CALL timing_start('ssh_atf')
[3294]249      !
[1438]250      IF( kt == nit000 ) THEN
251         IF(lwp) WRITE(numout,*)
[11057]252         IF(lwp) WRITE(numout,*) 'ssh_atf : Asselin time filter of sea surface height'
[1438]253         IF(lwp) WRITE(numout,*) '~~~~~~~ '
254      ENDIF
[6140]255      !              !==  Euler time-stepping: no filter, just swap  ==!
[11057]256      IF ( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN   ! Only do time filtering for leapfrog timesteps
257         !                                                  ! filtered "now" field
258         ssh(:,:,Kmm) = ssh(:,:,Kmm) + atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa) )
259         IF( .NOT.ln_linssh ) THEN                          ! "now" <-- with forcing removed
[6140]260            zcoef = atfp * rdt * r1_rau0
[11057]261            ssh(:,:,Kmm) = ssh(:,:,Kmm) - zcoef * (     emp_b(:,:) - emp   (:,:)   &
[7753]262               &                             -    rnf_b(:,:) + rnf   (:,:)   &
263               &                             + fwfisf_b(:,:) - fwfisf(:,:)   ) * ssmask(:,:)
[6140]264         ENDIF
[1438]265      ENDIF
266      !
[11057]267      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssh(:,:,Kmm), clinfo1=' ssh(:,:,Kmm)  - : ', mask1=tmask )
[2528]268      !
[11057]269      IF( ln_timing )   CALL timing_stop('ssh_atf')
[3294]270      !
[11057]271   END SUBROUTINE ssh_atf
[3]272
[10978]273   SUBROUTINE wAimp( kt, Kmm )
[10364]274      !!----------------------------------------------------------------------
275      !!                ***  ROUTINE wAimp  ***
276      !!                   
277      !! ** Purpose :   compute the Courant number and partition vertical velocity
278      !!                if a proportion needs to be treated implicitly
279      !!
280      !! ** Method  : -
281      !!
[10978]282      !! ** action  :   ww      : now vertical velocity (to be handled explicitly)
[10364]283      !!            :   wi      : now vertical velocity (for implicit treatment)
284      !!
285      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
286      !!----------------------------------------------------------------------
287      INTEGER, INTENT(in) ::   kt   ! time step
[10978]288      INTEGER, INTENT(in) ::   Kmm  ! time level index
[10364]289      !
290      INTEGER  ::   ji, jj, jk   ! dummy loop indices
291      REAL(wp)             ::   zCu, zcff, z1_e3w                     ! local scalars
292      REAL(wp) , PARAMETER ::   Cu_min = 0.15_wp                      ! local parameters
293      REAL(wp) , PARAMETER ::   Cu_max = 0.27                         ! local parameters
294      REAL(wp) , PARAMETER ::   Cu_cut = 2._wp*Cu_max - Cu_min        ! local parameters
295      REAL(wp) , PARAMETER ::   Fcu    = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters
296      !!----------------------------------------------------------------------
297      !
298      IF( ln_timing )   CALL timing_start('wAimp')
299      !
300      IF( kt == nit000 ) THEN
301         IF(lwp) WRITE(numout,*)
302         IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
303         IF(lwp) WRITE(numout,*) '~~~~~ '
304         !
305         Cu_adv(:,:,jpk) = 0._wp              ! bottom value : Cu_adv=0 (set once for all)
306      ENDIF
307      !
308      DO jk = 1, jpkm1            ! calculate Courant numbers
309         DO jj = 2, jpjm1
310            DO ji = 2, fs_jpim1   ! vector opt.
[10978]311               z1_e3w = 1._wp / e3w(ji,jj,jk,Kmm)
312               Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )    &
313                  &                      + ( MAX( e2u(ji  ,jj)*e3uw(ji  ,jj,jk,Kmm)*uu(ji  ,jj,jk,Kmm), 0._wp ) -   &
314                  &                          MIN( e2u(ji-1,jj)*e3uw(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) )   &
[10364]315                  &                        * r1_e1e2t(ji,jj)                                                  &
[10978]316                  &                      + ( MAX( e1v(ji,jj  )*e3vw(ji,jj  ,jk,Kmm)*vv(ji,jj  ,jk,Kmm), 0._wp ) -   &
317                  &                          MIN( e1v(ji,jj-1)*e3vw(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm), 0._wp ) )   &
[10364]318                  &                        * r1_e1e2t(ji,jj)                                                  &
319                  &                      ) * z1_e3w
320            END DO
321         END DO
322      END DO
323      !
324      CALL iom_put("Courant",Cu_adv)
325      !
326      wi(:,:,:) = 0._wp                                 ! Includes top and bottom values set to zero
327      IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere
328         DO jk = 1, jpkm1                               ! or scan Courant criterion and partition
329            DO jj = 2, jpjm1                            ! w where necessary
330               DO ji = 2, fs_jpim1   ! vector opt.
331                  !
332                  zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) )
333                  !
334                  IF( zCu < Cu_min ) THEN               !<-- Fully explicit
335                     zcff = 0._wp
336                  ELSEIF( zCu < Cu_cut ) THEN           !<-- Mixed explicit
337                     zcff = ( zCu - Cu_min )**2
338                     zcff = zcff / ( Fcu + zcff )
339                  ELSE                                  !<-- Mostly implicit
340                     zcff = ( zCu - Cu_max )/ zCu
341                  ENDIF
342                  zcff = MIN(1._wp, zcff)
343                  !
[10978]344                  wi(ji,jj,jk) =           zcff   * ww(ji,jj,jk)
345                  ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk)
[10364]346                  !
347                  Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient
348               END DO
349            END DO
350         END DO
351      ELSE
352         ! Fully explicit everywhere
353         Cu_adv = 0.0_wp                                ! Reuse array to output coefficient
354      ENDIF
355      CALL iom_put("wimp",wi) 
356      CALL iom_put("wi_cff",Cu_adv)
[10978]357      CALL iom_put("wexp",ww)
[10364]358      !
359      IF( ln_timing )   CALL timing_stop('wAimp')
360      !
361   END SUBROUTINE wAimp
[3]362   !!======================================================================
[1565]363END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.