Changeset 12822


Ignore:
Timestamp:
2020-04-28T11:10:38+02:00 (5 months ago)
Author:
gm
Message:

symmetric sterss tensor and half cell modifications (wet point only, ghost cells)

Location:
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE
Files:
1 added
1 edited
3 copied

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domain.F90

    r12614 r12822  
    4444   USE lbclnk         ! ocean lateral boundary condition (or mpp link) 
    4545   USE lib_mpp        ! distributed memory computing library 
    46  
     46!!an45 
     47!   USE usrdef_nam, ONLY : ln_45machin 
     48   ! 
    4749   IMPLICIT NONE 
    4850   PRIVATE 
     
    5052   PUBLIC   dom_init     ! called by nemogcm.F90 
    5153   PUBLIC   domain_cfg   ! called by nemogcm.F90 
    52  
     54#  include "do_loop_substitute.h90" 
    5355   !!------------------------------------------------------------------------- 
    5456   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    8183      INTEGER , DIMENSION(jpi,jpj) ::   ik_top , ik_bot       ! top and bottom ocean level 
    8284      REAL(wp), DIMENSION(jpi,jpj) ::   z1_hu_0, z1_hv_0 
     85      REAL(wp)::   zcoeff                                     ! local real 
     86 
    8387      !!---------------------------------------------------------------------- 
    8488      ! 
     
    154158         hf_0(:,:) = hf_0(:,:) + e3f_0(:,:,jk) * ssfmask(:,:)     ! CAUTION : only valid in SWE, not with bathymetry 
    155159      END DO 
     160      ! 
     161!!anhf hf_0 = mean(ht_0*tmask) so hf = mimj( ht0 + ssht) 
     162! ne pas combiner avec an45 tout de suite 
     163!      DO_2D_10_10 
     164!         hf_0(ji,jj) = 0.25_wp * (   ht_0(ji,jj+1) * tmask(ji,jj+1,1) + ht_0(ji+1,jj+1) * tmask(ji+1,jj+1,1)   & 
     165!            &                      + ht_0(ji,jj  ) * tmask(ji,jj  ,1) + ht_0(ji+1,jj  ) * tmask(ji+1,jj  ,1)   ) 
     166!      END_2D 
     167!      CALL lbc_lnk( 'domain', hf_0, 'F', 1. )      ! Lateral boundary conditions 
     168!!anhf 
    156169      !                                 ! Inverse of reference ocean thickness 
    157170      r1_ht_0(:,:) =  ssmask(:,:) / ( ht_0(:,:) + 1._wp -  ssmask(:,:) ) 
     
    159172      r1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) ) 
    160173      r1_hf_0(:,:) = ssfmask(:,:) / ( hf_0(:,:) + 1._wp - ssfmask(:,:) ) 
    161        
     174      ! 
     175!!an45 Ligne de cote a 45deg : e1e2t *= ( mi(umask) + mj(vmask) ) /2 
     176!!                             idem pour e1e2f 
     177!      DO_2D_10_10 
     178!      zcoeff = 0.25_wp * (   umask(ji,jj+1,1) + umask(ji+1,jj+1,1)   & 
     179!         &                 + vmask(ji,jj  ,1) + vmask(ji+1,jj  ,1)   ) 
     180!         IF ( zcoeff /= 0._wp )   THEN 
     181!               e1e2t(ji,jj) =    e1e2t(ji,jj) * zcoeff 
     182!            r1_e1e2t(ji,jj) = r1_e1e2t(ji,jj) / zcoeff 
     183!         ENDIF 
     184!      END_2D 
     185!      WRITE(numout,*) '   an45 half T cell e1e2t ' 
     186!      zcoeff = 0.25_wp * (   umask(ji,jj+1,1) + umask(ji+1,jj+1,1)   & 
     187!         &                 + vmask(ji,jj  ,1) + vmask(ji+1,jj  ,1)   ) 
     188!         IF ( zcoeff /= 0._wp )   THEN 
     189!               e1e2t(ji,jj) =    e1e2t(ji,jj) * zcoeff 
     190!            r1_e1e2t(ji,jj) = r1_e1e2t(ji,jj) / zcoeff 
     191!!an45 
    162192      !           !==  time varying part of coordinate system  ==! 
    163193      ! 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynadv.F90

    r12529 r12822  
    114114         WRITE(numout,*) '      linear dynamics : no momentum advection          ln_dynadv_OFF  = ', ln_dynadv_OFF 
    115115         WRITE(numout,*) '      Vector form: 2nd order centered scheme           ln_dynadv_vec  = ', ln_dynadv_vec 
    116          WRITE(numout,*) '         with Hollingsworth scheme (=1) or not (=0)       nn_dynkeg   = ', nn_dynkeg 
     116!!an45 
     117         WRITE(numout,*) '         with Hollingsworth scheme (=1) or not (=0,2)      nn_dynkeg  = ', nn_dynkeg 
     118!!an45 
    117119         WRITE(numout,*) '      flux form: 2nd order centred scheme              ln_dynadv_cen2 = ', ln_dynadv_cen2 
    118120         WRITE(numout,*) '                 3rd order UBS scheme                  ln_dynadv_ubs  = ', ln_dynadv_ubs 
     
    126128 
    127129      IF( ioptio /= 1 )   CALL ctl_stop( 'choose ONE and only ONE advection scheme' ) 
    128       IF( nn_dynkeg /= nkeg_C2 .AND. nn_dynkeg /= nkeg_HW )   CALL ctl_stop( 'KEG scheme wrong value of nn_dynkeg' ) 
    129  
     130!!an45 
     131      IF( nn_dynkeg /= nkeg_C2 .AND. nn_dynkeg /= nkeg_C2_wpo .AND. nn_dynkeg /= nkeg_HW )   &  
     132         &   CALL ctl_stop( 'KEG scheme wrong value of nn_dynkeg' ) 
     133!!an45 
    130134 
    131135      IF(lwp) THEN                    ! Print the choice 
     
    134138         CASE( np_LIN_dyn )   ;   WRITE(numout,*) '   ==>>>   linear dynamics : no momentum advection used' 
    135139         CASE( np_VEC_c2  )   ;   WRITE(numout,*) '   ==>>>   vector form : keg + zad + vor is used'  
    136             IF( nn_dynkeg == nkeg_C2  )   WRITE(numout,*) '              with Centered standard keg scheme' 
    137             IF( nn_dynkeg == nkeg_HW  )   WRITE(numout,*) '              with Hollingsworth keg scheme' 
     140            IF( nn_dynkeg == nkeg_C2     )   WRITE(numout,*) '              with Centered standard keg scheme' 
     141!!an45 
     142            IF( nn_dynkeg == nkeg_C2_wpo )   WRITE(numout,*) '              with Centered standard keg scheme (wet point only)' 
     143!!an45 
     144            IF( nn_dynkeg == nkeg_HW     )   WRITE(numout,*) '              with Hollingsworth keg scheme' 
    138145         CASE( np_FLX_c2  )   ;   WRITE(numout,*) '   ==>>>   flux form   : 2nd order scheme is used' 
    139146         CASE( np_FLX_ubs )   ;   WRITE(numout,*) '   ==>>>   flux form   : UBS       scheme is used' 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynldf_lap_blp.F90

    r12529 r12822  
    55   !!====================================================================== 
    66   !! History : 3.7  ! 2014-01  (G. Madec, S. Masson)  Original code, re-entrant laplacian 
     7   !!           4.0  ! 2020-04  (A. Nasser, G. Madec)  Add symmetric mixing tensor  
    78   !!---------------------------------------------------------------------- 
    89 
     
    2526   PUBLIC dyn_ldf_lap  ! called by dynldf.F90 
    2627   PUBLIC dyn_ldf_blp  ! called by dynldf.F90 
    27  
     28!!anSYM   
     29   INTEGER, PUBLIC, PARAMETER ::   np_dynldf_lap_rot  = 1         ! div-rot   laplacian  
     30   INTEGER, PUBLIC, PARAMETER ::   np_dynldf_lap_sym  = 2         ! symmetric laplacian (Griffies&Hallberg 2000) 
     31   INTEGER, PUBLIC, PARAMETER ::   np_dynldf_lap_symN = 3         ! symmetric laplacian (cartesian) 
     32    
     33   INTEGER, PUBLIC, PARAMETER ::   ln_dynldf_lap_typ = 1         ! choose type of laplacian (ideally from namelist) 
     34!!anSYM 
    2835   !! * Substitutions 
    2936#  include "do_loop_substitute.h90" 
     
    4451      !! ** Method  :   The Laplacian operator apply on horizontal velocity is  
    4552      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) )  
     53      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) )  
    4654      !! 
    4755      !! ** Action : - pu_rhs, pv_rhs increased by the harmonic operator applied on pu, pv. 
     56      !! 
     57      !! Reference : S.Griffies, R.Hallberg 2000 Mon.Wea.Rev., DOI:/  
    4858      !!---------------------------------------------------------------------- 
    4959      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index 
     
    5767      REAL(wp) ::   zua, zva     ! local scalars 
    5868      REAL(wp), DIMENSION(jpi,jpj) ::   zcur, zdiv 
    59       !!---------------------------------------------------------------------- 
    60       ! 
     69      REAL(wp), DIMENSION(jpi,jpj) ::   zten, zshe   ! tension (diagonal) and shearing (anti-diagonal) terms 
     70      !!---------------------------------------------------------------------- 
     71      ! 
     72!!anSYM TO BE ADDED : reading of laplacian operator (ln_dynldf_lap_typ -> to be written nn_) shall be added in dyn_ldf_init  
     73!!                 as the writing    
     74!!                 and an integer as np_dynldf_lap for instance taken as argument by dyn_ldf_lap call in dyn_ldf 
    6175      IF( kt == nit000 .AND. lwp ) THEN 
    6276         WRITE(numout,*) 
    6377         WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass 
    6478         WRITE(numout,*) '~~~~~~~ ' 
     79         WRITE(numout,*) '                  ln_dynldf_lap_typ = ', ln_dynldf_lap_typ 
     80         SELECT CASE( ln_dynldf_lap_typ )             ! print the choice of operator 
     81         CASE( np_dynldf_lap_rot )   ;   WRITE(numout,*) '   ==>>>   div-rot   laplacian' 
     82         CASE( np_dynldf_lap_sym )   ;   WRITE(numout,*) '   ==>>>   symmetric laplacian (covariant form)' 
     83         CASE( np_dynldf_lap_symN)   ;   WRITE(numout,*) '   ==>>>   symmetric laplacian (simple form)' 
     84         END SELECT 
    6585      ENDIF 
    6686      ! 
     
    6989      ENDIF 
    7090      ! 
    71       !                                                ! =============== 
    72       DO jk = 1, jpkm1                                 ! Horizontal slab 
    73          !                                             ! =============== 
    74          DO_2D_01_01 
    75             !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1) 
    76 !!gm open question here : e3f  at before or now ?    probably now... 
     91      SELECT CASE( ln_dynldf_lap_typ )   
     92         !               
     93         CASE ( np_dynldf_lap_rot )       !==  Vorticity-Divergence form  ==! 
     94            !    
     95            DO jk = 1, jpkm1                                 ! Horizontal slab 
     96               !                                              
     97               DO_2D_01_01 
     98               !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1) 
     99!!gm open question here : e3f  at before or now ?    probably now...  
    77100!!gm note that ahmf has already been multiplied by fmask 
    78             zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       & 
    79                &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  & 
    80                &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  ) 
    81             !                                      ! ahm * div        (computed from 2 to jpi/jpj) 
     101                  zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       & 
     102                     &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  & 
     103                     &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  ) 
     104                  !                                      ! ahm * div        (computed from 2 to jpi/jpj) 
    82105!!gm note that ahmt has already been multiplied by tmask 
    83             zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)                                         & 
    84                &     * (  e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk)  & 
    85                &        + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk)  ) 
    86          END_2D 
     106                  zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)                                         & 
     107                     &     * (  e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk)  & 
     108                     &        + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk)  ) 
     109               END_2D 
     110               ! 
     111               DO_2D_00_00 
     112                  pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * (                                             & 
     113                     &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   & 
     114                     &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                       ) 
     115                     ! 
     116                  pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * (                                             & 
     117                     &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm)   & 
     118                     &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                       ) 
     119               END_2D 
     120               ! 
     121            END DO                                           !   End of slab 
     122            ! 
     123         CASE ( np_dynldf_lap_sym )       !==  Symmetric form  ==!   (Griffies&Hallberg 2000) 
     124            ! 
     125            DO jk = 1, jpkm1                                 ! Horizontal slab 
     126               ! 
     127               DO_2D_01_01 
     128                  !                                      ! shearing stress component (F-point)   NB : ahmf has already been multiplied by fmask 
     129                  zshe(ji-1,jj-1) = ahmf(ji-1,jj-1,jk)                                                              & 
     130                     &     * (    e1f(ji-1,jj-1)    * r1_e2f(ji-1,jj-1)                                             & 
     131                     &         * ( pu(ji-1,jj  ,jk) * r1_e1u(ji-1,jj  )  - pu(ji-1,jj-1,jk) * r1_e1u(ji-1,jj-1) )   & 
     132                     &         +  e2f(ji-1,jj-1)    * r1_e1f(ji-1,jj-1)                                             & 
     133                     &         * ( pv(ji  ,jj-1,jk) * r1_e2v(ji  ,jj-1)  - pv(ji-1,jj-1,jk) * r1_e2v(ji-1,jj-1) )   )  
     134                  !                                      ! tension stress component (T-point)   NB : ahmt has already been multiplied by tmask 
     135                  zten(ji,jj)    = ahmt(ji,jj,jk)                                                       & 
     136                     &     * (    e2t(ji,jj)    * r1_e1t(ji,jj)                                         & 
     137                     &         * ( pu(ji,jj,jk) * r1_e2u(ji,jj)  - pu(ji-1,jj,jk) * r1_e2u(ji-1,jj) )   & 
     138                     &         -  e1t(ji,jj)    * r1_e2t(ji,jj)                                         & 
     139                     &         * ( pv(ji,jj,jk) * r1_e1v(ji,jj)  - pv(ji,jj-1,jk) * r1_e1v(ji,jj-1) )   )    
     140               END_2D 
     141               ! 
     142               DO_2D_00_00 
     143                  pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)                               & 
     144                     &    * (   (   zten(ji+1,jj  ) * e2t(ji+1,jj  )*e2t(ji+1,jj  ) * e3t(ji+1,jj  ,jk,Kmm)                       & 
     145                     &            - zten(ji  ,jj  ) * e2t(ji  ,jj  )*e2t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm) ) * r1_e2u(ji,jj)     &                                                     
     146                     &        + (   zshe(ji  ,jj  ) * e1f(ji  ,jj  )*e1f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                           & 
     147                     &            - zshe(ji  ,jj-1) * e1f(ji  ,jj-1)*e1f(ji  ,jj-1) * e3f(ji  ,jj-1,jk)     ) * r1_e1u(ji,jj) )    
     148                  ! 
     149                  pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                               & 
     150                     &    * (   (   zshe(ji  ,jj  ) * e2f(ji  ,jj  )*e2f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                           & 
     151                     &            - zshe(ji-1,jj  ) * e2f(ji-1,jj  )*e2f(ji-1,jj  ) * e3f(ji-1,jj  ,jk)     ) * r1_e2v(ji,jj)     & 
     152                     &        - (   zten(ji  ,jj+1) * e1t(ji  ,jj+1)*e1t(ji  ,jj+1) * e3t(ji  ,jj+1,jk,Kmm)                       & 
     153                     &            - zten(ji  ,jj  ) * e1t(ji  ,jj  )*e1t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm) ) * r1_e1v(ji,jj) ) 
     154                   ! 
     155               END_2D 
     156               ! 
     157            END DO                                           !   End of slab 
     158            ! 
     159         CASE ( np_dynldf_lap_symN )       !==  Symmetric form  ==!   (naive way) 
     160            ! 
     161            DO jk = 1, jpkm1                                 ! Horizontal slab 
     162               ! 
     163               DO_2D_01_01 
     164                  !                                      ! shearing stress component (F-point)   NB : ahmf has already been multiplied by fmask 
     165                  zshe(ji-1,jj-1) = ahmf(ji-1,jj-1,jk)                                           & 
     166                     &     * (   r1_e2f(ji-1,jj-1) * ( pu(ji-1,jj  ,jk) - pu(ji-1,jj-1,jk)  )   & 
     167                     &         + r1_e1f(ji-1,jj-1) * ( pv(ji  ,jj-1,jk) - pv(ji-1,jj-1,jk)  )   )  
     168                  !                                      ! tension stress component (T-point)   NB : ahmt has already been multiplied by tmask 
     169                  zten(ji,jj)    = ahmt(ji,jj,jk)                                       & 
     170                     &     * (   r1_e1t(ji,jj) * ( pu(ji,jj,jk) - pu(ji-1,jj  ,jk)  )   & 
     171                     &         - r1_e2t(ji,jj) * ( pv(ji,jj,jk) - pv(ji  ,jj-1,jk)  )   )    
     172               END_2D 
     173               ! 
     174               DO_2D_00_00 
     175                  pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   & 
     176                     &    * (   zten(ji+1,jj  ) * e2t(ji+1,jj  ) * e3t(ji+1,jj  ,jk,Kmm)              & 
     177                     &        - zten(ji  ,jj  ) * e2t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm)              & 
     178                     &        + zshe(ji  ,jj  ) * e1f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                  & 
     179                     &        - zshe(ji  ,jj-1) * e1f(ji  ,jj-1) * e3f(ji  ,jj-1,jk)                  )     
     180                  ! 
     181                  pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)   & 
     182                     &    * (   zshe(ji  ,jj  ) * e2f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                  & 
     183                     &        - zshe(ji-1,jj  ) * e2f(ji-1,jj  ) * e3f(ji-1,jj  ,jk)                  & 
     184                     &        - zten(ji  ,jj+1) * e1t(ji  ,jj+1) * e3t(ji  ,jj+1,jk,Kmm)              & 
     185                     &        + zten(ji  ,jj  ) * e1t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm)              ) 
     186                   ! 
     187               END_2D 
     188               ! 
     189            END DO                                           !   End of slab 
     190            !   
     191         CASE DEFAULT                                     ! error 
     192            CALL ctl_stop('STOP','dyn_ldf_lap: wrong value for ln_dynldf_lap_typ'  ) 
     193         END SELECT 
    87194         ! 
    88          DO_2D_00_00 
    89             pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * (                                                 & 
    90                &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   & 
    91                &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                     ) 
    92                ! 
    93             pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * (                                                 & 
    94                &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm)   & 
    95                &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                     ) 
    96          END_2D 
    97          !                                             ! =============== 
    98       END DO                                           !   End of slab 
    99       !                                                ! =============== 
    100195      ! 
    101196   END SUBROUTINE dyn_ldf_lap 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/ldfdyn.F90

    r12529 r12822  
    107107      INTEGER  ::   ioptio, ierr, inum, ios, inn   ! local integer 
    108108      REAL(wp) ::   zah0, zah_max, zUfac           ! local scalar 
     109      REAL(wp) ::   zsum                           ! local scalar 
    109110      CHARACTER(len=5) ::   cl_Units               ! units (m2/s or m4/s) 
    110111      !! 
     
    323324            IF(     ln_dynldf_lap ) THEN                 !   laplacian operator (mask only) 
    324325               ahmt(:,:,1:jpkm1) =       ahmt(:,:,1:jpkm1)   * tmask(:,:,1:jpkm1) 
     326               WRITE(numout,*) ' ahmt tmask ' 
     327!! mask tension at the coast (equivalent of ghostpoints at T) 
     328!              DO jk = 1, jpk 
     329!                 DO jj = 1, jpjm1 
     330!                    DO ji = 1, jpim1      ! NO vector opt. 
     331!                       ! si sum(fmask)==3 = mouillé (on touche pas) 
     332!                       ! si sum = 2 alors on met a 0 zsum = fmask + fmask + fmask,.. et si zsum < 2 -> 0 sinon = 1 
     333!                       zsum =   fmask(ji,jj  ,jk) + fmask(ji+1,jj  ,jk)   & 
     334!                          &   + fmask(ji,jj+1,jk) + fmask(ji+1,jj+1,jk) 
     335!                       IF ( zsum < 2._wp )   ahmt(ji,jj,jk) = 0 
     336!                       ! 
     337!                       !ahmt(ji,jj,jk) = ahmt(ji,jj,jk) * fmask(ji,jj  ,jk) * fmask(ji+1,jj  ,jk)   & 
     338!                       !   &                            * fmask(ji,jj+1,jk) * fmask(ji+1,jj+1,jk) 
     339!                    END DO 
     340!                 END DO 
     341!              END DO 
     342!              ahmt(jpi,:,1:jpkm1) =  0._wp 
     343!              ahmt(:,jpj,1:jpkm1) =  0._wp 
     344!              WRITE(numout,*) '   an45 ahmt x0' 
     345 
    325346               ahmf(:,:,1:jpkm1) =       ahmf(:,:,1:jpkm1)   * fmask(:,:,1:jpkm1) 
     347               WRITE(numout,*) ' ahmf fmask ' 
     348!!an apply no slip at the coast (ssfmask = 1 within the domain and at the coast contrary to fmask in free slip) 
     349!               DO jk = 1, jpkm1 
     350!                  ahmf(:,:,jk) =    ahmf(:,:,jk) * ( 2._wp * ssfmask(:,:) - fmask(:,:,jk) ) 
     351!               END DO 
     352!               WRITE(numout,*) '   an45 ahmf x2' 
     353 
    326354            ELSEIF( ln_dynldf_blp ) THEN                 ! bilaplacian operator (square root + mask) 
    327355               ahmt(:,:,1:jpkm1) = SQRT( ahmt(:,:,1:jpkm1) ) * tmask(:,:,1:jpkm1) 
Note: See TracChangeset for help on using the changeset viewer.