Changeset 6041 for branches/2015/dev_r5776_UKMO2_OBS_efficiency_improvs/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

Timestamp:

2015-12-14T10:06:06+01:00 (8 years ago)

Author:

timgraham

Message:

Merged head of trunk into branch

File:

• branches/2015/dev_r5776_UKMO2_OBS_efficiency_improvs/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 (modified) (10 diffs)

branches/2015/dev_r5776_UKMO2_OBS_efficiency_improvs/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

@@ -20 +20 @@
    USE sbc_oce         ! surface boundary condition: ocean
    USE domvvl          ! Variable volume
-   USE divcur          ! hor. divergence and curl      (div & cur routines)
-   USE restart         ! only for lrst_oce
-   USE in_out_manager  ! I/O manager
-   USE prtctl          ! Print control
-   USE phycst
-   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
-   USE lib_mpp         ! MPP library
+   USE divhor          ! horizontal divergence
+   USE phycst          ! physical constants
    USE bdy_oce
    USE bdy_par         
@@ -36 +31 @@
    USE asminc          ! Assimilation increment
 #endif
+   USE in_out_manager  ! I/O manager
+   USE restart         ! only for lrst_oce
+   USE prtctl          ! Print control
+   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
+   USE lib_mpp         ! MPP library
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
@@ -66 +66 @@
       !!      by the time step.
       !!
-      !! ** action  :   ssha    : after sea surface height
+      !! ** action  :   ssha, after sea surface height
       !!
       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
       !!----------------------------------------------------------------------
-      !
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zhdiv
-      INTEGER, INTENT(in) ::   kt                      ! time step
+      INTEGER, INTENT(in) ::   kt   ! time step
       ! 
-      INTEGER             ::   jk                      ! dummy loop indice
-      REAL(wp)            ::   z2dt, z1_rau0           ! local scalars
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt')
-      !
-      CALL wrk_alloc( jpi, jpj, zhdiv ) 
+      INTEGER  ::   jk            ! dummy loop indice
+      REAL(wp) ::   z2dt, zcoef   ! local scalars
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhdiv   ! 2D workspace
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('ssh_nxt')
+      !
+      CALL wrk_alloc( jpi,jpj,   zhdiv ) 
       !
       IF( kt == nit000 ) THEN
-         !
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
          IF(lwp) WRITE(numout,*) '~~~~~~~ '
-         !
-      ENDIF
-      !
-      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
+      ENDIF
+      !
+      CALL div_hor( kt )                              ! Horizontal divergence
       !
       z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
@@ -106 +103 @@
       ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
       ! 
-      z1_rau0 = 0.5_wp * r1_rau0
-      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
-
-#if ! defined key_dynspg_ts
-      ! These lines are not necessary with time splitting since
-      ! boundary condition on sea level is set during ts loop
-#if defined key_agrif
-      CALL agrif_ssh( kt )
-#endif
-#if defined key_bdy
-      IF (lk_bdy) THEN
-         CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary
-         CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries
-      ENDIF
-#endif
-#endif
+      zcoef = 0.5_wp * r1_rau0
+      ssha(:,:) = (  sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
+
+      IF ( .NOT.ln_dynspg_ts ) THEN
+         ! These lines are not necessary with time splitting since
+         ! boundary condition on sea level is set during ts loop
+# if defined key_agrif
+         CALL agrif_ssh( kt )
+# endif
+# if defined key_bdy
+         IF( lk_bdy ) THEN
+            CALL lbc_lnk( ssha, 'T', 1. )    ! Not sure that's necessary
+            CALL bdy_ssh( ssha )             ! Duplicate sea level across open boundaries
+         ENDIF
+# endif
+      ENDIF
 
 #if defined key_asminc
-      !                                                ! Include the IAU weighted SSH increment
-      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
+      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN     ! Include the IAU weighted SSH increment
          CALL ssh_asm_inc( kt )
          ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
       ENDIF
 #endif
-
       !                                           !------------------------------!
       !                                           !           outputs            !
@@ -159 +154 @@
       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
       !!----------------------------------------------------------------------
-      !
-      INTEGER, INTENT(in) ::   kt           ! time step
+      INTEGER, INTENT(in) ::   kt   ! time step
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   z1_2dt       ! local scalars
       REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
-      !
-      INTEGER             ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp)            ::   z1_2dt       ! local scalars
-      !!----------------------------------------------------------------------
-      
-      IF( nn_timing == 1 )  CALL timing_start('wzv')
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('wzv')
       !
       IF( kt == nit000 ) THEN
-         !
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
@@ -177 +170 @@
          !
          wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
-         !
       ENDIF
       !                                           !------------------------------!
@@ -193 +185 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zhdiv(ji,jj,jk) = r1_e12t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
+                  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) )
                END DO
             END DO
@@ -216 +208 @@
 
 #if defined key_bdy
-      IF (lk_bdy) THEN
+      IF( lk_bdy ) THEN
          DO jk = 1, jpkm1
             wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
@@ -224 +216 @@
       !
       IF( nn_timing == 1 )  CALL timing_stop('wzv')
-
-
+      !
    END SUBROUTINE wzv
+
 
    SUBROUTINE ssh_swp( kt )
@@ -258 +250 @@
       ENDIF
 
-# if defined key_dynspg_ts
-      IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN    !** Euler time-stepping: no filter
-# else
-      IF ( neuler == 0 .AND. kt == nit000 ) THEN   !** Euler time-stepping at first time-step : no filter
-#endif
+      IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ( ln_bt_fw .AND. ln_dynspg_ts ) ) THEN 
+                                                   !** Euler time-stepping: no filter
          sshb(:,:) = sshn(:,:)                           ! before <-- now
          sshn(:,:) = ssha(:,:)                           ! now    <-- after  (before already = now)
+         !
       ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
          sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered