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 1518 for trunk/NEMO/OPA_SRC/ZDF/zdftmx.F90 – NEMO

Ignore:
Timestamp:
2009-07-22T11:14:04+02:00 (15 years ago)
Author:
ctlod
Message:

optimisation finalisation, see ticket, #457

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/OPA_SRC/ZDF/zdftmx.F90

    r1496 r1518  
    2727   PRIVATE 
    2828 
    29    PUBLIC zdf_tmx          ! called in step module  
     29   PUBLIC   zdf_tmx    ! called in step module  
    3030 
    3131   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftmx = .TRUE.    !: tidal mixing flag 
    3232 
    33    !                                                    !!* Namelist  namtmx : tidal mixing * 
    34    REAL(wp) ::  rn_htmx  = 500.                          ! vertical decay scale for turbulence (meters) 
    35    REAL(wp) ::  rn_n2min = 1.e-8                         ! threshold of the Brunt-Vaisala frequency (s-1) 
    36    REAL(wp) ::  rn_tfe   = 1./3.                         ! tidal dissipation efficiency (St Laurent et al. 2002) 
    37    REAL(wp) ::  rn_tfe_itf  = 1.                         ! ITF tidal dissipation efficiency (St Laurent et al. 2002) 
    38    REAL(wp) ::  rn_me    = 0.2                           ! mixing efficiency (Osborn 1980) 
    39  
    40    REAL(wp), DIMENSION(jpi,jpj) ::   en_tmx              ! energy available for tidal mixing (W/m2) 
    41    REAL(wp), DIMENSION(jpi,jpj) ::   mask_itf            ! mask to use over Indonesian area 
    42  
    43    REAL(wp), DIMENSION(jpi,jpj,jpk) ::   az_tmx          ! coefficient used to evaluate the tidal induced Kz 
     33   !                                  !!* Namelist  namtmx : tidal mixing * 
     34   REAL(wp) ::  rn_htmx    = 500.      ! vertical decay scale for turbulence (meters) 
     35   REAL(wp) ::  rn_n2min   = 1.e-8     ! threshold of the Brunt-Vaisala frequency (s-1) 
     36   REAL(wp) ::  rn_tfe     = 1./3.     ! tidal dissipation efficiency (St Laurent et al. 2002) 
     37   REAL(wp) ::  rn_me      = 0.2       ! mixing efficiency (Osborn 1980) 
     38   LOGICAL  ::  ln_tmx_itf = .TRUE.    ! Indonesian Through Flow (ITF): Koch-Larrouy et al. (2007) parameterization 
     39   REAL(wp) ::  rn_tfe_itf = 1.        ! ITF tidal dissipation efficiency (St Laurent et al. 2002) 
     40 
     41   REAL(wp), DIMENSION(jpi,jpj)     ::   en_tmx     ! energy available for tidal mixing (W/m2) 
     42   REAL(wp), DIMENSION(jpi,jpj)     ::   mask_itf   ! mask to use over Indonesian area 
     43   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   az_tmx     ! coefficient used to evaluate the tidal induced Kz 
    4444 
    4545   !! * Substitutions 
     
    9797      IF( kt == nit000  )   CALL zdf_tmx_init      ! Initialization (first time-step only) 
    9898 
    99       !                     ! ----------------------- ! 
    100       !                     !  Standard tidal mixing  !  (compute av_tide) 
    101       !                     ! ----------------------- ! 
     99      !                                        ! ----------------------- ! 
     100      !                                        !  Standard tidal mixing  !  (compute av_tide) 
     101      !                                        ! ----------------------- ! 
    102102      !                             !* First estimation (with n2 bound by rn_n2min) bounded by 60 cm2/s 
    103103      av_tide(:,:,:) = MIN(  60.e-4, az_tmx(:,:,:) / MAX( rn_n2min, rn2(:,:,:) )  ) 
     
    133133      ENDIF 
    134134        
    135       !                     ! ----------------------- ! 
    136       CALL tmx_itf( kt )    !    ITF  tidal mixing    !  (update av_tide) 
    137       !                     ! ----------------------- ! 
    138  
    139       !                     ! ----------------------- ! 
    140       !                     !   Update  mixing coefs  !                           
    141       !                     ! ----------------------- ! 
     135      !                                        ! ----------------------- ! 
     136      IF( ln_tmx_itf )   CALL tmx_itf( kt )    !    ITF  tidal mixing    !  (update av_tide) 
     137      !                                        ! ----------------------- ! 
     138 
     139      !                                        ! ----------------------- ! 
     140      !                                        !   Update  mixing coefs  !                           
     141      !                                        ! ----------------------- ! 
    142142      DO jk = 2, jpkm1              !* update momentum & tracer diffusivity with tidal mixing 
    143143         avt(:,:,jk) = avt(:,:,jk) + av_tide(:,:,jk) 
     
    180180      INTEGER  ::   ji, jj, jk    ! dummy loop indices 
    181181      REAL(wp) ::   zcoef, ztpc   ! temporary scalar 
    182       REAL(wp), DIMENSION(jpi,jpj) ::   zkz   ! temporary 2D workspace 
    183       REAL(wp), DIMENSION(jpi,jpj) ::   zsum1 , zsum2 , zsum     
    184       REAL(wp), DIMENSION(jpi,jpj) ::   z1sum1, z1sum2, z1sum   
    185       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zempba_3d_1, zempba_3d_2   ! 3D workspace 
    186       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zempba_3d  , zdn2dz        !  -      - 
    187       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zavt_itf                   !  -      - 
     182      REAL(wp), DIMENSION(jpi,jpj)     ::   zkz                        ! 2D workspace 
     183      REAL(wp), DIMENSION(jpi,jpj)     ::   zsum1 , zsum2 , zsum       !  -      - 
     184      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zempba_3d_1, zempba_3d_2   ! 3D workspace 
     185      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zempba_3d  , zdn2dz        !  -      - 
     186      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zavt_itf                   !  -      - 
    188187      !!---------------------------------------------------------------------- 
    189188 
    190189      !                             ! compute the form function using N2 at each time step 
    191 !!gm better/faster coding  (not tmask as rn2 is already masked) 
    192 !      zempba_3d_1(:,:,jpk) = 0.e0 
    193 !      zempba_3d_2(:,:,jpk) = 0.e0 
    194 !      DO jk = 1, jpkm1              
    195 !         zdn2dz     (:,:,jk) = rn2(:,:,jk) - rn2(:,:,jk+1)           ! Vertical profile of dN2/dz 
     190      zempba_3d_1(:,:,jpk) = 0.e0 
     191      zempba_3d_2(:,:,jpk) = 0.e0 
     192      DO jk = 1, jpkm1              
     193         zdn2dz     (:,:,jk) = rn2(:,:,jk) - rn2(:,:,jk+1)           ! Vertical profile of dN2/dz 
    196194!CDIR NOVERRCHK 
    197 !         zempba_3d_1(:,:,jk) = SQRT(  MAX( O.e0, rn2(:,:,jk) )  )    !    -        -    of N 
    198 !         zempba_3d_2(:,:,jk) =        MAX( O.e0, rn2(:,:,jk) )       !    -        -    of N^2 
    199 !      END DO 
    200 !!gm replacing: 
    201       zempba_3d_1(:,:,:) = 0.e0 
    202       zempba_3d_2(:,:,:) = 0.e0 
    203     
    204       DO jk = 1, jpk-1              ! Vertical profile dN2/dz 
    205          DO jj = 1, jpj                
    206             DO ji = 1, jpi 
    207                zdn2dz(ji,jj,jk) =(rn2(ji,jj,jk) - rn2(ji,jj,jk+1) ) 
    208                IF( rn2(ji,jj,jk) > 0 ) THEN  
    209                   zempba_3d_1(ji,jj,jk)=SQRT(rn2(ji,jj,jk))*tmask(ji,jj,jk) 
    210                   zempba_3d_2(ji,jj,jk)=rn2(ji,jj,jk)*tmask(ji,jj,jk) 
    211                ENDIF 
    212             END DO 
    213          END DO 
    214       END DO 
    215 !!gm end 
    216  
    217 !!gm better/faster coding  (suppression of z1sum, z1sum1 and z1sum2) 
    218 !      zsum (:,:) = 0.e0 
    219 !      zsum1(:,:) = 0.e0 
    220 !      zsum2(:,:) = 0.e0 
    221 !      DO jk= 2, jpk 
    222 !         zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * fse3w(:,:,jk) 
    223 !         zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * fse3w(:,:,jk)                 
    224 !      END DO 
    225 !      DO jj = 1, jpj 
    226 !         DO ji = 1, jpi 
    227 !            IF( zsum1(ji,jj) /= 0.e0 )   zsum1(ji,jj) = 1.e0 / zsum1(ji,jj) 
    228 !            IF( zsum2(ji,jj) /= 0.e0 )   zsum2(ji,jj) = 1.e0 / zsum2(ji,jj)                 
    229 !         END DO 
    230 !      END DO 
    231 ! 
    232 !      DO jk= 1, jpk 
    233 !         DO jj = 1, jpj 
    234 !            DO ji = 1, jpi 
    235 !               zcoef = 0.5 - SIGN( 0.5, zdn2dz(ji,jj,jk) )       ! =0 if dN2/dz > 0, =1 otherwise  
    236 !               ztpc  = zempba_3d_1(ji,jj,jk) * z1sum1(ji,jj) *        zcoef     & 
    237 !                  &  + zempba_3d_2(ji,jj,jk) * z1sum2(ji,jj) * ( 1. - zcoef ) 
    238 !               ! 
    239 !               zempba_3d(ji,jj,jk) =               ztpc  
    240 !               zsum     (ji,jj)    = zsum(ji,jj) + ztpc * fse3w(ji,jj,jk) 
    241 !            END DO 
    242 !         END DO 
    243 !       END DO 
    244 !       DO jj = 1, jpj 
    245 !          DO ji = 1, jpi 
    246 !             IF( zsum(ji,jj) > 0.e0 )   zsum(ji,jj) = en_tmx(:,:) / zsum(ji,jj)                 
    247 !          END DO 
    248 !       END DO 
    249 ! 
    250 !      !                             ! first estimation bounded by 10 cm2/s (with n2 bounded by rn_n2min)  
    251 !      zcoef = rn_tfe_itf / ( rn_tfe * rau0 ) 
    252 !      DO jk = 1, jpk 
    253 !         zavt_itf(:,:,jk) = MIN(  10.e-4, zcoef * zsum(:,:) * zempba_3d(:,:,jk)   & 
    254 !            &                                      / MAX( rn_n2min, rn2(:,:,jk) ) * tmask(:,:,jk)  ) 
    255 !      END DO            
    256 !!gm replacing: 
    257  
    258       zsum1(:,:)=0.e0 
    259       zsum2(:,:)=0.e0 
     195         zempba_3d_1(:,:,jk) = SQRT(  MAX( 0.e0, rn2(:,:,jk) )  )    !    -        -    of N 
     196         zempba_3d_2(:,:,jk) =        MAX( 0.e0, rn2(:,:,jk) )       !    -        -    of N^2 
     197      END DO 
     198      ! 
     199      zsum (:,:) = 0.e0 
     200      zsum1(:,:) = 0.e0 
     201      zsum2(:,:) = 0.e0 
    260202      DO jk= 2, jpk 
    261          DO jj = 1, jpj 
    262             DO ji = 1, jpi 
    263                zsum1(ji,jj)=zsum1(ji,jj)+zempba_3d_1(ji,jj,jk)*fse3w(ji,jj,jk) 
    264                zsum2(ji,jj)=zsum2(ji,jj)+zempba_3d_2(ji,jj,jk)*fse3w(ji,jj,jk)                 
    265             END DO 
    266          END DO 
    267       END DO 
    268  
    269       z1sum1(:,:)=0.e0 
    270       z1sum2(:,:)=0.e0 
     203         zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * fse3w(:,:,jk) 
     204         zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * fse3w(:,:,jk)                 
     205      END DO 
    271206      DO jj = 1, jpj 
    272         DO ji = 1, jpi 
    273          IF (zsum1(ji,jj) /= 0) z1sum1(ji,jj)=1.0/zsum1(ji,jj) 
    274          IF (zsum2(ji,jj) /= 0) z1sum2(ji,jj)=1.0/zsum2(ji,jj)                 
    275         END DO 
    276       END DO 
    277  
    278       zsum (:,:)=0.e0 
    279       z1sum(:,:)=0.e0 
     207         DO ji = 1, jpi 
     208            IF( zsum1(ji,jj) /= 0.e0 )   zsum1(ji,jj) = 1.e0 / zsum1(ji,jj) 
     209            IF( zsum2(ji,jj) /= 0.e0 )   zsum2(ji,jj) = 1.e0 / zsum2(ji,jj)                 
     210         END DO 
     211      END DO 
     212 
    280213      DO jk= 1, jpk 
    281214         DO jj = 1, jpj 
    282215            DO ji = 1, jpi 
    283                zempba_3d(ji,jj,jk)=zempba_3d_1(ji,jj,jk)*z1sum1(ji,jj)*0.5*(1-SIGN(1.0,zdn2dz(ji,jj,jk)))       & 
    284                   &               +zempba_3d_2(ji,jj,jk)*z1sum2(ji,jj)*0.5*(1+SIGN(1.0,zdn2dz(ji,jj,jk))) 
    285                zsum(ji,jj)=zsum(ji,jj)+zempba_3d(ji,jj,jk)*fse3w(ji,jj,jk) 
     216               zcoef = 0.5 - SIGN( 0.5, zdn2dz(ji,jj,jk) )       ! =0 if dN2/dz > 0, =1 otherwise  
     217               ztpc  = zempba_3d_1(ji,jj,jk) * zsum1(ji,jj) *        zcoef     & 
     218                  &  + zempba_3d_2(ji,jj,jk) * zsum2(ji,jj) * ( 1. - zcoef ) 
     219               ! 
     220               zempba_3d(ji,jj,jk) =               ztpc  
     221               zsum     (ji,jj)    = zsum(ji,jj) + ztpc * fse3w(ji,jj,jk) 
    286222            END DO 
    287223         END DO 
    288224       END DO 
    289  
    290225       DO jj = 1, jpj 
    291226          DO ji = 1, jpi 
    292              IF( zsum(ji,jj) > 0 ) z1sum(ji,jj)=1.0/zsum(ji,jj)                 
     227             IF( zsum(ji,jj) > 0.e0 )   zsum(ji,jj) = 1.e0 / zsum(ji,jj)                 
    293228          END DO 
    294229       END DO 
     
    296231      !                             ! first estimation bounded by 10 cm2/s (with n2 bounded by rn_n2min)  
    297232      zcoef = rn_tfe_itf / ( rn_tfe * rau0 ) 
    298       DO jk= 1, jpk 
    299          zavt_itf(:,:,jk) = MIN(  10.e-4, zcoef * en_tmx(:,:) * z1sum(:,:) * zempba_3d(:,:,jk)   & 
    300             &                                   / MAX( rn_n2min, rn2(:,:,jk) ) * tmask(:,:,jk)  ) 
     233      DO jk = 1, jpk 
     234         zavt_itf(:,:,jk) = MIN(  10.e-4, zcoef * en_tmx(:,:) * zsum(:,:) * zempba_3d(:,:,jk)   & 
     235            &                                      / MAX( rn_n2min, rn2(:,:,jk) ) * tmask(:,:,jk)  ) 
    301236      END DO            
    302 !!gm end 
    303237 
    304238      zkz(:,:) = 0.e0               ! Associated potential energy consummed over the whole water column 
     
    375309      !!              Koch-Larrouy et al. 2007, GRL. 
    376310      !!---------------------------------------------------------------------- 
    377       INTEGER ::   ji, jj, jk                           ! dummy loop indices 
    378       INTEGER ::   inum                                 ! temporary logical unit 
    379       REAL(wp) ::   ztpc, ze_z                          ! total power consumption 
    380       REAL(wp), DIMENSION(jpi,jpj) ::  zem2, zek1       ! read M2 and K1 tidal energy 
    381       REAL(wp), DIMENSION(jpi,jpj) ::  zkz              ! total M2, K1 and S2 tidal energy 
    382       REAL(wp), DIMENSION(jpi,jpj) ::  zfact            ! used for vertical structure function 
    383       REAL(wp), DIMENSION(jpi,jpj) ::  zhdep                ! Ocean depth  
    384       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zpc          ! power consumption 
    385       !! 
    386       NAMELIST/namtmx/ rn_htmx, rn_n2min, rn_tfe, rn_tfe_itf, rn_me 
     311      INTEGER ::   ji, jj, jk    ! dummy loop indices 
     312      INTEGER ::   inum          ! temporary logical unit 
     313      REAL(wp) ::   ztpc, ze_z   ! total power consumption 
     314      REAL(wp), DIMENSION(jpi,jpj) ::  zem2, zek1   ! read M2 and K1 tidal energy 
     315      REAL(wp), DIMENSION(jpi,jpj) ::  zkz          ! total M2, K1 and S2 tidal energy 
     316      REAL(wp), DIMENSION(jpi,jpj) ::  zfact        ! used for vertical structure function 
     317      REAL(wp), DIMENSION(jpi,jpj) ::  zhdep        ! Ocean depth  
     318      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zpc      ! power consumption 
     319      !! 
     320      NAMELIST/namtmx/ rn_htmx, rn_n2min, rn_tfe, rn_me, ln_tmx_itf, rn_tfe_itf 
    387321      !!---------------------------------------------------------------------- 
    388322 
     
    404338         WRITE(numout,*) '             Brunt-Vaisala frequency threshold     = ', rn_n2min 
    405339         WRITE(numout,*) '             Tidal dissipation efficiency          = ', rn_tfe 
     340         WRITE(numout,*) '             Mixing efficiency                     = ', rn_me 
     341         WRITE(numout,*) '             ITF specific parameterisation         = ', ln_tmx_itf 
    406342         WRITE(numout,*) '             ITF tidal dissipation efficiency      = ', rn_tfe_itf 
    407          WRITE(numout,*) '             Mixing efficiency                     = ', rn_me 
    408343         WRITE(numout,*) 
    409344      ENDIF 
    410345 
    411       ! read mask kz banda 
    412       CALL iom_open('mask_itf',inum) 
    413       CALL iom_get (inum, jpdom_data, 'tmaskitf',mask_itf,1) !  
    414       CALL iom_close(inum) 
     346      IF( ln_tmx_itf ) THEN   ! read the Indonesian Through Flow mask 
     347         CALL iom_open('mask_itf',inum) 
     348         CALL iom_get (inum, jpdom_data, 'tmaskitf',mask_itf,1) !  
     349         CALL iom_close(inum) 
     350      ENDIF 
    415351 
    416352      ! read M2 tidal energy flux : W/m2  ( zem2 < 0 ) 
Note: See TracChangeset for help on using the changeset viewer.