Changeset 3259


Ignore:
Timestamp:
2012-01-13T16:24:15+01:00 (9 years ago)
Author:
hliu
Message:

re-implement the semi-implicit bottom friction in DYNZDF_IMP.F90 to optimize the performance of this piece of code, Thanks for Italo's profiling results

Location:
branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90

    r3231 r3259  
    55   !!============================================================================== 
    66   !! History :  3.2  ! 2008-11  (A. C. Coward)  Original code 
     7   !!            3.4  ! 2011-09  (H. Liu) Make it consistent with semi-implicit 
     8   !!                            Bottom friction (ln_bfrimp = .true.)  
    79   !!---------------------------------------------------------------------- 
    810 
     
    5759      IF( .NOT.ln_bfrimp) THEN     ! only for explicit bottom friction form 
    5860                                    ! implicit bfr is implemented in dynzdf_imp 
    59                                     ! H. Liu,  Sept. 2011 
    6061 
    6162        zm1_2dt = - 1._wp / ( 2._wp * rdt ) 
  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r3231 r3259  
    88   !!   NEMO     0.5  !  2002-08  (G. Madec)  F90: Free form and module 
    99   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps 
     10   !!            3.4  !  2012-01  (H. Liu) Semi-implicit bottom friction 
    1011   !!---------------------------------------------------------------------- 
    1112 
     
    2021   USE in_out_manager  ! I/O manager 
    2122   USE lib_mpp         ! MPP library 
    22    USE zdfbfr          ! bottom friction setup 
     23   USE zdfbfr          ! Bottom friction setup 
    2324   USE wrk_nemo        ! Memory Allocation 
    2425   USE timing          ! Timing 
     
    6061      REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step 
    6162      !! 
    62       INTEGER  ::   ji, jj, jk     ! dummy loop indices 
    63       INTEGER  ::   ikbum1, ikbvm1 ! local variable 
    64       REAL(wp) ::   z1_p2dt, z2dtf, zcoef, zzwi, zzws, zrhs ! local scalars 
    65  
    66       !! * Local variables for implicit bottom friction.    H. Liu 
    67       REAL(wp) ::   zbfru, zbfrv  
    68       REAL(wp) ::   zbfr_imp = 0._wp                        ! toggle (SAVE'd by assignment)  
     63      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     64      INTEGER  ::   ikbu, ikbv   ! local integers 
     65      REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars 
     66      !!---------------------------------------------------------------------- 
     67 
    6968      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwd, zws 
    7069      !!---------------------------------------------------------------------- 
     
    7877         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator' 
    7978         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
    80          IF(ln_bfrimp) zbfr_imp = 1._wp 
    8179      ENDIF 
    8280 
     
    8583      z1_p2dt = 1._wp / p2dt      ! inverse of the timestep 
    8684 
    87       ! 1. Vertical diffusion on u 
    88  
    89       ! Vertical diffusion on u&v 
     85      ! 1. Apply semi-implicit bottom friction 
     86      ! -------------------------------------- 
     87      ! Only needed for semi-implicit bottom friction setup. The explicit 
     88      ! bottom friction has been included in "u(v)a" which act as the R.H.S 
     89      ! column vector of the tri-diagonal matrix equation 
     90      ! 
     91      IF( ln_bfrimp ) THEN 
     92# if defined key_vectopt_loop 
     93      DO jj = 1, 1 
     94         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
     95# else 
     96      DO jj = 2, jpjm1 
     97         DO ji = 2, jpim1 
     98# endif 
     99            ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
     100            ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
     101            avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)  
     102            avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1) 
     103         END DO 
     104      END DO 
     105      ENDIF 
     106 
     107      ! 2. Vertical diffusion on u 
    90108      ! --------------------------- 
    91109      ! Matrix and second member construction 
    92       !! bottom boundary condition: both zwi and zws must be masked as avmu can take 
    93       !! non zero value at the ocean bottom depending on the bottom friction 
    94       !! used but the bottom velocities have already been updated with the bottom 
    95       !! friction velocity in dyn_bfr using values from the previous timestep. There 
    96       !! is no need to include these in the implicit calculation. 
    97  
    98       ! The code has been modified here to implicitly implement bottom 
    99       ! friction: u(v)mask is not necessary here anymore.  
    100       ! H. Liu, April 2010. 
    101  
    102       ! 1. Vertical diffusion on u 
    103       DO jj = 2, jpjm1  
    104          DO ji = fs_2, fs_jpim1   ! vector opt. 
    105             ikbum1 = mbku(ji,jj) 
    106                zbfru = bfrua(ji,jj) 
    107  
    108             DO jk = 1, ikbum1 
     110      ! bottom boundary condition: both zwi and zws must be masked as avmu can take 
     111      ! non zero value at the ocean bottom depending on the bottom friction used. 
     112      ! 
     113      DO jk = 1, jpkm1        ! Matrix 
     114         DO jj = 2, jpjm1  
     115            DO ji = fs_2, fs_jpim1   ! vector opt. 
    109116               zcoef = - p2dt / fse3u(ji,jj,jk) 
    110                zwi(ji,jj,jk) = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  ) 
    111                zws(ji,jj,jk) = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
    112                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk) 
    113             END DO 
    114  
    115       ! Surface boundary conditions 
     117               zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  ) 
     118               zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk) 
     119               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
     120               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1) 
     121               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 
     122            END DO 
     123         END DO 
     124      END DO 
     125      DO jj = 2, jpjm1        ! Surface boudary conditions 
     126         DO ji = fs_2, fs_jpim1   ! vector opt. 
    116127            zwi(ji,jj,1) = 0._wp 
    117128            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
    118  
    119       ! Bottom boundary conditions  ! H. Liu, May, 2010 
    120 !           !commented out to be consistent with v3.2, h.liu 
    121 !           z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 2.0_wp * zbfr_imp 
    122             z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 1.0_wp * zbfr_imp 
    123             zws(ji,jj,ikbum1) = 0._wp 
    124             zwd(ji,jj,ikbum1) = 1._wp - zwi(ji,jj,ikbum1) - z2dtf  
     129         END DO 
     130      END DO 
    125131 
    126132      ! Matrix inversion starting from the first level 
     
    138144      !   The solution (the after velocity) is in ua 
    139145      !----------------------------------------------------------------------- 
    140  
    141       ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k) 
    142             DO jk = 2, ikbum1 
     146      ! 
     147      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
     148         DO jj = 2, jpjm1    
     149            DO ji = fs_2, fs_jpim1   ! vector opt. 
    143150               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
    144151            END DO 
    145  
    146       ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1 
    147             z2dtf = 0.5_wp * p2dt / ( fse3u(ji,jj,1) * rau0 ) 
    148             ua(ji,jj,1) = ub(ji,jj,1) + p2dt * ua(ji,jj,1) + z2dtf * (utau_b(ji,jj) + utau(ji,jj)) 
    149            DO jk = 2, ikbum1    
     152         END DO 
     153      END DO 
     154      ! 
     155      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
     156         DO ji = fs_2, fs_jpim1   ! vector opt. 
     157            ua(ji,jj,1) = ub(ji,jj,1) + p2dt * (  ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
     158               &                                                       / ( fse3u(ji,jj,1) * rau0       )  ) 
     159         END DO 
     160      END DO 
     161      DO jk = 2, jpkm1 
     162         DO jj = 2, jpjm1    
     163            DO ji = fs_2, fs_jpim1   ! vector opt. 
    150164               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side 
    151165               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 
    152166            END DO 
    153  
    154  
    155       ! third recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk 
    156             ua(ji,jj,ikbum1) = ua(ji,jj,ikbum1) / zwd(ji,jj,ikbum1) 
    157             DO jk = ikbum1-1, 1, -1 
    158                ua(ji,jj,jk) =( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
     167         END DO 
     168      END DO 
     169      ! 
     170      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  == 
     171         DO ji = fs_2, fs_jpim1   ! vector opt. 
     172            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
     173         END DO 
     174      END DO 
     175      DO jk = jpk-2, 1, -1 
     176         DO jj = 2, jpjm1    
     177            DO ji = fs_2, fs_jpim1   ! vector opt. 
     178               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
    159179            END DO 
    160180         END DO 
     
    171191 
    172192 
    173       ! 2. Vertical diffusion on v 
     193      ! 3. Vertical diffusion on v 
    174194      ! --------------------------- 
    175  
    176       DO ji = fs_2, fs_jpim1   ! vector opt. 
    177          DO jj = 2, jpjm1    
    178             ikbvm1 = mbkv(ji,jj) 
    179                zbfrv = bfrva(ji,jj) 
    180  
    181             DO jk = 1, ikbvm1 
     195      ! Matrix and second member construction 
     196      ! bottom boundary condition: both zwi and zws must be masked as avmv can take 
     197      ! non zero value at the ocean bottom depending on the bottom friction used 
     198      ! 
     199      DO jk = 1, jpkm1        ! Matrix 
     200         DO jj = 2, jpjm1    
     201            DO ji = fs_2, fs_jpim1   ! vector opt. 
    182202               zcoef = -p2dt / fse3v(ji,jj,jk) 
    183                zwi(ji,jj,jk) = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  ) 
    184                zws(ji,jj,jk) = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 
    185                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk) 
    186             END DO 
    187  
    188       ! Surface boundary conditions 
     203               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  ) 
     204               zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk) 
     205               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 
     206               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1) 
     207               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 
     208            END DO 
     209         END DO 
     210      END DO 
     211      DO jj = 2, jpjm1        ! Surface boudary conditions 
     212         DO ji = fs_2, fs_jpim1   ! vector opt. 
    189213            zwi(ji,jj,1) = 0._wp 
    190214            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
    191  
    192       ! Bottom boundary conditions  ! H. Liu, May, 2010 
    193 !           !commented out to be consistent with v3.2, h.liu 
    194 !           z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 2.0_wp * zbfr_imp 
    195             z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 1.0_wp * zbfr_imp 
    196             zws(ji,jj,ikbvm1) = 0._wp 
    197             zwd(ji,jj,ikbvm1) = 1._wp - zwi(ji,jj,ikbvm1) - z2dtf  
     215         END DO 
     216      END DO 
    198217 
    199218      ! Matrix inversion 
     
    207226      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk ) 
    208227      ! 
    209       !   m is decomposed in the product of an upper and lower triangular 
    210       !   matrix 
     228      !   m is decomposed in the product of an upper and lower triangular matrix 
    211229      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 
    212230      !   The solution (after velocity) is in 2d array va 
    213231      !----------------------------------------------------------------------- 
    214  
    215       ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k) 
    216             DO jk = 2, ikbvm1 
     232      ! 
     233      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
     234         DO jj = 2, jpjm1    
     235            DO ji = fs_2, fs_jpim1   ! vector opt. 
    217236               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
    218237            END DO 
    219  
    220       ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1 
    221             z2dtf = 0.5_wp * p2dt / ( fse3v(ji,jj,1)*rau0 ) 
    222             va(ji,jj,1) = vb(ji,jj,1) + p2dt * va(ji,jj,1) + z2dtf * (vtau_b(ji,jj) + vtau(ji,jj)) 
    223             DO jk = 2, ikbvm1 
     238         END DO 
     239      END DO 
     240      ! 
     241      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
     242         DO ji = fs_2, fs_jpim1   ! vector opt. 
     243            va(ji,jj,1) = vb(ji,jj,1) + p2dt * (  va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
     244               &                                                       / ( fse3v(ji,jj,1) * rau0       )  ) 
     245         END DO 
     246      END DO 
     247      DO jk = 2, jpkm1 
     248         DO jj = 2, jpjm1 
     249            DO ji = fs_2, fs_jpim1   ! vector opt. 
    224250               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side 
    225251               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 
    226252            END DO 
    227  
    228       ! third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk 
    229             va(ji,jj,ikbvm1) = va(ji,jj,ikbvm1) / zwd(ji,jj,ikbvm1) 
    230  
    231             DO jk = ikbvm1-1, 1, -1 
    232                va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
    233             END DO 
    234  
     253         END DO 
     254      END DO 
     255      ! 
     256      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  == 
     257         DO ji = fs_2, fs_jpim1   ! vector opt. 
     258            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
     259         END DO 
     260      END DO 
     261      DO jk = jpk-2, 1, -1 
     262         DO jj = 2, jpjm1    
     263            DO ji = fs_2, fs_jpim1   ! vector opt. 
     264               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
     265            END DO 
    235266         END DO 
    236267      END DO 
  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

    r3229 r3259  
    123123            &                       tab2d_2=bfrva, clinfo2=       ' v: ', mask2=vmask,ovlap=1 ) 
    124124      ENDIF 
     125 
    125126      ! 
    126127      IF( nn_timing == 1 )  CALL timing_stop('zdf_bfr') 
     
    170171            WRITE(numout,*) 'Implicit bottom friction can only be used when ln_zdfexp=.false.' 
    171172            WRITE(numout,*) '         but you set: ln_bfrimp=.true. and ln_zdfexp=.true.!!!!' 
    172             WRITE(ctmp1,*)  '         bad ln_bfrimp value = .true.'  
     173            WRITE(ctmp1,*)  '         set either "ln_zdfexp = .false" or "ln_bfrimp = .false." 
    173174            CALL ctl_stop( ctmp1 ) 
    174175         END IF 
     
    272273         CALL mpp_max( zmaxbfr ) 
    273274      ENDIF 
     275      IF( .NOT.ln_bfrimp) THEN 
    274276      IF( lwp .AND. ictu + ictv > 0 ) THEN 
    275277         WRITE(numout,*) ' Bottom friction stability check failed at ', ictu, ' U-points ' 
     
    278280         WRITE(numout,*) ' Bottom friction coefficient will be reduced where necessary' 
    279281      ENDIF 
     282      ENDIF 
    280283      ! 
    281284      IF( nn_timing == 1 )  CALL timing_stop('zdf_bfr_init') 
Note: See TracChangeset for help on using the changeset viewer.