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 12777 – NEMO

Changeset 12777


Ignore:
Timestamp:
2020-04-20T13:31:01+02:00 (4 years ago)
Author:
davestorkey
Message:

UKMO/NEMO_4.0.1_biharmonic_GM : Remove duplicate dyn_ldf_lap_no_ahm routine. This version bit compares with previous version.

Location:
NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN/dynldf.F90

    r12776 r12777  
    6565      SELECT CASE ( nldf_dyn )                   ! compute lateral mixing trend and add it to the general trend 
    6666      ! 
    67       CASE ( np_lap   )    ;   CALL dyn_ldf_lap( kt, ub, vb, ua, va, 1 )      ! iso-level    laplacian 
     67      CASE ( np_lap   )    ;   CALL dyn_ldf_lap( kt, ahmt, ahmf, ub, vb, ua, va, 1 )      ! iso-level    laplacian 
    6868      CASE ( np_lap_i )    ;   CALL dyn_ldf_iso( kt )                         ! rotated      laplacian 
    6969      CASE ( np_blp   )    ;   CALL dyn_ldf_blp( kt, ub, vb, ua, va    )      ! iso-level bi-laplacian 

  • NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN/dynldf_lap_blp.F90

    r12775 r12777  
    2828   PUBLIC dyn_ldf_blp  ! called by dynldf.F90 
    2929   PUBLIC dyn_ldf_bgm  ! called by dynldf.F90 
    30    PRIVATE dyn_ldf_lap_no_ahm !needed by dyn_ldf_bgm 
    3130 
    3231   !! * Substitutions 
    … …  
    3938CONTAINS 
    4039 
    41    SUBROUTINE dyn_ldf_lap( kt, pub, pvb, pua, pva, kpass ) 
     40   SUBROUTINE dyn_ldf_lap( kt, pahmt, pahmf, pub, pvb, pua, pva, kpass ) 
    4241      !!---------------------------------------------------------------------- 
    4342      !!                     ***  ROUTINE dyn_ldf_lap  *** 
    … …  
    5352      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index 
    5453      INTEGER                         , INTENT(in   ) ::   kpass      ! =1/2 first or second passage 
     54      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pahmt, pahmf ! viscosity coefficients 
    5555      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity  [m/s] 
    5656      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! velocity trend   [m/s2] 
    … …  
    7979               !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1) 
    8080!!gm open question here : e3f  at before or now ?    probably now... 
    81 !!gm note that ahmf has already been multiplied by fmask 
    82                zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f_n(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       & 
     81!!gm note that pahmf has already been multiplied by fmask 
     82               zcur(ji-1,jj-1) = pahmf(ji-1,jj-1,jk) * e3f_n(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       & 
    8383                  &     * (  e2v(ji  ,jj-1) * pvb(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk)  & 
    8484                  &        - e1u(ji-1,jj  ) * pub(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk)  ) 
    8585               !                                      ! ahm * div        (computed from 2 to jpi/jpj) 
    86 !!gm note that ahmt has already been multiplied by tmask 
    87                zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk)                                         & 
     86!!gm note that pahmt has already been multiplied by tmask 
     87               zdiv(ji,jj)     = pahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk)                                         & 
    8888                  &     * (  e2u(ji,jj)*e3u_b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk)  & 
    8989                  &        + e1v(ji,jj)*e3v_b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*e3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk)  ) 
    … …  
    108108   END SUBROUTINE dyn_ldf_lap 
    109109 
    110 !CW: new subroutine for the Laplacian of velocity, copied from dyn_ldf_lap above, but without eddy viscosity ahmf, ahmt 
    111  
    112    SUBROUTINE dyn_ldf_lap_no_ahm( kt, pub, pvb, pulap, pvlap, kpass ) 
    113       !!---------------------------------------------------------------------- 
    114       !!                     ***  ROUTINE dyn_ldf_lap_no_ahm  *** 
    115       !!                        
    116       !! ** Purpose :   Companion subroutine to dyn_ldf_bgm to calculate  
    117       !!      the before horizontal momentum Laplacian  
    118       !!      trend and add it to the general trend of momentum equation.   
    119       !!      Note: mimics dyn_ldf_lap but without eddy viscosity ahmf, ahmt,  
    120       !!      because biharmonic GM eddy diffusivity is applied in dyn_bgm. 
    121       !! 
    122       !! ** Method  :   The Laplacian operator apply on horizontal velocity is  
    123       !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) )  
    124       !! 
    125       !! ** Action : - pua, pva increased by the harmonic operator applied on pub, pvb. 
    126       !!---------------------------------------------------------------------- 
    127       INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index 
    128       INTEGER                         , INTENT(in   ) ::   kpass      ! =1/2 first or second passage 
    129       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity  [m/s] 
    130       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(out) ::   pulap, pvlap ! velocity trend   [m/s2] 
    131       ! 
    132       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    133       REAL(wp) ::   zsign        ! local scalars 
    134       REAL(wp) ::   zua, zva     ! local scalars 
    135       REAL(wp), DIMENSION(jpi,jpj) ::   zcur, zdiv 
    136       !!---------------------------------------------------------------------- 
    137       ! 
    138       IF( kt == nit000 .AND. lwp ) THEN 
    139          WRITE(numout,*) 
    140          WRITE(numout,*) 'dyn_ldf_lap_no_ahm : companion iso-level harmonic (laplacian) operator to dyn_bgm, pass=', kpass 
    141          WRITE(numout,*) '~~~~~~~ ' 
    142       ENDIF 
    143       ! 
    144       IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign 
    145       ELSE                    ;   zsign = -1._wp      !  (eddy viscosity coef. >0) 
    146       ENDIF 
    147       ! 
    148       !                                                ! =============== 
    149       DO jk = 1, jpkm1                                 ! Horizontal slab 
    150          !                                             ! =============== 
    151          DO jj = 2, jpj 
    152             DO ji = fs_2, jpi   ! vector opt. 
    153                !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1) 
    154 !!gm open question here : e3f  at before or now ?    probably now... 
    155  
    156 !!CW: note that the removed ahmf had already been multiplied by fmask, so we multiply by fmask explicitly here, 
    157 !!    with the i,j,k of fmask aligning with those of ahmf, as defined in ldfdyn.F90. 
    158                zcur(ji-1,jj-1) = fmask(ji-1,jj-1,jk) * e3f_n(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)      & 
    159                   &     * (  e2v(ji  ,jj-1) * pvb(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk)  & 
    160                   &        - e1u(ji-1,jj  ) * pub(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk)  ) 
    161                !                                      ! ahm * div        (computed from 2 to jpi/jpj) 
    162 !!CW: note that the removed ahmt had already been multiplied by tmask, so we multiply by tmask explicitly here, 
    163 !!    with the i,j,k of tmask aligning with those of ahmt, as defined in ldfdyn.F90 
    164                zdiv(ji,jj)     = tmask(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk)                                         & 
    165                   &     * (  e2u(ji,jj)*e3u_b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk)  & 
    166                   &        + e1v(ji,jj)*e3v_b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*e3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk)  ) 
    167             END DO   
    168          END DO   
    169          ! 
    170 !CW: multiply pulap, pvlap by umask, vmask 
    171          DO jj = 2, jpjm1                             ! - curl( curl) + grad( div ) 
    172             DO ji = fs_2, fs_jpim1   ! vector opt. 
    173                pulap(ji,jj,jk) = umask(ji,jj,jk) * zsign * (                                            & 
    174                   &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk)   & 
    175                   &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                     ) 
    176                   ! 
    177                pvlap(ji,jj,jk) = vmask(ji,jj,jk) * zsign * (                                             & 
    178                   &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk)   & 
    179                   &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                     ) 
    180             END DO 
    181          END DO 
    182          !                                             ! =============== 
    183       END DO                                           !   End of slab 
    184       !                                                ! =============== 
    185       ! 
    186    END SUBROUTINE dyn_ldf_lap_no_ahm 
    187  
    188  
    189  
    190  
    191 !CW: new subroutine for biharmonic GM, copied from SUBROUTINE dyn_ldf_blp below 
    192  
    193110   SUBROUTINE dyn_ldf_bgm( kt, pub, pvb, pua, pva ) 
    194111      !!---------------------------------------------------------------------- 
    … …  
    223140      ENDIF 
    224141      ! 
    225    
     142      ! Calculate (del2 u) by calling dyn_ldf_lap with "viscosity" set to 1.0. 
     143      ! Normally we pass ahmt and ahmf to dyn_ldf_lap, which have been multiplied by tmask and fmask respectively.  
     144      ! So here we just pass tmask and fmask. 
     145      zulap(:,:,:) = 0.0_wp ; zvlap(:,:,:) = 0.0_wp 
     146      CALL dyn_ldf_lap( kt, tmask, fmask, pub, pvb, zulap, zvlap, 1 )    
     147      ! 
    226148      ! Calculate ratio of bhm and avm for stabilising correction terms 
    227149      ! and add to avm to be included in the vertical diffusion calculation later. 
     150      zmu(:,:,:) = 0.0_wp 
    228151      DO jk = 2, jpkm1  
    229152         DO jj = 2, jpjm1 
    … …  
    236159      CALL lbc_lnk_multi( 'dyn_ldf_bgm', zmu, 'W', 1. , avm, 'W', 1. )  
    237160       
    238       ! Calculate (del2 u) 
    239       CALL dyn_ldf_lap_no_ahm( kt, pub, pvb, zulap, zvlap, 1 )    
    240  
    241161      ! =============== 
    242162      !CW: calculate -bhm * d/dz(del^2 u)   
    … …  
    315235      zvlap(:,:,:) = 0._wp 
    316236      ! 
    317       CALL dyn_ldf_lap( kt, pub, pvb, zulap, zvlap, 1 )   ! rotated laplacian applied to ptb (output in zlap) 
     237      CALL dyn_ldf_lap( kt, ahmt, ahmf, pub, pvb, zulap, zvlap, 1 )   ! rotated laplacian applied to ptb (output in zlap) 
    318238      ! 
    319239      CALL lbc_lnk_multi( 'dynldf_lap_blp', zulap, 'U', -1., zvlap, 'V', -1. )             ! Lateral boundary conditions 
    320240      ! 
    321       CALL dyn_ldf_lap( kt, zulap, zvlap, pua, pva, 2 )   ! rotated laplacian applied to zlap (output in pta) 
     241      CALL dyn_ldf_lap( kt, ahmt, ahmf, zulap, zvlap, pua, pva, 2 )   ! rotated laplacian applied to zlap (output in pta) 
    322242      ! 
    323243   END SUBROUTINE dyn_ldf_blp 

  • NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/LDF/ldfdyn.F90

    r12776 r12777  
    180180      IF( ln_dynldf_blp ) THEN   ;                              ioptio = ioptio + 5   ;   ENDIF 
    181181      ! 
    182       IF( (ioptio < 1).OR.(ioptio > 6) )   CALL ctl_stop( 'dyn_ldf_init: use at least ONE of the 4 operator options (NONE/lap/blp) but not (lap+blp)' ) 
     182      IF( (ioptio < 1).OR.(ioptio > 6) )   CALL ctl_stop( 'dyn_ldf_init: use at least ONE of the 3 operator options (NONE/lap/blp) but not (lap+blp)' ) 
    183183      ! 
    184184 
Note: See TracChangeset for help on using the changeset viewer.