Changeset 4312 for branches

Timestamp:

2013-11-22T12:21:18+01:00 (10 years ago)

Author:

cetlod

Message:

time-splitting corrections

Location:

branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC

Files:

• DYN/dynnxt.F90 (modified) (7 diffs)
• step.F90 (modified) (2 diffs)

branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

@@ -106 +106 @@
       REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec   ! local scalars
       REAL(wp) ::   zve3a, zve3n, zve3b, zvf        !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:)   ::  zua, zva
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zua, zva, zhura, zhvra
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f 
       !!----------------------------------------------------------------------
@@ -113 +113 @@
       !
       CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
-      IF ( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zua, zva )
+      IF ( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zua, zva, zhura, zhvra )
       !
       IF( kt == nit000 ) THEN
@@ -162 +162 @@
       zva(:,:) = 0._wp
       IF (lk_vvl) THEN
+         zhura(:,:) = 0._wp
+         zhvra(:,:) = 0._wp
          DO jk = 1, jpkm1
             zua(:,:) = zua(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
-            zva(:,:) = zva(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)    
-         END DO
-         DO jk = 1, jpkm1
-            ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * hur_e(:,:) + ua_b(:,:) ) * umask(:,:,jk)
-            va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * hvr_e(:,:) + va_b(:,:) ) * vmask(:,:,jk)
+            zva(:,:) = zva(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
+            zhura(:,:) = zhura(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)
+            zhvra(:,:) = zhvra(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk)  
+         END DO
+         zhura(:,:) = umask(:,:,1) / ( zhura(:,:) + 1._wp - umask(:,:,1) ) 
+         zhvra(:,:) = vmask(:,:,1) / ( zhvra(:,:) + 1._wp - vmask(:,:,1) )  
+         DO jk = 1, jpkm1
+            ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * zhura(:,:) + ua_b(:,:) ) * umask(:,:,jk)
+            va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * zhvra(:,:) + va_b(:,:) ) * vmask(:,:,jk)
          END DO
       ELSE
@@ -184 +190 @@
          ! Remove advective velocity from "now velocities" 
          ! prior to asselin filtering     
-         ! In the forward case, this is done below after asselin filtering    
+         ! In the forward case, this is done below after asselin filtering   
+         ! so that asselin contribution is removed at the same time 
          DO jk = 1, jpkm1
             un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
@@ -308 +315 @@
                DO jk = 1, jpkm1
                   DO jj = 1, jpj
-                     DO ji = 1, jpi                   !                              ---------------------------
+                     DO ji = 1, jpi                  
                         zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
                         zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
@@ -333 +340 @@
          !
          IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
-         ! Remove asselin filtering of barotropic velocities if forward time splitting
-         ! note that we replace barotropic velocities by advective velocities       
+            ! Revert "before" velocities to time split estimate
+            ! Doing it here also means that asselin filter contribution is removed  
             zua(:,:) = 0._wp
             zva(:,:) = 0._wp
@@ -360 +367 @@
       ! 
       CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
-      IF ( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zua, zva )
+      IF ( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zua, zva, zhura, zhvra )
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')

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

@@ -132 +132 @@
                                   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)
+                                  ua_sv(:,:,:) = ua(:,:,:)     ! Save trends (barotropic trend has been fully updated at this stage)
                                   va_sv(:,:,:) = va(:,:,:)
 
@@ -140 +138 @@
       ENDIF
       IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors
-                         CALL wzv           ( kstp )  ! now cross-level velocity (original)
+                         CALL wzv           ( kstp )  ! now cross-level velocity 
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>