Changeset 4288

Timestamp:

2013-11-20T16:02:06+01:00 (10 years ago)

Author:

cbricaud

Message:

correction for timesplitting

Location:

branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC

Files:

• DYN/dynspg_ts.F90 (modified) (5 diffs)
• DYN/sshwzv.F90 (modified) (2 diffs)
• step.F90 (modified) (5 diffs)

branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

@@ -56 +56 @@
 
    ! Potential namelist parameters below to be read in dyn_spg_ts_init
-   LOGICAL,  PUBLIC,  PARAMETER :: ln_bt_fw=.true.        !: Forward integration of barotropic sub-stepping
+   LOGICAL,  PUBLIC,  PARAMETER :: ln_bt_fw=.TRUE.        !: Forward integration of barotropic sub-stepping
    LOGICAL,  PRIVATE, PARAMETER :: ln_bt_av=.TRUE.        !: Time averaging of barotropic variables
    LOGICAL,  PRIVATE, PARAMETER :: ln_bt_nn_auto=.FALSE.  !: Set number of iterations automatically
@@ -228 +228 @@
       IF ( kt == nit000 .OR. lk_vvl ) THEN
          IF ( ln_dynvor_een ) THEN
-            zht(:,:) = 0.
-            DO jk = 1, jpk
-               zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
-            END DO
+            ! JC: Simplification needed below: define ht_0 even when volume is fixed
+            IF (lk_vvl) THEN
+               zht(:,:) = (ht_0(:,:) + sshn(:,:)) * tmask(:,:,1) 
+            ELSE
+               zht(:,:) = 0.
+               DO jk = 1, jpkm1
+                  zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
+               END DO
+            ENDIF
+
             DO jj = 1, jpjm1
                DO ji = 1, jpim1
@@ -1049 +1055 @@
       INTEGER         , INTENT(in) ::   kt         ! ocean time-step
       !
-      INTEGER  :: ji ,jj
+      INTEGER  :: ji ,jj, jk
       REAL(wp) :: zxr2, zyr2, zcmax
-      REAL(wp), POINTER, DIMENSION(:,:) :: zcu
+      REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zht
       !!
 !      NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &
@@ -1060 +1066 @@
       !         ! Max courant number for ext. grav. waves
       !
-      CALL wrk_alloc( jpi, jpj, zcu )
-      !
+      CALL wrk_alloc( jpi, jpj, zcu, zht )
+      !
+      ! JC: Simplification needed below: define ht_0 even when volume is fixed
+      IF (lk_vvl) THEN 
+         zht(:,:) = ht_0(:,:) * tmask(:,:,1)
+      ELSE
+         zht(:,:) = 0.
+         DO jk = 1, jpkm1
+            zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
+         END DO
+      ENDIF
+
       DO jj = 1, jpj
-         DO ji =1, jpi 
+         DO ji =1, jpi
             zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))
             zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))
-            zcu(ji,jj) = sqrt(grav*ht_0(ji,jj)*(zxr2 + zyr2) )
+            zcu(ji,jj) = sqrt(grav*zht(ji,jj)*(zxr2 + zyr2) )
          END DO
       END DO
@@ -1120 +1136 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj, zcu )
+      CALL wrk_dealloc( jpi, jpj, zcu, zht )
       !
    END SUBROUTINE dyn_spg_ts_init

branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

@@ -36 +36 @@
    USE iom
    USE sbcrnf, ONLY: h_rnf, nk_rnf, sbc_rnf_div   ! River runoff 
-   USE dynspg_oce,  ONLY: un_b, vn_b, un_adv, vn_adv
    USE dynspg_ts,   ONLY: ln_bt_fw
    USE dynspg_oce, ONLY: lk_dynspg_ts
@@ -288 +287 @@
       ENDIF
 
-      IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN    !** Euler time-stepping at first time-step : no filter
+# 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
          sshb(:,:) = sshn(:,:)                           ! before <-- now
          sshn(:,:) = ssha(:,:)                           ! now    <-- after  (before already = now)

branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step.F90

@@ -99 +99 @@
       !  Ocean dynamics : hdiv, rot, ssh, e3, wn
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+                         CALL zdf_bfr( kstp )         ! bottom friction (if quadratic)
                          CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_cur)
       IF( lk_dynspg_ts ) THEN
@@ -108 +109 @@
           IF( ln_zps      )       CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &  ! zps: now hor. derivative
                 &                                          rhd, gru , grv  )     ! of t, s, rd at the last ocean level
-                                  CALL zdf_bfr( kstp )         ! bottom friction (if quadratic)
 
                                   ua(:,:,:) = 0.e0             ! set dynamics trends to zero
@@ -125 +125 @@
                                   CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure
                                   CALL dyn_spg( kstp, indic )  ! surface pressure gradient
+
+                                  hdivb(:,:,:) = hdivn(:,:,:)  ! Store now divergence and rot temporarly, revert to these below 
+                                  rotb(:,:,:)  = rotn(:,:,:)     
+                                  ua_sv(:,:,:) = ua(:,:,:)     ! Save trends (barotropic trend has been fully updated)
+                                  va_sv(:,:,:) = va(:,:,:)
+
                                   CALL div_cur( kstp )         ! Horizontal divergence & Relative vorticity (2nd call in time-split case)
-
-                                  ua_sv(:,:,:) = ua(:,:,:)     ! save next velocities (not trends !)
-                                  va_sv(:,:,:) = va(:,:,:)
       ENDIF
       IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors
@@ -140 +143 @@
       !
       !  VERTICAL PHYSICS
-      IF( .NOT. lk_dynspg_ts ) CALL zdf_bfr( kstp )         ! bottom friction
-
       !                                               ! Vertical eddy viscosity and diffusivity coefficients
       IF( lk_zdfric  )   CALL zdf_ric( kstp )            ! Richardson number dependent Kz
@@ -252 +253 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       IF( lk_dynspg_ts   )  THEN
-! revert to previously computed tendencies:
+! revert to previously computed momentum tendencies
 ! (not using ua, va as temporary arrays during tracers' update could avoid that)
                                ua(:,:,:) = ua_sv(:,:,:)
                                va(:,:,:) = va_sv(:,:,:)
+! Revert now divergence and rotational to previously computed ones 
+!(needed because of the time swap in div_cur, at the beginning of each time step)
+                               hdivn(:,:,:) = hdivb(:,:,:)
+                               rotn(:,:,:)  = rotb(:,:,:) 
+
                                CALL dyn_bfr( kstp )         ! bottom friction
                                CALL dyn_zdf( kstp )         ! vertical diffusion