Changeset 13696


Ignore:
Timestamp:
2020-10-28T19:09:27+01:00 (6 months ago)
Author:
techene
Message:

#2555 use same e3f definition as in EEN in ENS and ENE instead of previous ugly fix and change time splitting accordingly change namelist vorticity scheme param nn_een_e3f into nn_e3f_typ

Location:
NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/domzgr_substitute.h90

    r13427 r13696  
    1616#   define  e3v(i,j,k,t)   (e3v_0(i,j,k)*(1._wp+r3v(i,j,t)*vmask(i,j,k))) 
    1717#   define  e3f(i,j,k)     (e3f_0(i,j,k)*(1._wp+r3f(i,j)*fmask(i,j,k))) 
     18#   define  e3f_vor(i,j,k) (e3f_0vor(i,j,k)*(1._wp+r3f(i,j)*fmask(i,j,k))) 
    1819#   define  e3w(i,j,k,t)   (e3w_0(i,j,k)*(1._wp+r3t(i,j,t))) 
    1920#   define  e3uw(i,j,k,t)  (e3uw_0(i,j,k)*(1._wp+r3u(i,j,t))) 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynspg_ts.F90

    r13427 r13696  
    11191119      !! although they should be updated in the variable volume case. Not a big approximation. 
    11201120      !! To remove this approximation, copy lines below inside barotropic loop 
    1121       !! and update depths at T-F points (ht and zhf resp.) at each barotropic time step 
     1121      !! and update depths at T- points (ht) at each barotropic time step 
    11221122      !! 
    11231123      !! Compute zwz = f / ( height of the water colomn ) 
     
    11261126      INTEGER  ::   ji ,jj, jk              ! dummy loop indices 
    11271127      REAL(wp) ::   z1_ht 
    1128       REAL(wp), DIMENSION(jpi,jpj) :: zhf 
    11291128      !!---------------------------------------------------------------------- 
    11301129      ! 
    11311130      SELECT CASE( nvor_scheme ) 
    1132       CASE( np_EEN )                != EEN scheme using e3f energy & enstrophy scheme 
    1133          SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point 
     1131      CASE( np_EEN, np_ENE, np_ENS , np_MIX )   !=  schemes using the same e3f definition 
     1132         SELECT CASE( nn_e3f_typ )                  !* ff_f/e3 at F-point 
    11341133         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4) 
    1135             DO_2D( 1, 0, 1, 0 ) 
     1134            DO_2D( 0, 0, 0, 0 ) 
    11361135               zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1)   & 
    11371136                    &       + ht(ji,jj  ) + ht(ji+1,jj  ) ) * 0.25_wp   
     
    11391138            END_2D 
    11401139         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask) 
    1141             DO_2D( 1, 0, 1, 0 ) 
     1140            DO_2D( 0, 0, 0, 0 ) 
    11421141               zwz(ji,jj) =     (    ht(ji,jj+1) +     ht(ji+1,jj+1)      & 
    11431142                    &            +   ht(ji,jj  ) +     ht(ji+1,jj  )  )   & 
     
    11481147         END SELECT 
    11491148         CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp ) 
     1149      END SELECT 
     1150      ! 
     1151      SELECT CASE( nvor_scheme ) 
     1152      CASE( np_EEN ) 
    11501153         ! 
    11511154         ftne(1,:) = 0._wp   ;   ftnw(1,:) = 0._wp   ;   ftse(1,:) = 0._wp   ;   ftsw(1,:) = 0._wp 
     
    11571160         END_2D 
    11581161         ! 
    1159       CASE( np_EET )                  != EEN scheme using e3t energy conserving scheme 
     1162      CASE( np_EET )                            != EEN scheme using e3t energy conserving scheme 
    11601163         ftne(1,:) = 0._wp   ;   ftnw(1,:) = 0._wp   ;   ftse(1,:) = 0._wp   ;   ftsw(1,:) = 0._wp 
    11611164         DO_2D( 0, 1, 0, 1 ) 
     
    11671170         END_2D 
    11681171         ! 
    1169       CASE( np_ENE, np_ENS , np_MIX )  != all other schemes (ENE, ENS, MIX) except ENT ! 
    1170          ! 
    1171          zwz(:,:) = 0._wp 
    1172          zhf(:,:) = 0._wp 
    1173           
    1174          !!gm  assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed  
    1175 !!gm    A priori a better value should be something like : 
    1176 !!gm          zhf(i,j) = masked sum of  ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)  
    1177 !!gm                     divided by the sum of the corresponding mask  
    1178 !!gm  
    1179 !!             
    1180          IF( .NOT.ln_sco ) THEN 
    1181    
    1182    !!gm  agree the JC comment  : this should be done in a much clear way 
    1183    
    1184    ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 
    1185    !     Set it to zero for the time being  
    1186    !              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level 
    1187    !              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth 
    1188    !              ENDIF 
    1189    !              zhf(:,:) = gdepw_0(:,:,jk+1) 
    1190             ! 
    1191          ELSE 
    1192             ! 
    1193             !zhf(:,:) = hbatf(:,:) 
    1194             DO_2D( 1, 0, 1, 0 ) 
    1195                zhf(ji,jj) =    (   ht_0  (ji,jj  ) + ht_0  (ji+1,jj  )          & 
    1196                     &            + ht_0  (ji,jj+1) + ht_0  (ji+1,jj+1)   )      & 
    1197                     &     / MAX(   ssmask(ji,jj  ) + ssmask(ji+1,jj  )          & 
    1198                     &            + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp  ) 
    1199             END_2D 
    1200          ENDIF 
    1201          ! 
    1202          DO jj = 1, jpjm1   ! keep only the value at the coast if sco 
    1203             zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 
    1204          END DO 
    1205          ! 
    1206          DO jk = 1, jpkm1   ! ocean point : sum of masked e3f 
    1207             DO jj = 1, jpjm1 
    1208                zhf(:,jj) = zhf(:,jj) + e3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 
    1209             END DO 
    1210          END DO 
    1211          CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp ) 
    1212          ! JC: TBC. hf should be greater than 0  
    1213          DO_2D( 1, 1, 1, 1 ) 
    1214             IF( zhf(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zhf(ji,jj) 
    1215          END_2D 
    1216          zwz(:,:) = ff_f(:,:) * zwz(:,:) 
    12171172      END SELECT 
    12181173       
    12191174   END SUBROUTINE dyn_cor_2d_init 
    1220  
    12211175 
    12221176 
     
    14121366   END SUBROUTINE wad_spg 
    14131367      
    1414  
    14151368 
    14161369   SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v ) 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynvor.F90

    r13644 r13696  
    6060   LOGICAL, PUBLIC ::   ln_dynvor_eeT   !: t-point energy conserving scheme     (EET) 
    6161   LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy & enstrophy conserving scheme (EEN) 
    62    INTEGER, PUBLIC ::      nn_een_e3f      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 
    6362   LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme                         (MIX) 
    6463   LOGICAL, PUBLIC ::   ln_dynvor_msk   !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 
     64   INTEGER, PUBLIC ::   nn_e3f_typ      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 
    6565 
    6666   INTEGER, PUBLIC ::   nvor_scheme     !: choice of the type of advection scheme 
     
    8484   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       -  
    8585   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2v_2e1e2f   ! = di(e2u)/(2*e1e2f)  used in F-point metric term calculation 
    86    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1v)/(2*e1e2f)   -        -      -       -  
     86   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1v)/(2*e1e2f)   -        -      -       - 
     87   ! 
     88   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   e3f_0vor   ! e3f used in EEN, ENE and ENS cases (key_qco only) 
    8789    
    8890   REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4 
     
    335337      ! 
    336338      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    337       REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars 
     339      REAL(wp) ::   zx1, zy1, zx2, zy2, ze3f, zmsk   ! local scalars 
    338340      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz   ! 2D workspace 
    339341      !!---------------------------------------------------------------------- 
     
    387389         ! 
    388390#if defined key_qco 
    389          !                                   !==  horizontal fluxes and potential vorticity  ==! 
     391         DO_2D( 1, 0, 1, 0 )                 !==  potential vorticity  ==!   (key_qco) 
     392            zwz(ji,jj) = zwz(ji,jj) / e3f_vor(ji,jj,jk) 
     393         END_2D 
     394#else 
     395         SELECT CASE( nn_e3f_typ  )           !==  potential vorticity  ==! 
     396         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4) 
     397            DO_2D( 1, 0, 1, 0 ) 
     398               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   & 
     399                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   & 
     400                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   & 
     401                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  ) 
     402               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * 4._wp / ze3f 
     403               ELSE                       ;   zwz(ji,jj) = 0._wp 
     404               ENDIF 
     405            END_2D 
     406         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask) 
     407            DO_2D( 1, 0, 1, 0 ) 
     408               ze3f = (   e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   & 
     409                  &     + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   & 
     410                  &     + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   & 
     411                  &     + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)   ) 
     412               zmsk = (   tmask(ji,jj+1,jk)   + tmask(ji+1,jj+1,jk)   & 
     413                  &     + tmask(ji,jj  ,jk)   + tmask(ji+1,jj  ,jk)   ) 
     414               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * zmsk / ze3f 
     415               ELSE                       ;   zwz(ji,jj) = 0._wp 
     416               ENDIF 
     417            END_2D 
     418         END SELECT 
     419#endif 
     420         !                                   !==  horizontal fluxes  ==! 
    390421         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
    391422         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
    392          zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
    393 #else 
    394          !                                   !==  horizontal fluxes and potential vorticity  ==! 
    395          IF( ln_sco ) THEN                        ! weighted by e3 
    396             zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
    397             zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
    398             zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
    399          ELSE                                     ! ugly fix: unweighted by e3 
    400             zwx(:,:) = e2u(:,:) * pu(:,:,jk) 
    401             zwy(:,:) = e1v(:,:) * pv(:,:,jk) 
    402          ENDIF 
    403 #endif 
    404423         ! 
    405424         !                                   !==  compute and add the vorticity term trend  =! 
     
    445464      ! 
    446465      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    447       REAL(wp) ::   zuav, zvau   ! local scalars 
     466      REAL(wp) ::   zuav, zvau, ze3f, zmsk   ! local scalars 
    448467      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz, zww   ! 2D workspace 
    449468      !!---------------------------------------------------------------------- 
     
    497516         ! 
    498517#if defined key_qco 
    499          !                                   !==  horizontal fluxes and potential vorticity  ==! 
     518         DO_2D( 1, 0, 1, 0 )                 !==  potential vorticity  ==!   (key_qco) 
     519            zwz(ji,jj) = zwz(ji,jj) / e3f_vor(ji,jj,jk) 
     520         END_2D 
     521#else 
     522         SELECT CASE( nn_e3f_typ )           !==  potential vorticity  ==! 
     523         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4) 
     524            DO_2D( 1, 0, 1, 0 ) 
     525               ze3f = (  e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   & 
     526                  &    + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   & 
     527                  &    + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   & 
     528                  &    + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  ) 
     529               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * 4._wp / ze3f 
     530               ELSE                       ;   zwz(ji,jj) = 0._wp 
     531               ENDIF 
     532            END_2D 
     533         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask) 
     534            DO_2D( 1, 0, 1, 0 ) 
     535               ze3f = (   e3t(ji  ,jj+1,jk,Kmm)*tmask(ji  ,jj+1,jk)   & 
     536                  &     + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   & 
     537                  &     + e3t(ji  ,jj  ,jk,Kmm)*tmask(ji  ,jj  ,jk)   & 
     538                  &     + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)   ) 
     539               zmsk = (   tmask(ji,jj+1,jk)   + tmask(ji+1,jj+1,jk)   & 
     540                  &     + tmask(ji,jj  ,jk)   + tmask(ji+1,jj  ,jk)   ) 
     541               IF( ze3f /= 0._wp ) THEN   ;   zwz(ji,jj) = zwz(ji,jj) * zmsk / ze3f 
     542               ELSE                       ;   zwz(ji,jj) = 0._wp 
     543               ENDIF 
     544            END_2D 
     545         END SELECT 
     546#endif 
     547         !                                   !==  horizontal fluxes  ==! 
    500548         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
    501549         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
    502          zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
    503 #else 
    504          !                                   !==  horizontal fluxes and potential vorticity  ==! 
    505          IF( ln_sco ) THEN                        ! weighted by e3 
    506             zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
    507             zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
    508             zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
    509          ELSE                                     ! ugly fix: unweighted by e3 
    510             zwx(:,:) = e2u(:,:) * pu(:,:,jk) 
    511             zwy(:,:) = e1v(:,:) * pv(:,:,jk) 
    512          ENDIF 
    513 #endif 
    514550         ! 
    515551         !                                   !==  compute and add the vorticity term trend  =! 
     
    570606         !                                             ! =============== 
    571607         ! 
    572          SELECT CASE( nn_een_e3f )           ! == reciprocal of e3 at F-point 
     608#if defined key_qco 
     609         DO_2D( 1, 0, 1, 0 )                 ! == reciprocal of e3 at F-point (key_qco) 
     610            z1_e3f(ji,jj) = 1._wp / e3f_vor(ji,jj,jk) 
     611         END_2D 
     612#else 
     613         SELECT CASE( nn_e3f_typ )           ! == reciprocal of e3 at F-point 
    573614         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4) 
    574615            DO_2D( 1, 0, 1, 0 ) 
     
    594635            END_2D 
    595636         END SELECT 
     637#endif 
    596638         ! 
    597639         SELECT CASE( kvor )                 !==  vorticity considered  ==! 
     
    801843      INTEGER ::   ji, jj, jk    ! dummy loop indices 
    802844      INTEGER ::   ioptio, ios   ! local integer 
     845      REAL(wp) ::   zmsk    ! local scalars 
    803846      !! 
    804847      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT,   & 
    805          &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk 
     848         &                 ln_dynvor_een, nn_e3f_typ   , ln_dynvor_mix, ln_dynvor_msk 
    806849      !!---------------------------------------------------------------------- 
    807850      ! 
     
    825868         WRITE(numout,*) '      energy conserving scheme  (een using e3t)      ln_dynvor_eeT = ', ln_dynvor_eeT 
    826869         WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een 
    827          WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f 
     870         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_e3f_typ = ', nn_e3f_typ 
    828871         WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix 
    829872         WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk 
     
    891934         ! 
    892935      END SELECT 
    893        
     936!!st ADD !! pense a recalculer le hf_0 
     937#if defined key_qco 
     938      SELECT CASE( nvor_scheme )    ! qco case: pre-computed a specific e3f_0 for some vorticity schemes 
     939      CASE( np_ENS , np_ENE , np_EEN , np_MIX ) 
     940         ! 
     941         ALLOCATE( e3f_0vor(jpi,jpj,jpk) ) 
     942         ! 
     943         SELECT CASE( nn_e3f_typ ) 
     944         CASE ( 0 )                        ! original formulation  (masked averaging of e3t divided by 4) 
     945            DO_3D( 0, 0, 0, 0, 1, jpk ) 
     946               e3f_0vor(ji,jj,jk) = (   e3t_0(ji  ,jj+1,jk)*tmask(ji  ,jj+1,jk)   & 
     947                  &                   + e3t_0(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
     948                  &                   + e3t_0(ji  ,jj  ,jk)*tmask(ji  ,jj  ,jk)   & 
     949                  &                   + e3t_0(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)   ) * 0.25_wp 
     950            END_3D 
     951         CASE ( 1 )                        ! new formulation  (masked averaging of e3t divided by the sum of mask) 
     952            DO_3D( 0, 0, 0, 0, 1, jpk ) 
     953               zmsk = (tmask(ji,jj+1,jk) +tmask(ji+1,jj+1,jk)   & 
     954                  &  + tmask(ji,jj  ,jk) +tmask(ji+1,jj  ,jk)  ) 
     955               ! 
     956               IF( zmsk /= 0._wp ) THEN  
     957                  e3f_0vor(ji,jj,jk) = (   e3t_0(ji  ,jj+1,jk)*tmask(ji  ,jj+1,jk)   & 
     958                     &                   + e3t_0(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
     959                     &                   + e3t_0(ji  ,jj  ,jk)*tmask(ji  ,jj  ,jk)   & 
     960                     &                   + e3t_0(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)   ) / zmsk 
     961               ENDIF 
     962            END_3D 
     963         END SELECT 
     964         ! 
     965         CALL lbc_lnk( 'dynvor', e3f_0vor, 'F', 1._wp ) 
     966         !                                 ! insure e3f_0vor /= 0 
     967         WHERE( e3f_0vor(:,:,:) == 0._wp )   e3f_0vor(:,:,:) = e3f_0(:,:,:) 
     968         ! 
     969      END SELECT 
     970      ! 
     971#endif 
     972!!st END       
    894973      IF(lwp) THEN                   ! Print the choice 
    895974         WRITE(numout,*) 
Note: See TracChangeset for help on using the changeset viewer.