New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 4194 – NEMO

Changeset 4194


Ignore:
Timestamp:
2013-11-14T12:18:20+01:00 (11 years ago)
Author:
acc
Message:

Branch 2013/dev_r3858_NOC_ZTC, #863. Best endeavour merge of Z-tilde branch and Mercator time-stepping changes. Still issues with key_dynspg_ts that need all party involvement to resolve.

Location:
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r4178 r4194  
    172172      REAL(wp) ::   zx1, zy1, zx2, zy2                        !   -      - 
    173173      REAL(wp) ::   z1_12, z1_8, z1_4, z1_2 
    174       REAL(wp) ::   za1, za2, za3                             !   -      - 
     174      REAL(wp) ::   za0, za1, za2, za3                        !   -      - 
    175175      REAL(wp) ::   zu_spg, zu_cor, zu_sld, zu_asp            !   -      - 
    176176      REAL(wp) ::   zv_spg, zv_cor, zv_sld, zv_asp            !   -      - 
     
    182182      REAL(wp), POINTER, DIMENSION(:,:) :: zhu_b, zhv_b 
    183183      REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zsshp2_e 
     184      REAL(wp), POINTER, DIMENSION(:,:) :: zhust_e, zhvst_e 
    184185      !!---------------------------------------------------------------------- 
    185186      ! 
     
    191192      CALL wrk_alloc( jpi, jpj, zhu_b, zhv_b                                     ) 
    192193      CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zsshp2_e                       ) 
     194      CALL wrk_alloc( jpi, jpj, zhust_e, zhvst_e                                 ) 
    193195      ! 
    194196 
     
    385387      zcoef = -1._wp * z1_2dt_b 
    386388 
    387 !--------> MERGED TO HERE 
    388389      IF(ln_bfrimp) THEN 
    389390      !                                   ! Remove the bottom stress trend from 3-D sea surface level gradient 
     
    455456      ! Initialize depths: 
    456457      IF ( lk_vvl.AND.(.NOT.ln_bt_fw) ) THEN 
    457 !        hu_e   (:,:) = hu_0(:,:) + sshu_b(:,:)  
    458 !        hv_e   (:,:) = hv_0(:,:) + sshv_b(:,:) 
    459 !        hur_e  (:,:) = umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 
    460 !        hvr_e  (:,:) = vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 
    461          hu_e   (:,:) = hu   (:,:)        
    462          hv_e   (:,:) = hv   (:,:)  
    463          hur_e  (:,:) = hur  (:,:)     
    464          hvr_e  (:,:) = hvr  (:,:) 
     458         hur_e  (:,:) = zhu_b  (:,:)     
     459         hvr_e  (:,:) = zhv_b  (:,:) 
     460         hu_e(:,:) = umask(:,:,1) / ( zhu_b(:,:) + 1. - umask(:,:,1) ) 
     461         hv_e(:,:) = vmask(:,:,1) / ( zhv_b(:,:) + 1. - vmask(:,:,1) ) 
    465462      ELSE 
    466463         hu_e   (:,:) = hu   (:,:)        
     
    526523 
    527524         ! Extrapolate barotropic velocities at step jit+0.5: 
    528          ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 
    529          va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 
     525         ua_e(:,:) = za1 * zun_e(:,:) + za2 * zub_e(:,:) + za3 * ubb_e(:,:) 
     526         va_e(:,:) = za1 * zvn_e(:,:) + za2 * zvb_e(:,:) + za3 * vbb_e(:,:) 
    530527 
    531528         IF( lk_vvl ) THEN                                !* Update ocean depth (variable volume case only) 
    532529            !                                             !  ------------------ 
    533530            ! Extrapolate Sea Level at step jit+0.5: 
    534             zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 
     531            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * zsshb_e(:,:) + za3 * sshbb_e(:,:) 
    535532            ! 
    536533            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points 
     
    554551         ! considering fluxes below: 
    555552         ! 
    556          zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5 
     553         zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)   ! fluxes at jn+0.5  ! ACC INSTANT FAILURE IF THESE ARE USED (VERY LARGE NEGATIVE SALINITY AFTER 1 TS) 
    557554         zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 
     555!        zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:)     ! fluxes at jn+0.5  ! ACC WORKS BETTER BUT FAILS AFTER O(100) TIMESTEPS 
     556!        zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:) 
    558557         DO jj = 2, jpjm1                                  
    559558            DO ji = fs_2, fs_jpim1   ! vector opt. 
    560                zhdiv(ji,jj) = (   e2u(ji  ,jj) * zun_e(ji  ,jj) * hu_e(ji  ,jj)     & 
    561                   &             - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj)     & 
    562                   &             + e1v(ji,jj  ) * zvn_e(ji,jj  ) * hv_e(ji,jj  )     & 
    563                   &             - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1)   ) / ( e1t(ji,jj) * e2t(ji,jj) ) 
     559               zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   & 
     560                  &             + zwy(ji,jj) - zwy(ji,jj-1)   & 
     561                  &           ) / ( e1t(ji,jj) * e2t(ji,jj) ) 
    564562            END DO 
    565563         END DO 
     
    582580            END DO 
    583581         END DO 
    584  
    585          !                                                !* after barotropic velocities (vorticity scheme dependent) 
    586          !                                                !  ---------------------------   
    587          zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:)     ! now_e transport 
    588          zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:) 
     582#if defined key_bdy 
     583         ! Duplicate sea level across open boundaries (this is only cosmetic if lk_vvl=.false.) 
     584         IF (lk_bdy) CALL bdy_ssh( ssha_e ) 
     585#endif 
     586         !                                  
     587         ! Half-step back interpolation of SSH for surface pressure computation: 
     588         !---------------------------------------------------------------------- 
     589         IF ((jn==1).AND.ll_init) THEN 
     590           za0=1._wp                        ! Forward-backward 
     591           za1=0._wp                            
     592           za2=0._wp 
     593           za3=0._wp 
     594         ELSEIF ((jn==2).AND.ll_init) THEN  ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12 
     595           za0= 1.0833333333333_wp          ! za0 = 1-gam-eps 
     596           za1=-0.1666666666666_wp          ! za1 = gam 
     597           za2= 0.0833333333333_wp          ! za2 = eps 
     598           za3= 0._wp               
     599         ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880  
     600           za0=0.614_wp                     ! za0 = 1/2 +   gam + 2*eps     
     601           za1=0.285_wp                     ! za1 = 1/2 - 2*gam - 3*eps  
     602           za2=0.088_wp                     ! za2 = gam 
     603           za3=0.013_wp                     ! za3 = eps 
     604         ENDIF 
     605 
     606         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:) & 
     607          &            + za2 *  zsshb_e(:,:) + za3 *  sshbb_e(:,:) 
     608 
     609         ! 
     610         ! Compute associated depths at U and V points: 
     611         IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN       
     612            !                                         
     613            DO jj = 2, jpjm1                             
     614               DO ji = 2, jpim1 
     615                  zx1 = z1_2 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       & 
     616                    &        * ( e1t(ji  ,jj) * e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  & 
     617                    &        +   e1t(ji+1,jj) * e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )                  
     618                  zy1 = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       & 
     619                     &       * ( e1t(ji,jj  ) * e2t(ji,jj  ) * zsshp2_e(ji,jj  )  & 
     620                     &       +   e1t(ji,jj+1) * e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) 
     621                  zhust_e(ji,jj) = hu_0(ji,jj) + zx1  
     622                  zhvst_e(ji,jj) = hv_0(ji,jj) + zy1 
     623               END DO 
     624            END DO 
     625         ENDIF 
     626         ! 
     627         ! Add Coriolis trend: 
     628         ! zwz array below or triads normally depend on sea level with key_vvl and should be updated 
     629         ! at each time step. We however keep them constant here for optimization. 
     630         ! Recall that zwx and zwy arrays hold fluxes at this stage: 
     631         ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)   ! fluxes at jn+0.5 
     632         ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 
    589633         ! 
    590634         IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      !==  energy conserving or mixed scheme  ==! 
     
    761805         !                                                !* Time filter and swap 
    762806         !                                                !  -------------------- 
     807         sshbb_e(:,:) = zsshb_e(:,:) 
     808         ubb_e  (:,:) = zub_e  (:,:) 
     809         vbb_e  (:,:) = zvb_e  (:,:) 
    763810         IF( jn == 1 ) THEN                                     ! Swap only (1st Euler time step) 
    764811            zsshb_e(:,:) = sshn_e(:,:) 
     
    826873      ! 
    827874      !                                   !* Time average ==> after barotropic u, v, ssh 
    828       zcoef =  1._wp / ( 2 * nn_baro + 1 )  
     875      zcoef =  1._wp / ( icycle + 1 )  
    829876      zu_sum(:,:) = zcoef * zu_sum  (:,:)  
    830877      zv_sum(:,:) = zcoef * zv_sum  (:,:)  
     
    855902      CALL wrk_dealloc( jpi, jpj, zhu_b, zhv_b                                     ) 
    856903      CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zsshp2_e                       ) 
     904      CALL wrk_dealloc( jpi, jpj, zhust_e, zhvst_e                                 ) 
    857905      ! 
    858906      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts') 
     
    10521100            sshn_b(:,:) = sshb(:,:)   ! if not in restart set previous time mean to current baroclinic before value    
    10531101         ENDIF  
     1102 
     1103         IF( .NOT.ln_bt_av .AND. iom_varid( numror, 'sshbb_e', ldstop = .FALSE. ) > 0) THEN 
     1104            CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:) )    
     1105            CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:) )    
     1106            CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:) ) 
     1107            CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:) )  
     1108            CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:) )    
     1109            CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:) ) 
     1110         ELSE 
     1111            sshbb_e = sshn_b                                                ! ACC GUESS WORK 
     1112            ubb_e   = ub_b 
     1113            vbb_e   = vb_b 
     1114            sshb_e  = sshn_b 
     1115            ub_e    = ub_b 
     1116            vb_e    = vb_b 
     1117         ENDIF 
    10541118      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 
    10551119         CALL iom_rstput( kt, nitrst, numrow, 'un_b'   , un_b  (:,:) )   ! external velocity and ssh 
    10561120         CALL iom_rstput( kt, nitrst, numrow, 'vn_b'   , vn_b  (:,:) )   ! issued from barotropic loop 
    10571121         CALL iom_rstput( kt, nitrst, numrow, 'sshn_b' , sshn_b(:,:) )   !  
     1122         IF (.NOT.ln_bt_av) THEN 
     1123            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:) )  
     1124            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:) ) 
     1125            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:) ) 
     1126            CALL iom_rstput( kt, nitrst, numrow, 'sshb_e'   ,  sshb_e(:,:) ) 
     1127            CALL iom_rstput( kt, nitrst, numrow, 'ub_e'     ,    ub_e(:,:) ) 
     1128            CALL iom_rstput( kt, nitrst, numrow, 'vb_e'     ,    vb_e(:,:) ) 
     1129         ENDIF 
    10581130      ENDIF 
    10591131      ! 
     
    11541226 
    11551227#else 
    1156    !!---------------------------------------------------------------------- 
    1157    !!   Default case :   Empty module   No standart free surface cst volume 
    1158    !!---------------------------------------------------------------------- 
     1228   !!--------------------------------------------------------------------------- 
     1229   !!   Default case :   Empty module   No standard free surface constant volume 
     1230   !!--------------------------------------------------------------------------- 
    11591231 
    11601232   USE par_kind 
  • branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r4105 r4194  
    326326      END DO 
    327327      ! 
    328       DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  == 
     328      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  == 
    329329         DO ji = fs_2, fs_jpim1   ! vector opt. 
    330330            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
Note: See TracChangeset for help on using the changeset viewer.