Changeset 4338 for branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

Timestamp:

2013-12-17T17:33:09+01:00 (10 years ago)

Author:

acc

Message:

Branch dev_MERGE_2013. #1209. Fix to compute after scale factors before wn computation in time-splitting case. dom_vvl_sf_nxt is now called twice but correctly computes the baroclinic contribution only once. Still need to sort out the asselin filtering of the separate components

File:

• branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 (modified) (9 diffs)

branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

@@ -61 +61 @@
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                            ! low frequency part of hz divergence
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n           ! baroclinic scale factors
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a                        ! baroclinic scale factors
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a          ! baroclinic scale factors
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                        ! retoring period for scale factors
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                        ! retoring period for low freq. divergence
@@ -82 +82 @@
       IF( ln_vvl_zstar ) dom_vvl_alloc = 0
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
-         ALLOCATE( tilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   &
-            &      un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     , STAT = dom_vvl_alloc        )
+         ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   &
+            &      dtilde_e3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   &
+            &      STAT = dom_vvl_alloc        )
          IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
          IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
@@ -215 +216 @@
                END DO
             END DO
-            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN
+            IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN
                ii0 = 103   ;   ii1 = 111        ! Suppress ztilde in the Foxe Basin for ORCA2
                ij0 = 128   ;   ij1 = 135   ;   
@@ -229 +230 @@
 
 
-   SUBROUTINE dom_vvl_sf_nxt( kt ) 
+   SUBROUTINE dom_vvl_sf_nxt( kt, kcall ) 
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE dom_vvl_sf_nxt  ***
@@ -255 +256 @@
       !! * Arguments
       INTEGER, INTENT( in )                  :: kt                    ! time step
+      INTEGER, INTENT( in ), OPTIONAL        :: kcall                 ! optional argument indicating call sequence
       !! * Local declarations
       INTEGER                                :: ji, jj, jk            ! dummy loop indices
@@ -260 +262 @@
       REAL(wp)                               :: z2dt                  ! temporary scalars
       REAL(wp)                               :: z_tmin, z_tmax        ! temporary scalars
+      LOGICAL                                :: ll_do_bclinic         ! temporary logical
       !!----------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_nxt')
@@ -271 +274 @@
       ENDIF
 
+      ll_do_bclinic = .TRUE.
+      IF( PRESENT(kcall) ) THEN
+         IF ( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE.
+      ENDIF
+
       ! ******************************* !
       ! After acale factors at t-points !
       ! ******************************* !
 
-      !                                                ! ----------------- !
-      IF( ln_vvl_zstar ) THEN                          ! z_star coordinate !
-         !                                             ! ----------------- !
-
-         z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
-         DO jk = 1, jpkm1
-            fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
-         END DO
-
-      !                                                ! --------------------------- !
-      ELSEIF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN   ! z_tilde or layer coordinate !
-         !                                             ! --------------------------- !
+      !                                                ! --------------------------------------------- !
+                                                       ! z_star coordinate and barotropic z-tilde part !
+      !                                                ! --------------------------------------------- !
+
+      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
+      DO jk = 1, jpkm1
+         ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0)
+         fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
+      END DO
+
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
+         !                                                            ! ------baroclinic part------ !
 
          ! I - initialization
@@ -423 +431 @@
          tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
 
-         ! Add "tilda" part to the after scale factor
+         !
+         ! "tilda" change in the after scale factor
          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         fse3t_a(:,:,:) = e3t_0(:,:,:) + tilde_e3t_a(:,:,:)
-
+         DO jk = 1, jpkm1
+            dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
+         END DO
          ! III - Barotropic repartition of the sea surface height over the baroclinic profile
          ! ==================================================================================
-         ! add e3t(n-1) "star" Asselin-filtered
-         DO jk = 1, jpkm1
-            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_b(:,:,jk) - e3t_0(:,:,jk) - tilde_e3t_b(:,:,jk)
-         END DO
-         ! add ( ssh increment + "baroclinicity error" ) proportionnaly to e3t(n)
+         ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)
          ! - ML - baroclinicity error should be better treated in the future
          !        i.e. locally and not spread over the water column.
@@ -441 +447 @@
             zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
          END DO
-         z_scale(:,:) = ( ssha(:,:) - sshb(:,:) - zht(:,:) ) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
-         DO jk = 1, jpkm1
-            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
-         END DO
-
-      ENDIF
-
-      IF( ln_vvl_dbg ) THEN   ! - ML - test: control prints for debuging
+         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
+         DO jk = 1, jpkm1
+            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
+         END DO
+
+      ENDIF
+
+      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
+      !                                           ! ---baroclinic part--------- !
+         DO jk = 1, jpkm1
+            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk)
+         END DO
+      ENDIF
+
+      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging
          !
          IF( lwp ) WRITE(numout, *) 'kt =', kt