Ignore:
Timestamp:
2011-09-28T10:18:52+02:00 (9 years ago)
Author:
hliu
Message:

1) added the semi-implicit bottom friction code (3 in DYN, 1 in ZDF), 2) added par_AMM7/12.h90, 3) added two arch files for NOCL mobius and ubuntu linux

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r2715 r2872  
    2020   USE in_out_manager  ! I/O manager 
    2121   USE lib_mpp         ! MPP library 
     22   USE zdfbfr          ! bottom friction setup 
    2223 
    2324   IMPLICIT NONE 
     
    6162      REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step 
    6263      !! 
    63       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    64       REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars 
     64      INTEGER  ::   ji, jj, jk     ! dummy loop indices 
     65      INTEGER  ::   ikbum1, ikbvm1 ! local variable 
     66      REAL(wp) ::   z1_p2dt, z2dtf, zcoef, zzwi, zzws, zrhs   ! local scalars 
     67 
     68      !! * Local variables for implicit bottom friction.    H. Liu 
     69      REAL(wp) ::   zbfru, zbfrv  
     70      REAL(wp) ::   bfr_imp = 1._wp 
     71      !!---------------------------------------------------------------------- 
    6572      !!---------------------------------------------------------------------- 
    6673 
     
    7380         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator' 
    7481         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
     82         IF(.NOT.ln_bfrimp) bfr_imp = 0._wp 
    7583      ENDIF 
    7684 
     
    8088 
    8189      ! 1. Vertical diffusion on u 
     90 
     91      ! Vertical diffusion on u&v 
    8292      ! --------------------------- 
    8393      ! Matrix and second member construction 
    84       ! bottom boundary condition: both zwi and zws must be masked as avmu can take 
    85       ! non zero value at the ocean bottom depending on the bottom friction 
    86       ! used but the bottom velocities have already been updated with the bottom 
    87       ! friction velocity in dyn_bfr using values from the previous timestep. There 
    88       ! is no need to include these in the implicit calculation. 
    89       ! 
    90       DO jk = 1, jpkm1        ! Matrix 
    91          DO jj = 2, jpjm1  
    92             DO ji = fs_2, fs_jpim1   ! vector opt. 
     94      !! bottom boundary condition: both zwi and zws must be masked as avmu can take 
     95      !! non zero value at the ocean bottom depending on the bottom friction 
     96      !! used but the bottom velocities have already been updated with the bottom 
     97      !! friction velocity in dyn_bfr using values from the previous timestep. There 
     98      !! is no need to include these in the implicit calculation. 
     99 
     100      ! The code has been modified here to implicitly implement bottom 
     101      ! friction: u(v)mask is not necessary here anymore.  
     102      ! H. Liu, April 2010. 
     103 
     104      ! 1. Vertical diffusion on u 
     105      DO jj = 2, jpjm1  
     106         DO ji = fs_2, fs_jpim1   ! vector opt. 
     107            ikbum1 = mbku(ji,jj) 
     108            ! Apply stability criteria on absolute value  : Min abs(bfr) => Max (bfr) 
     109!               zbfru = MAX( bfrua(ji,jj), fse3u(ji,jj,ikbum1)*zinv ) 
     110               zbfru = bfrua(ji,jj) 
     111 
     112            DO jk = 1, ikbum1 
    93113               zcoef = - p2dt / fse3u(ji,jj,jk) 
    94                zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  ) 
    95                zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk) 
    96                zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
    97                zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1) 
    98                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 
    99             END DO 
    100          END DO 
    101       END DO 
    102       DO jj = 2, jpjm1        ! Surface boudary conditions 
    103          DO ji = fs_2, fs_jpim1   ! vector opt. 
     114               zwi(ji,jj,jk) = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  ) 
     115               zws(ji,jj,jk) = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
     116               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk) 
     117            END DO 
     118 
     119      ! Surface boudary conditions 
    104120            zwi(ji,jj,1) = 0._wp 
    105121            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
    106          END DO 
    107       END DO 
     122 
     123      ! Bottom boudary conditions  ! H. Liu, May, 2010 
     124!            !commented out to be consisent with v3.2, h.liu 
     125!            z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 2.0 * bfr_imp 
     126            z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 1.0_wp * bfr_imp 
     127            zws(ji,jj,ikbum1) = 0._wp 
     128            zwd(ji,jj,ikbum1) = 1._wp - zwi(ji,jj,ikbum1) - z2dtf  
    108129 
    109130      ! Matrix inversion starting from the first level 
     
    121142      !   The solution (the after velocity) is in ua 
    122143      !----------------------------------------------------------------------- 
    123       ! 
    124       DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    125          DO jj = 2, jpjm1    
    126             DO ji = fs_2, fs_jpim1   ! vector opt. 
     144 
     145      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k) 
     146            DO jk = 2, ikbum1 
    127147               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
    128148            END DO 
    129          END DO 
    130       END DO 
    131       ! 
    132       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
    133          DO ji = fs_2, fs_jpim1   ! vector opt. 
    134             ua(ji,jj,1) = ub(ji,jj,1) + p2dt * (  ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    135                &                                                       / ( fse3u(ji,jj,1) * rau0       )  ) 
    136          END DO 
    137       END DO 
    138       DO jk = 2, jpkm1 
    139          DO jj = 2, jpjm1    
    140             DO ji = fs_2, fs_jpim1   ! vector opt. 
     149 
     150      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1 
     151            z2dtf = 0.5_wp * p2dt / ( fse3u(ji,jj,1) * rau0 ) 
     152            ua(ji,jj,1) = ub(ji,jj,1) + p2dt * ua(ji,jj,1) + z2dtf * (utau_b(ji,jj) + utau(ji,jj)) 
     153           DO jk = 2, ikbum1    
    141154               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side 
    142155               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 
    143156            END DO 
    144          END DO 
    145       END DO 
    146       ! 
    147       DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  == 
    148          DO ji = fs_2, fs_jpim1   ! vector opt. 
    149             ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
    150          END DO 
    151       END DO 
    152       DO jk = jpk-2, 1, -1 
    153          DO jj = 2, jpjm1    
    154             DO ji = fs_2, fs_jpim1   ! vector opt. 
    155                ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
     157 
     158 
     159      ! thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk 
     160            ua(ji,jj,ikbum1) = ua(ji,jj,ikbum1) / zwd(ji,jj,ikbum1) 
     161            DO jk = ikbum1-1, 1, -1 
     162               ua(ji,jj,jk) =( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
    156163            END DO 
    157164         END DO 
     
    170177      ! 2. Vertical diffusion on v 
    171178      ! --------------------------- 
    172       ! Matrix and second member construction 
    173       ! bottom boundary condition: both zwi and zws must be masked as avmv can take 
    174       ! non zero value at the ocean bottom depending on the bottom friction 
    175       ! used but the bottom velocities have already been updated with the bottom 
    176       ! friction velocity in dyn_bfr using values from the previous timestep. There 
    177       ! is no need to include these in the implicit calculation. 
    178       ! 
    179       DO jk = 1, jpkm1        ! Matrix 
     179 
     180      DO ji = fs_2, fs_jpim1   ! vector opt. 
    180181         DO jj = 2, jpjm1    
    181             DO ji = fs_2, fs_jpim1   ! vector opt. 
     182            ikbvm1 = mbkv(ji,jj) 
     183            ! Apply stability criteria on absolute value  : Min abs(bfr) => Max (bfr) 
     184!               zbfrv = MAX( bfrva(ji,jj), fse3v(ji,jj,ikbvm1)*zinv ) 
     185               zbfrv = bfrva(ji,jj) 
     186 
     187            DO jk = 1, ikbvm1 
    182188               zcoef = -p2dt / fse3v(ji,jj,jk) 
    183                zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  ) 
    184                zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk) 
    185                zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 
    186                zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1) 
    187                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 
    188             END DO 
    189          END DO 
    190       END DO 
    191       DO jj = 2, jpjm1        ! Surface boudary conditions 
    192          DO ji = fs_2, fs_jpim1   ! vector opt. 
     189               zwi(ji,jj,jk) = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  ) 
     190               zws(ji,jj,jk) = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 
     191               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk) 
     192            END DO 
     193 
     194      ! Surface boudary conditions 
    193195            zwi(ji,jj,1) = 0._wp 
    194196            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
    195          END DO 
    196       END DO 
     197 
     198      ! Bottom boudary conditions   
     199      ! Bottom boudary conditions  ! H. Liu, May, 2010 
     200!            !commented out to be consisent with v3.2, h.liu 
     201!            z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 2.0 * bfr_imp 
     202            z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 1.0_wp * bfr_imp 
     203            zws(ji,jj,ikbvm1) = 0._wp 
     204            zwd(ji,jj,ikbvm1) = 1._wp - zwi(ji,jj,ikbvm1) - z2dtf  
    197205 
    198206      ! Matrix inversion 
     
    206214      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk ) 
    207215      ! 
    208       !   m is decomposed in the product of an upper and lower triangular matrix 
     216      !   m is decomposed in the product of an upper and lower triangular 
     217      !   matrix 
    209218      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 
    210219      !   The solution (after velocity) is in 2d array va 
    211220      !----------------------------------------------------------------------- 
    212       ! 
    213       DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    214          DO jj = 2, jpjm1    
    215             DO ji = fs_2, fs_jpim1   ! vector opt. 
     221 
     222      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k) 
     223            DO jk = 2, ikbvm1 
    216224               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
    217225            END DO 
    218          END DO 
    219       END DO 
    220       ! 
    221       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
    222          DO ji = fs_2, fs_jpim1   ! vector opt. 
    223             va(ji,jj,1) = vb(ji,jj,1) + p2dt * (  va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    224                &                                                       / ( fse3v(ji,jj,1) * rau0       )  ) 
    225          END DO 
    226       END DO 
    227       DO jk = 2, jpkm1 
    228          DO jj = 2, jpjm1 
    229             DO ji = fs_2, fs_jpim1   ! vector opt. 
     226 
     227      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1 
     228            z2dtf = 0.5_wp * p2dt / ( fse3v(ji,jj,1)*rau0 ) 
     229            va(ji,jj,1) = vb(ji,jj,1) + p2dt * va(ji,jj,1) + z2dtf * (vtau_b(ji,jj) + vtau(ji,jj)) 
     230            DO jk = 2, ikbvm1 
    230231               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side 
    231232               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 
    232233            END DO 
    233          END DO 
    234       END DO 
    235       ! 
    236       DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  == 
    237          DO ji = fs_2, fs_jpim1   ! vector opt. 
    238             va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
    239          END DO 
    240       END DO 
    241       DO jk = jpk-2, 1, -1 
    242          DO jj = 2, jpjm1    
    243             DO ji = fs_2, fs_jpim1   ! vector opt. 
    244                va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
    245             END DO 
     234 
     235      ! thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk 
     236            va(ji,jj,ikbvm1) = va(ji,jj,ikbvm1) / zwd(ji,jj,ikbvm1) 
     237 
     238            DO jk = ikbvm1-1, 1, -1 
     239               va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
     240            END DO 
     241 
    246242         END DO 
    247243      END DO 
     
    258254      IF( wrk_not_released(3, 3) )   CALL ctl_stop('dyn_zdf_imp: failed to release workspace array') 
    259255      ! 
     256 
    260257   END SUBROUTINE dyn_zdf_imp 
    261258 
Note: See TracChangeset for help on using the changeset viewer.