Ignore:
Timestamp:
2013-12-17T17:33:09+01:00 (7 years ago)
Author:
acc
Message:

Branch dev_MERGE_2013. #1209. Fix to compute after scale factors before wn computation in time-splitting case. dom_vvl_sf_nxt is now called twice but correctly computes the baroclinic contribution only once. Still need to sort out the asselin filtering of the separate components

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r4296 r4338  
    6161   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                            ! low frequency part of hz divergence 
    6262   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n           ! baroclinic scale factors 
    63    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a                        ! baroclinic scale factors 
     63   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a          ! baroclinic scale factors 
    6464   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                        ! retoring period for scale factors 
    6565   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                        ! retoring period for low freq. divergence 
     
    8282      IF( ln_vvl_zstar ) dom_vvl_alloc = 0 
    8383      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
    84          ALLOCATE( tilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   & 
    85             &      un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     , STAT = dom_vvl_alloc        ) 
     84         ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   & 
     85            &      dtilde_e3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   & 
     86            &      STAT = dom_vvl_alloc        ) 
    8687         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc ) 
    8788         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 
     
    215216               END DO 
    216217            END DO 
    217             IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN 
     218            IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN 
    218219               ii0 = 103   ;   ii1 = 111        ! Suppress ztilde in the Foxe Basin for ORCA2 
    219220               ij0 = 128   ;   ij1 = 135   ;    
     
    229230 
    230231 
    231    SUBROUTINE dom_vvl_sf_nxt( kt )  
     232   SUBROUTINE dom_vvl_sf_nxt( kt, kcall )  
    232233      !!---------------------------------------------------------------------- 
    233234      !!                ***  ROUTINE dom_vvl_sf_nxt  *** 
     
    255256      !! * Arguments 
    256257      INTEGER, INTENT( in )                  :: kt                    ! time step 
     258      INTEGER, INTENT( in ), OPTIONAL        :: kcall                 ! optional argument indicating call sequence 
    257259      !! * Local declarations 
    258260      INTEGER                                :: ji, jj, jk            ! dummy loop indices 
     
    260262      REAL(wp)                               :: z2dt                  ! temporary scalars 
    261263      REAL(wp)                               :: z_tmin, z_tmax        ! temporary scalars 
     264      LOGICAL                                :: ll_do_bclinic         ! temporary logical 
    262265      !!---------------------------------------------------------------------- 
    263266      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_nxt') 
     
    271274      ENDIF 
    272275 
     276      ll_do_bclinic = .TRUE. 
     277      IF( PRESENT(kcall) ) THEN 
     278         IF ( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE. 
     279      ENDIF 
     280 
    273281      ! ******************************* ! 
    274282      ! After acale factors at t-points ! 
    275283      ! ******************************* ! 
    276284 
    277       !                                                ! ----------------- ! 
    278       IF( ln_vvl_zstar ) THEN                          ! z_star coordinate ! 
    279          !                                             ! ----------------- ! 
    280  
    281          z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) ) 
    282          DO jk = 1, jpkm1 
    283             fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 
    284          END DO 
    285  
    286       !                                                ! --------------------------- ! 
    287       ELSEIF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN   ! z_tilde or layer coordinate ! 
    288          !                                             ! --------------------------- ! 
     285      !                                                ! --------------------------------------------- ! 
     286                                                       ! z_star coordinate and barotropic z-tilde part ! 
     287      !                                                ! --------------------------------------------- ! 
     288 
     289      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) ) 
     290      DO jk = 1, jpkm1 
     291         ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0) 
     292         fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 
     293      END DO 
     294 
     295      IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate ! 
     296         !                                                            ! ------baroclinic part------ ! 
    289297 
    290298         ! I - initialization 
     
    423431         tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) ) 
    424432 
    425          ! Add "tilda" part to the after scale factor 
     433         ! 
     434         ! "tilda" change in the after scale factor 
    426435         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
    427          fse3t_a(:,:,:) = e3t_0(:,:,:) + tilde_e3t_a(:,:,:) 
    428  
     436         DO jk = 1, jpkm1 
     437            dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk) 
     438         END DO 
    429439         ! III - Barotropic repartition of the sea surface height over the baroclinic profile 
    430440         ! ================================================================================== 
    431          ! add e3t(n-1) "star" Asselin-filtered 
    432          DO jk = 1, jpkm1 
    433             fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_b(:,:,jk) - e3t_0(:,:,jk) - tilde_e3t_b(:,:,jk) 
    434          END DO 
    435          ! add ( ssh increment + "baroclinicity error" ) proportionnaly to e3t(n) 
     441         ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n) 
    436442         ! - ML - baroclinicity error should be better treated in the future 
    437443         !        i.e. locally and not spread over the water column. 
     
    441447            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
    442448         END DO 
    443          z_scale(:,:) = ( ssha(:,:) - sshb(:,:) - zht(:,:) ) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) ) 
    444          DO jk = 1, jpkm1 
    445             fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 
    446          END DO 
    447  
    448       ENDIF 
    449  
    450       IF( ln_vvl_dbg ) THEN   ! - ML - test: control prints for debuging 
     449         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) ) 
     450         DO jk = 1, jpkm1 
     451            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 
     452         END DO 
     453 
     454      ENDIF 
     455 
     456      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate ! 
     457      !                                           ! ---baroclinic part--------- ! 
     458         DO jk = 1, jpkm1 
     459            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) 
     460         END DO 
     461      ENDIF 
     462 
     463      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging 
    451464         ! 
    452465         IF( lwp ) WRITE(numout, *) 'kt =', kt 
Note: See TracChangeset for help on using the changeset viewer.