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 12928 for NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/ZDF/zdfdrg.F90 – NEMO

Ignore:
Timestamp:
2020-05-14T21:46:00+02:00 (4 years ago)
Author:
smueller
Message:

Synchronizing with /NEMO/trunk@12925 (ticket #2170)

Location:
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser

    • Property svn:externals
      •  

        old new  
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
         8 
         9# SETTE 
         10^/utils/CI/sette@HEAD         sette 
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/ZDF/zdfdrg.F90

    r12178 r12928  
    7373 
    7474   !! * Substitutions 
    75 #  include "vectopt_loop_substitute.h90" 
     75#  include "do_loop_substitute.h90" 
    7676   !!---------------------------------------------------------------------- 
    7777   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    8181CONTAINS 
    8282 
    83    SUBROUTINE zdf_drg( kt, k_mk, pCdmin, pCdmax, pz0, pke0, pCd0,   &   ! <<== in  
     83   SUBROUTINE zdf_drg( kt, Kmm, k_mk, pCdmin, pCdmax, pz0, pke0, pCd0,   &   ! <<== in  
    8484      &                                                     pCdU )      ! ==>> out : bottom drag [m/s] 
    8585      !!---------------------------------------------------------------------- 
     
    9999      !!---------------------------------------------------------------------- 
    100100      INTEGER                 , INTENT(in   ) ::   kt       ! ocean time-step index 
     101      INTEGER                 , INTENT(in   ) ::   Kmm      ! ocean time level index 
    101102      !                       !               !!         !==  top or bottom variables  ==! 
    102103      INTEGER , DIMENSION(:,:), INTENT(in   ) ::   k_mk     ! wet level (1st or last) 
     
    114115      ! 
    115116      IF( l_log_not_linssh ) THEN     !==  "log layer"  ==!   compute Cd and -Cd*|U| 
    116          DO jj = 2, jpjm1 
    117             DO ji = 2, jpim1 
    118                imk = k_mk(ji,jj)          ! ocean bottom level at t-points 
    119                zut = un(ji,jj,imk) + un(ji-1,jj,imk)     ! 2 x velocity at t-point 
    120                zvt = vn(ji,jj,imk) + vn(ji,jj-1,imk) 
    121                zzz = 0.5_wp * e3t_n(ji,jj,imk)           ! altitude below/above (top/bottom) the boundary 
    122                ! 
     117         DO_2D_00_00 
     118            imk = k_mk(ji,jj)          ! ocean bottom level at t-points 
     119            zut = uu(ji,jj,imk,Kmm) + uu(ji-1,jj,imk,Kmm)     ! 2 x velocity at t-point 
     120            zvt = vv(ji,jj,imk,Kmm) + vv(ji,jj-1,imk,Kmm) 
     121            zzz = 0.5_wp * e3t(ji,jj,imk,Kmm)           ! altitude below/above (top/bottom) the boundary 
     122            ! 
    123123!!JC: possible WAD implementation should modify line below if layers vanish 
    124                zcd = (  vkarmn / LOG( zzz / pz0 )  )**2 
    125                zcd = pCd0(ji,jj) * MIN(  MAX( pCdmin , zcd ) , pCdmax  )   ! here pCd0 = mask*boost 
    126                pCdU(ji,jj) = - zcd * SQRT(  0.25 * ( zut*zut + zvt*zvt ) + pke0  ) 
    127             END DO 
    128          END DO 
     124            zcd = (  vkarmn / LOG( zzz / pz0 )  )**2 
     125            zcd = pCd0(ji,jj) * MIN(  MAX( pCdmin , zcd ) , pCdmax  )   ! here pCd0 = mask*boost 
     126            pCdU(ji,jj) = - zcd * SQRT(  0.25 * ( zut*zut + zvt*zvt ) + pke0  ) 
     127         END_2D 
    129128      ELSE                                            !==  standard Cd  ==! 
    130          DO jj = 2, jpjm1 
    131             DO ji = 2, jpim1 
    132                imk = k_mk(ji,jj)    ! ocean bottom level at t-points 
    133                zut = un(ji,jj,imk) + un(ji-1,jj,imk)     ! 2 x velocity at t-point 
    134                zvt = vn(ji,jj,imk) + vn(ji,jj-1,imk) 
    135                !                                                           ! here pCd0 = mask*boost * drag 
    136                pCdU(ji,jj) = - pCd0(ji,jj) * SQRT(  0.25 * ( zut*zut + zvt*zvt ) + pke0  ) 
    137             END DO 
    138          END DO 
    139       ENDIF 
    140       ! 
    141       IF(ln_ctl)   CALL prt_ctl( tab2d_1=pCdU, clinfo1=' Cd*U ') 
     129         DO_2D_00_00 
     130            imk = k_mk(ji,jj)    ! ocean bottom level at t-points 
     131            zut = uu(ji,jj,imk,Kmm) + uu(ji-1,jj,imk,Kmm)     ! 2 x velocity at t-point 
     132            zvt = vv(ji,jj,imk,Kmm) + vv(ji,jj-1,imk,Kmm) 
     133            !                                                           ! here pCd0 = mask*boost * drag 
     134            pCdU(ji,jj) = - pCd0(ji,jj) * SQRT(  0.25 * ( zut*zut + zvt*zvt ) + pke0  ) 
     135         END_2D 
     136      ENDIF 
     137      ! 
     138      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pCdU, clinfo1=' Cd*U ') 
    142139      ! 
    143140   END SUBROUTINE zdf_drg 
    144141 
    145142 
    146    SUBROUTINE zdf_drg_exp( kt, pub, pvb, pua, pva ) 
     143   SUBROUTINE zdf_drg_exp( kt, Kmm, pub, pvb, pua, pva ) 
    147144      !!---------------------------------------------------------------------- 
    148145      !!                  ***  ROUTINE zdf_drg_exp  *** 
     
    157154      !!--------------------------------------------------------------------- 
    158155      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index 
     156      INTEGER                         , INTENT(in   ) ::   Kmm        ! time level indices 
    159157      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pub, pvb   ! the two components of the before velocity 
    160158      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! the two components of the velocity tendency 
     
    167165      !!--------------------------------------------------------------------- 
    168166      ! 
    169 !!gm bug : time step is only rdt (not 2 rdt if euler start !) 
    170       zm1_2dt = - 1._wp / ( 2._wp * rdt ) 
     167!!gm bug : time step is only rn_Dt (not 2 rn_Dt if euler start !) 
     168      zm1_2dt = - 1._wp / ( 2._wp * rn_Dt ) 
    171169 
    172170      IF( l_trddyn ) THEN      ! trends: store the input trends 
     
    176174      ENDIF 
    177175 
    178       DO jj = 2, jpjm1 
    179          DO ji = 2, jpim1 
    180             ikbu = mbku(ji,jj)          ! deepest wet ocean u- & v-levels 
    181             ikbv = mbkv(ji,jj) 
     176      DO_2D_00_00 
     177         ikbu = mbku(ji,jj)          ! deepest wet ocean u- & v-levels 
     178         ikbv = mbkv(ji,jj) 
     179         ! 
     180         ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 
     181         zCdu = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / e3u(ji,jj,ikbu,Kmm) 
     182         zCdv = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / e3v(ji,jj,ikbv,Kmm) 
     183         ! 
     184         pua(ji,jj,ikbu) = pua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * pub(ji,jj,ikbu) 
     185         pva(ji,jj,ikbv) = pva(ji,jj,ikbv) + MAX(  zCdv , zm1_2dt  ) * pvb(ji,jj,ikbv) 
     186      END_2D 
     187      ! 
     188      IF( ln_isfcav ) THEN        ! ocean cavities 
     189         DO_2D_00_00 
     190            ikbu = miku(ji,jj)          ! first wet ocean u- & v-levels 
     191            ikbv = mikv(ji,jj) 
    182192            ! 
    183193            ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 
    184             zCdu = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / e3u_n(ji,jj,ikbu) 
    185             zCdv = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / e3v_n(ji,jj,ikbv) 
     194            zCdu = 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / e3u(ji,jj,ikbu,Kmm)    ! NB: Cdtop masked 
     195            zCdv = 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / e3v(ji,jj,ikbv,Kmm) 
    186196            ! 
    187197            pua(ji,jj,ikbu) = pua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * pub(ji,jj,ikbu) 
    188198            pva(ji,jj,ikbv) = pva(ji,jj,ikbv) + MAX(  zCdv , zm1_2dt  ) * pvb(ji,jj,ikbv) 
    189          END DO 
    190       END DO 
    191       ! 
    192       IF( ln_isfcav ) THEN        ! ocean cavities 
    193          DO jj = 2, jpjm1 
    194             DO ji = 2, jpim1 
    195                ikbu = miku(ji,jj)          ! first wet ocean u- & v-levels 
    196                ikbv = mikv(ji,jj) 
    197                ! 
    198                ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 
    199                zCdu = 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / e3u_n(ji,jj,ikbu)    ! NB: Cdtop masked 
    200                zCdv = 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / e3v_n(ji,jj,ikbv) 
    201                ! 
    202                pua(ji,jj,ikbu) = pua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * pub(ji,jj,ikbu) 
    203                pva(ji,jj,ikbv) = pva(ji,jj,ikbv) + MAX(  zCdv , zm1_2dt  ) * pvb(ji,jj,ikbv) 
    204            END DO 
    205          END DO 
     199         END_2D 
    206200      ENDIF 
    207201      ! 
     
    209203         ztrdu(:,:,:) = pua(:,:,:) - ztrdu(:,:,:) 
    210204         ztrdv(:,:,:) = pva(:,:,:) - ztrdv(:,:,:) 
    211          CALL trd_dyn( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_bfr, kt ) 
     205         CALL trd_dyn( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_bfr, kt, Kmm ) 
    212206         DEALLOCATE( ztrdu, ztrdv ) 
    213207      ENDIF 
    214208      !                                          ! print mean trends (used for debugging) 
    215       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pua, clinfo1=' bfr  - Ua: ', mask1=umask,               & 
    216          &                       tab3d_2=pva, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     209      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pua, clinfo1=' bfr  - Ua: ', mask1=umask,               & 
     210         &                                  tab3d_2=pva, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    217211      ! 
    218212   END SUBROUTINE zdf_drg_exp 
     
    236230      !                     !==  drag nature  ==! 
    237231      ! 
    238       REWIND( numnam_ref )                   ! Namelist namdrg in reference namelist 
    239232      READ  ( numnam_ref, namdrg, IOSTAT = ios, ERR = 901) 
    240233901   IF( ios /= 0 )   CALL ctl_nam( ios , 'namdrg in reference namelist' ) 
    241       REWIND( numnam_cfg )                   ! Namelist namdrg in configuration namelist 
    242234      READ  ( numnam_cfg, namdrg, IOSTAT = ios, ERR = 902 ) 
    243235902   IF( ios >  0 )   CALL ctl_nam( ios , 'namdrg in configuration namelist' ) 
     
    335327      !                          !==  read namlist  ==! 
    336328      ! 
    337       REWIND( numnam_ref )                   ! Namelist cl_namdrg in reference namelist 
    338329      IF(ll_top)   READ  ( numnam_ref, namdrg_top, IOSTAT = ios, ERR = 901) 
    339330      IF(ll_bot)   READ  ( numnam_ref, namdrg_bot, IOSTAT = ios, ERR = 901) 
    340331901   IF( ios /= 0 )   CALL ctl_nam( ios , TRIM(cl_namref) ) 
    341       REWIND( numnam_cfg )                   ! Namelist cd_namdrg in configuration namelist 
    342332      IF(ll_top)   READ  ( numnam_cfg, namdrg_top, IOSTAT = ios, ERR = 902 ) 
    343333      IF(ll_bot)   READ  ( numnam_cfg, namdrg_bot, IOSTAT = ios, ERR = 902 ) 
     
    431421            l_log_not_linssh = .FALSE.    !- don't update Cd at each time step 
    432422            ! 
    433             DO jj = 1, jpj                   ! pCd0 = mask (and boosted) logarithmic drag coef.  
    434                DO ji = 1, jpi 
    435                   zzz =  0.5_wp * e3t_0(ji,jj,k_mk(ji,jj)) 
    436                   zcd = (  vkarmn / LOG( zzz / rn_z0 )  )**2 
    437                   pCd0(ji,jj) = zmsk_boost(ji,jj) * MIN(  MAX( rn_Cd0 , zcd ) , rn_Cdmax  )  ! rn_Cd0 < Cd0 < rn_Cdmax 
    438                END DO 
    439             END DO 
     423            DO_2D_11_11 
     424               zzz =  0.5_wp * e3t_0(ji,jj,k_mk(ji,jj)) 
     425               zcd = (  vkarmn / LOG( zzz / rn_z0 )  )**2 
     426               pCd0(ji,jj) = zmsk_boost(ji,jj) * MIN(  MAX( rn_Cd0 , zcd ) , rn_Cdmax  )  ! rn_Cd0 < Cd0 < rn_Cdmax 
     427            END_2D 
    440428         ELSE                       !* Cd updated at each time-step ==> pCd0 = mask * boost 
    441429            IF(lwp) WRITE(numout,*) 
Note: See TracChangeset for help on using the changeset viewer.