Changeset 8842 for branches/2017/dev_r8624_ENHANCE4_FREESURFACE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

Timestamp:

2017-11-29T10:42:39+01:00 (7 years ago)

Author:

jchanut

Message:

Enhance4-freesurface. step 3: diffusive barotropic time stepping - #1959

File:

• branches/2017/dev_r8624_ENHANCE4_FREESURFACE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 (modified) (5 diffs)

branches/2017/dev_r8624_ENHANCE4_FREESURFACE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

@@ -16 +16 @@
    !!             3.7  ! 2015-11  (J. Chanut) free surface simplification
    !!              -   ! 2016-12  (G. Madec, E. Clementi) update for Stoke-Drift divergence
+   !!              -   ! 2017-11  (J. Chanut) add optional temporal diffusion 
    !!---------------------------------------------------------------------
 
@@ -150 +151 @@
       REAL(wp) ::   zhura, zhvra          !   -      -
       REAL(wp) ::   za0, za1, za2, za3    !   -      -
+      REAL(wp) ::   zepsilon, zgamma            !   -      -
       !
       REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e
@@ -758 +760 @@
            za3= 0._wp              
          ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 
-           za0=0.614_wp                     ! za0 = 1/2 +   gam + 2*eps    
-           za1=0.285_wp                     ! za1 = 1/2 - 2*gam - 3*eps 
-           za2=0.088_wp                     ! za2 = gam
-           za3=0.013_wp                     ! za3 = eps
+            IF (rn_bt_alpha==0._wp) THEN
+               za0=0.614_wp                 ! za0 = 1/2 +   gam + 2*eps
+               za1=0.285_wp                 ! za1 = 1/2 - 2*gam - 3*eps
+               za2=0.088_wp                 ! za2 = gam
+               za3=0.013_wp                 ! za3 = eps
+            ELSE
+               zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
+               zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
+               za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
+               za1 = 1._wp - za0 - zgamma - zepsilon
+               za2 = zgamma
+               za3 = zepsilon
+            ENDIF
          ENDIF
          !
@@ -1021 +1032 @@
       zwx(:,:) = un_adv(:,:)
       zwy(:,:) = vn_adv(:,:)
+      IF  ( kt == nit000 .AND. neuler==0 ) un_bf(:,:) = 0._wp ; vn_bf(:,:) = 0._wp
+      !     
       IF( ( kt == nit000 .AND. neuler==0 ) .OR. .NOT.ln_bt_fw ) THEN     
          un_adv(:,:) = zwx(:,:) * r1_hu_n(:,:)
          vn_adv(:,:) = zwy(:,:) * r1_hv_n(:,:)
-         !
-         un_bf(:,:) = 0._wp
-         vn_bf(:,:) = 0._wp 
       ELSE
          un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) ) * r1_hu_n(:,:)
@@ -1323 +1333 @@
       IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax
       !
+      IF(lwp) WRITE(numout,*)    '     Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
+
+      IF (rn_bt_alpha < 0._wp ) THEN
+         CALL ctl_stop( 'dynspg_ts ERROR: rn_bt_alpha < 0, it must be >= 0.' )
+      ENDIF
+
+      IF ( ln_bt_av.AND.(rn_bt_alpha > 0._wp) ) THEN
+         CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
+      ENDIF
+      !
       IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
          CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )