Changeset 12095


Ignore:
Timestamp:
2019-12-06T15:34:21+01:00 (8 months ago)
Author:
deazer
Message:

First iteration on correcting u or v depth points when using enveloping bathy
This needs tobe done also for cases other than vec
Code is very ineffiecint as written all in dynspg_ts , this should be rewritten
to only calculate some variables at run start.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0.1_MIRROR_WAD_ZENV/src/OCE/DYN/dynspg_ts.F90

    r12083 r12095  
    152152      REAL(wp) ::   r1_2dt_b, z1_hu, z1_hv          ! local scalars 
    153153      REAL(wp) ::   za0, za1, za2, za3              !   -      - 
     154!CEOD 
     155      REAL(wp), DIMENSION(jpi,jpj) :: e3u_0_tot, e3v_0_tot 
     156      REAL(wp), DIMENSION(jpi,jpj) :: sum_e3t_0_min_ik 
     157      REAL(wp), DIMENSION(jpi,jpj) :: sum_e3t_0_min_jk 
     158      REAL(wp), DIMENSION(jpi,jpj) :: scaled_e3u_0_tot, scaled_e3v_0_tot 
     159      REAL(wp), DIMENSION(jpi,jpj) :: k_bot 
     160      REAL(wp), DIMENSION(jpi,jpj) :: k_bot_i_min 
     161      REAL(wp), DIMENSION(jpi,jpj) :: k_bot_j_min 
     162      REAL(wp), DIMENSION(jpi,jpj) :: zdep_u, zdep_v 
     163!CEOD 
    154164      REAL(wp) ::   zmdi, zztmp, zldg               !   -      - 
    155165      REAL(wp) ::   zhu_bck, zhv_bck, zhdiv         !   -      - 
     
    181191!     zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp 
    182192      !                                         ! inverse of baroclinic time step  
     193!CEOD 
     194         e3u_0_tot(:,:)=0.0 
     195         e3v_0_tot(:,:)=0.0 
     196!CEOD 
     197         DO jk=1,jpk 
     198            e3u_0_tot(:,:) = e3u_0_tot(:,:) + e3u_0(:,:,jk) 
     199            e3v_0_tot(:,:) = e3v_0_tot(:,:) + e3v_0(:,:,jk) 
     200         ENDDO  
     201         k_bot(:,:) = 1. 
     202         DO jk=1,jpkm1 
     203          DO jj = 1, jpj 
     204            DO ji = 1, jpi             ! SPG with the application of W/D gravity filters 
     205              IF ( tmask(ji,jj,jk) .ne. 0 ) THEN ! If its a wet point we keep going down 
     206                 k_bot(ji,jj) = jk 
     207              ENDIF 
     208            ENDDO 
     209           ENDDO 
     210         ENDDO  
     211!work out k_bot min 
     212         DO jj = 1, jpjm1 
     213            DO ji = 1, jpim1                ! SPG with the application of W/D gravity filters 
     214                 k_bot_i_min(ji,jj) = MIN( k_bot(ji,jj), k_bot(ji+1,jj  )) 
     215                 k_bot_j_min(ji,jj) = MIN( k_bot(ji,jj), k_bot(ji,  jj+1)) 
     216            ENDDO 
     217        ENDDO 
     218        sum_e3t_0_min_ik(:,:) = 0 
     219        sum_e3t_0_min_jk(:,:) = 0 
     220 
     221        !Get sum of e3t_0s sown to local min 
     222        DO jj = 1, jpjm1 
     223           DO ji = 1, jpim1                
     224             DO jk = 1, k_bot_i_min(ji,jj) 
     225                sum_e3t_0_min_ik(ji,jj) =  sum_e3t_0_min_ik(ji,jj) + e3t_0(ji,jj,jk) 
     226             ENDDO 
     227            
     228             DO jk = 1, k_bot_j_min(ji,jj) 
     229                sum_e3t_0_min_jk(ji,jj) =  sum_e3t_0_min_jk(ji,jj) + e3t_0(ji,jj,jk) 
     230             ENDDO 
     231             !CEOD Now scale that against ht_0 will be one in the case it matches of course(Potential for ht_0 = zero) 
     232             scaled_e3u_0_tot(ji,jj) = sum_e3t_0_min_ik(ji,jj)/( ht_0(ji,jj) + 1._wp - ssmask(ji,jj)) 
     233             scaled_e3v_0_tot(ji,jj) = sum_e3t_0_min_jk(ji,jj)/( ht_0(ji,jj) + 1._wp - ssmask(ji,jj)) 
     234!CEOD             WRITE(numout,*)'scaled values', ji,jj ,scaled_e3u_0_tot(ji,jj) , sum_e3t_0_min_ik(ji,jj), ht_0(ji,jj) + 1._wp - ssmask(ji,jj) 
     235          ENDDO 
     236       ENDDO 
     237!CEO to join things up 
     238         CALL lbc_lnk_multi( 'dynspg_ts', scaled_e3u_0_tot, 'U', -1._wp ) 
     239         CALL lbc_lnk_multi( 'dynspg_ts', scaled_e3v_0_tot, 'V', -1._wp ) 
     240               
     241 
    183242      IF( kt == nit000 .AND. neuler == 0 ) THEN   ;   r1_2dt_b = 1._wp / (         rdt ) 
    184243      ELSE                                        ;   r1_2dt_b = 1._wp / ( 2._wp * rdt ) 
     
    463522            DO jj = 1, jpj 
    464523               DO ji = 1, jpim1   ! not jpi-column 
    465                   zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj)                        & 
    466                        &                              * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  & 
    467                        &                                 + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj) 
     524!CEOD                  zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj)                        & 
     525!CEOD                       &                              * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  & 
     526!CEOD                       &                                 + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj) 
     527                  zhup2_e(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji  ,jj) * zhtp2_e(ji  ,jj)*scaled_e3u_0_tot(ji  ,jj) & 
     528                                                 &         +  e1e2t(ji+1,jj) * zhtp2_e(ji+1,jj)*scaled_e3u_0_tot(ji+1,jj) ) *ssumask(ji,jj) 
    468529               END DO 
    469530            END DO 
    470531            DO jj = 1, jpjm1        ! not jpj-row 
    471532               DO ji = 1, jpi 
    472                   zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj)                        & 
    473                        &                              * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  & 
    474                        &                                 + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj) 
     533!CEOD                  zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj)                        & 
     534!CEOD                       &                              * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  & 
     535!CEOD                       &                                 + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj) 
     536                  zhvp2_e(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj  ) * zhtp2_e(ji,  jj)*scaled_e3v_0_tot(ji,jj  ) & 
     537                                                     &     +  e1e2t(ji,jj+1) * zhtp2_e(ji,jj+1)*scaled_e3v_0_tot(ji,jj+1) ) *ssvmask(ji,jj) 
    475538               END DO 
    476539            END DO 
     
    651714                  !                    ! hu_e, hv_e hold depth at jn,  zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2 
    652715                  !                    ! backward interpolated depth used in spg terms at jn+1/2 
    653                   zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)    & 
    654                        &                                          + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj) 
    655                   zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )    & 
    656                        &                                          + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj) 
     716                  !CEODzhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)    & 
     717                  zhu_bck = (hu_0(ji,jj)/e3u_0_tot(ji,jj))*(e3u_0_tot(ji,jj)+ r1_2*r1_e1e2u(ji,jj) * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)    & 
     718                       &                                          + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj) ) 
     719                  !CEODzhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )    & 
     720                  zhv_bck = (hv_0(ji,jj)/e3v_0_tot(ji,jj))*(e3v_0_tot(ji,jj) + r1_2*r1_e1e2v(ji,jj) * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )    & 
     721                       &                                          + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj) ) 
    657722                  !                    ! inverse depth at jn+1 
    658                   z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) ) 
    659                   z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) ) 
     723                  !CEODz1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) ) 
     724                  !CEODz1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) ) 
     725                  z1_hu = ssumask(ji,jj) / ( (hu_0(ji,jj)/e3u_0_tot(ji,jj))*(e3u_0_tot(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj)) )  
     726                  z1_hv = ssvmask(ji,jj) / ( (hv_0(ji,jj)/e3v_0_tot(ji,jj))*(e3v_0_tot(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj)) ) 
    660727                  ! 
    661728                  ua_e(ji,jj) = (               hu_e  (ji,jj) *   un_e (ji,jj)      &  
     
    682749        
    683750         IF( .NOT.ln_linssh ) THEN   !* Update ocean depth (variable volume case only) 
    684             hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1) 
     751!CEOD TBD 
     752         DO jj = 1, jpjm1 
     753            DO ji = 1, jpim1      ! NO Vector Opt. 
     754               zdep_u(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji,jj  ) * (ssha_e(ji,  jj)+ht_0(ji,  jj))*scaled_e3u_0_tot(ji  ,jj) & 
     755                                                 &     +  e1e2t(ji+1,jj) * (ssha_e(ji+1,jj)+ht_0(ji+1,jj))*scaled_e3u_0_tot(ji+1,jj) )*ssumask(ji,jj) 
     756               zdep_v(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj  ) * (ssha_e(ji,  jj)+ht_0(ji,  jj))*scaled_e3v_0_tot(ji,jj  ) & 
     757                                                 &     +  e1e2t(ji,jj+1) * (ssha_e(ji,jj+1)+ht_0(ji,jj+1))*scaled_e3v_0_tot(ji,jj+1) )*ssvmask(ji,jj) 
     758 
     759!CEOD             WRITE(numout,*)'zdep_u', ji,jj ,zdep_u(ji,jj) ,scaled_e3u_0_tot(ji,jj),ssha_e(ji,jj),jn 
     760            ENDDO 
     761        ENDDO 
     762         CALL lbc_lnk_multi( 'dynspg_ts', zdep_u, 'U', -1._wp ) 
     763         CALL lbc_lnk_multi( 'dynspg_ts', zdep_v, 'V', -1._wp ) 
     764 
     765            !CEODhu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1) 
     766            hu_e (2:jpim1,2:jpjm1) = zdep_u(2:jpim1,2:jpjm1)  
    685767            hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) ) 
    686             hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1) 
     768            !CEODhv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1) 
     769            hv_e (2:jpim1,2:jpjm1) = zdep_v(2:jpim1,2:jpjm1) 
    687770            hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) ) 
    688771            CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp  & 
     
    782865                  &              * ( e1e2t(ji,jj  ) * ssha(ji,jj  )      & 
    783866                  &              +   e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 
    784             END DO 
    785          END DO 
     867               zdep_u(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji,jj  ) * (ssha(ji,  jj)+ht_0(ji,  jj))*scaled_e3u_0_tot(ji  ,jj) & 
     868                                                 &     +  e1e2t(ji+1,jj) * (ssha(ji+1,jj)+ht_0(ji+1,jj))*scaled_e3u_0_tot(ji+1,jj) ) 
     869               zdep_v(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj  ) * (ssha(ji,  jj)+ht_0(ji,  jj))*scaled_e3v_0_tot(ji,jj  ) & 
     870                                                 &     +  e1e2t(ji,jj+1) * (ssha(ji,jj+1)+ht_0(ji,jj+1))*scaled_e3v_0_tot(ji,jj+1) ) 
     871            END DO 
     872         END DO 
     873         CALL lbc_lnk_multi( 'dynspg_ts', zdep_u, 'U', -1._wp ) 
     874         CALL lbc_lnk_multi( 'dynspg_ts', zdep_v, 'V', -1._wp ) 
     875 
    786876         CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 
    787877         ! 
     
    791881         END DO 
    792882         ! Save barotropic velocities not transport: 
    793          ua_b(:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 
    794          va_b(:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 
     883!CEOD         ua_b(:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 
     884!CEOD         va_b(:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 
     885         ua_b(:,:) =  ua_b(:,:) / ( zdep_u(:,:) + 1._wp - ssumask(:,:) ) 
     886         va_b(:,:) =  va_b(:,:) / ( zdep_v(:,:) + 1._wp - ssvmask(:,:) ) 
    795887      ENDIF 
    796888 
Note: See TracChangeset for help on using the changeset viewer.