Changeset 4338


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

Location:
branches/2013/dev_MERGE_2013/NEMOGCM
Files:
1 added
5 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_2/limdia_2.F90

    r4325 r4338  
    242242 
    243243       ! opening  "ice_evolu" file 
    244        CALL ctl_opn( numevo_ice, 'ice_evolu', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 
     244       CALL ctl_opn( numevo_ice, 'ice_evolu', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 
    245245 
    246246       !- ecriture de 2 lignes d''entete : 
  • 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 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

    r4328 r4338  
    252252            ! ---------------------------------------------------- 
    253253            IF (lk_dynspg_ts.AND.ln_bt_fw) THEN 
    254                ! Remove asselin filtering on thicknesses if forward time splitting 
     254               ! No asselin filtering on thicknesses if forward time splitting 
    255255                  fse3t_b(:,:,:) = fse3t_n(:,:,:) 
    256256            ELSE 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r4328 r4338  
    110110      END DO 
    111111      !                                                ! Sea surface elevation time stepping 
    112       ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used 
    113       ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp 
     112      ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to 
     113      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 
     114      !  
    114115      z1_rau0 = 0.5_wp * r1_rau0 
    115116      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1) 
     
    197198         DO jk = 1, jpkm1 
    198199            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t) 
    199             ! - ML - note: computation allready done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way) 
     200            ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way) 
    200201            DO jj = 2, jpjm1 
    201202               DO ji = fs_2, fs_jpim1   ! vector opt. 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/step.F90

    r4328 r4338  
    7272      INTEGER ::   jk       ! dummy loop indice 
    7373      INTEGER ::   indic    ! error indicator if < 0 
     74      INTEGER ::   kcall    ! optional integer argument (dom_vvl_sf_nxt) 
    7475      !! --------------------------------------------------------------------- 
    7576 
     
    106107                         CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_cur) 
    107108      IF( lk_dynspg_ts ) THEN 
     109          IF( lk_vvl     )        CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors  
    108110                                  CALL wzv           ( kstp )  ! now cross-level velocity  
    109111          ! In case the time splitting case, update almost all momentum trends here: 
    110112          ! Note that the computation of vertical velocity above, hence "after" sea level 
    111113          ! is necessary to compute momentum advection for the rhs of barotropic loop: 
    112                                   CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )                 ! now in situ density for hpg computation 
    113           IF( ln_zps      )       CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &  ! zps: now hor. derivative 
    114                 &                                          rhd, gru , grv  )     ! of t, s, rd at the last ocean level 
     114                                  CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 
     115          IF( ln_zps      )       CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &   ! zps: now hor. derivative 
     116                &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level 
    115117 
    116118                                  ua(:,:,:) = 0.e0             ! set dynamics trends to zero 
     
    135137                                  CALL div_cur( kstp )         ! Horizontal divergence & Relative vorticity (2nd call in time-split case) 
    136138      ENDIF 
    137       IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors 
     139      IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors 
    138140                         CALL wzv           ( kstp )  ! now cross-level velocity  
    139141 
Note: See TracChangeset for help on using the changeset viewer.