Changeset 13197


Ignore:
Timestamp:
2020-07-01T16:09:00+02:00 (4 months ago)
Author:
gsamson
Message:

merge with trunk@r13136 with a more recent svn version; pass all SETTE tests; results identical to trunk@r13136; ticket #2419

Location:
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement
Files:
1 deleted
47 edited
1 copied

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/cfgs/C1D_PAPA/MY_SRC/usrdef_zgr.F90

    r12377 r13197  
    3030   PUBLIC   usr_def_zgr        ! called by domzgr.F90 
    3131 
     32   !! * Substitutions 
     33#  include "do_loop_substitute.h90" 
    3234   !!---------------------------------------------------------------------- 
    3335   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    157159         pe3vw(:,:,jk) = pe3w_1d (jk) 
    158160      END DO 
    159       DO jj = 1, jpj                      ! bottom scale factors and depth at T- and W-points 
    160          DO ji = 1, jpi 
    161             ik = k_bot(ji,jj) 
    162             pdepw(ji,jj,ik+1) = MIN( zht(ji,jj) , pdepw_1d(ik+1) ) 
    163             pe3t (ji,jj,ik  ) = pdepw(ji,jj,ik+1) - pdepw(ji,jj,ik) 
    164             pe3t (ji,jj,ik+1) = pe3t (ji,jj,ik  )  
    165             ! 
    166             pdept(ji,jj,ik  ) = pdepw(ji,jj,ik  ) + pe3t (ji,jj,ik  ) * 0.5_wp 
    167             pdept(ji,jj,ik+1) = pdepw(ji,jj,ik+1) + pe3t (ji,jj,ik+1) * 0.5_wp 
    168             pe3w (ji,jj,ik+1) = pdept(ji,jj,ik+1) - pdept(ji,jj,ik)              ! = pe3t (ji,jj,ik  ) 
    169          END DO 
    170       END DO          
     161      ! bottom scale factors and depth at T- and W-points 
     162      DO_2D_11_11 
     163         ik = k_bot(ji,jj) 
     164         pdepw(ji,jj,ik+1) = MIN( zht(ji,jj) , pdepw_1d(ik+1) ) 
     165         pe3t (ji,jj,ik  ) = pdepw(ji,jj,ik+1) - pdepw(ji,jj,ik) 
     166         pe3t (ji,jj,ik+1) = pe3t (ji,jj,ik  )  
     167         ! 
     168         pdept(ji,jj,ik  ) = pdepw(ji,jj,ik  ) + pe3t (ji,jj,ik  ) * 0.5_wp 
     169         pdept(ji,jj,ik+1) = pdepw(ji,jj,ik+1) + pe3t (ji,jj,ik+1) * 0.5_wp 
     170         pe3w (ji,jj,ik+1) = pdept(ji,jj,ik+1) - pdept(ji,jj,ik)              ! = pe3t (ji,jj,ik  ) 
     171      END_2D         
    171172      !                                   ! bottom scale factors and depth at  U-, V-, UW and VW-points 
    172173      !                                   ! usually Computed as the minimum of neighbooring scale factors 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk_algo_coare3p0.F90

    r13159 r13197  
    395395      ! 
    396396      DO_2D_11_11 
    397          ! 
    398          zw = pwnd(ji,jj)   ! wind speed 
    399          ! 
    400          ! Charnock's constant, increases with the wind : 
    401          zgt10 = 0.5 + SIGN(0.5_wp,(zw - 10))  ! If zw<10. --> 0, else --> 1 
    402          zgt18 = 0.5 + SIGN(0.5_wp,(zw - 18.)) ! If zw<18. --> 0, else --> 1 
    403          ! 
    404          alfa_charn_3p0(ji,jj) =  (1. - zgt10)*0.011    &    ! wind is lower than 10 m/s 
    405             &     + zgt10*((1. - zgt18)*(0.011 + (0.018 - 0.011) & 
    406             &      *(zw - 10.)/(18. - 10.)) + zgt18*( 0.018 ) )    ! Hare et al. (1999) 
    407          ! 
     397      ! 
     398      zw = pwnd(ji,jj)   ! wind speed 
     399      ! 
     400      ! Charnock's constant, increases with the wind : 
     401      zgt10 = 0.5 + SIGN(0.5_wp,(zw - 10))  ! If zw<10. --> 0, else --> 1 
     402      zgt18 = 0.5 + SIGN(0.5_wp,(zw - 18.)) ! If zw<18. --> 0, else --> 1 
     403      ! 
     404      alfa_charn_3p0(ji,jj) =  (1. - zgt10)*0.011    &    ! wind is lower than 10 m/s 
     405         &     + zgt10*((1. - zgt18)*(0.011 + (0.018 - 0.011) & 
     406         &      *(zw - 10.)/(18. - 10.)) + zgt18*( 0.018 ) )    ! Hare et al. (1999) 
     407      ! 
    408408      END_2D 
    409409      ! 
     
    431431      ! 
    432432      DO_2D_11_11 
    433          ! 
    434          zta = pzeta(ji,jj) 
    435          ! 
    436          zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable 
    437          ! 
    438          zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   & 
    439             & - 2.*ATAN(zphi_m) + 0.5*rpi 
    440          ! 
    441          zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective 
    442          ! 
    443          zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 
    444             &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 
    445          ! 
    446          zf = zta*zta 
    447          zf = zf/(1. + zf) 
    448          zc = MIN(50._wp, 0.35_wp*zta) 
    449          zstab = 0.5 + SIGN(0.5_wp, zta) 
    450          ! 
    451          psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0) 
    452             &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0) 
    453             &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     " 
    454          ! 
     433      ! 
     434      zta = pzeta(ji,jj) 
     435      ! 
     436      zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable 
     437      ! 
     438      zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   & 
     439         & - 2.*ATAN(zphi_m) + 0.5*rpi 
     440      ! 
     441      zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective 
     442      ! 
     443      zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 
     444         &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 
     445      ! 
     446      zf = zta*zta 
     447      zf = zf/(1. + zf) 
     448      zc = MIN(50._wp, 0.35_wp*zta) 
     449      zstab = 0.5 + SIGN(0.5_wp, zta) 
     450      ! 
     451      psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0) 
     452         &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0) 
     453         &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     " 
     454      ! 
    455455      END_2D 
    456456      ! 
     
    482482      ! 
    483483      DO_2D_11_11 
    484          ! 
    485          zta = pzeta(ji,jj) 
    486          ! 
    487          zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable) 
    488          ! 
    489          zpsi_k = 2.*LOG((1. + zphi_h)/2.) 
    490          ! 
    491          zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective 
    492          ! 
    493          zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 
    494             &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 
    495          ! 
    496          zf = zta*zta 
    497          zf = zf/(1. + zf) 
    498          zc = MIN(50._wp,0.35_wp*zta) 
    499          zstab = 0.5 + SIGN(0.5_wp, zta) 
    500          ! 
    501          psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & 
    502             &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     & 
    503             &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 ) 
    504          ! 
     484      ! 
     485      zta = pzeta(ji,jj) 
     486      ! 
     487      zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable) 
     488      ! 
     489      zpsi_k = 2.*LOG((1. + zphi_h)/2.) 
     490      ! 
     491      zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective 
     492      ! 
     493      zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 
     494         &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 
     495      ! 
     496      zf = zta*zta 
     497      zf = zf/(1. + zf) 
     498      zc = MIN(50._wp,0.35_wp*zta) 
     499      zstab = 0.5 + SIGN(0.5_wp, zta) 
     500      ! 
     501      psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & 
     502         &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     & 
     503         &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 ) 
     504      ! 
    505505      END_2D 
    506506      ! 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk_algo_coare3p6.F90

    r13159 r13197  
    431431      ! 
    432432      DO_2D_11_11 
    433          ! 
    434          zta = pzeta(ji,jj) 
    435          ! 
    436          zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable 
    437          ! 
    438          zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   & 
    439             & - 2.*ATAN(zphi_m) + 0.5*rpi 
    440          ! 
    441          zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective 
    442          ! 
    443          zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 
    444             &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 
    445          ! 
    446          zf = zta*zta 
    447          zf = zf/(1. + zf) 
    448          zc = MIN(50._wp, 0.35_wp*zta) 
    449          zstab = 0.5 + SIGN(0.5_wp, zta) 
    450          ! 
    451          psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0) 
    452             &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0) 
    453             &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     " 
    454          ! 
     433      ! 
     434      zta = pzeta(ji,jj) 
     435      ! 
     436      zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable 
     437      ! 
     438      zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   & 
     439         & - 2.*ATAN(zphi_m) + 0.5*rpi 
     440      ! 
     441      zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective 
     442      ! 
     443      zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 
     444         &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 
     445      ! 
     446      zf = zta*zta 
     447      zf = zf/(1. + zf) 
     448      zc = MIN(50._wp, 0.35_wp*zta) 
     449      zstab = 0.5 + SIGN(0.5_wp, zta) 
     450      ! 
     451      psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0) 
     452         &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0) 
     453         &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     " 
     454      ! 
    455455      END_2D 
    456456      ! 
     
    482482      ! 
    483483      DO_2D_11_11 
    484          ! 
    485          zta = pzeta(ji,jj) 
    486          ! 
    487          zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable) 
    488          ! 
    489          zpsi_k = 2.*LOG((1. + zphi_h)/2.) 
    490          ! 
    491          zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective 
    492          ! 
    493          zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 
    494             &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 
    495          ! 
    496          zf = zta*zta 
    497          zf = zf/(1. + zf) 
    498          zc = MIN(50._wp,0.35_wp*zta) 
    499          zstab = 0.5 + SIGN(0.5_wp, zta) 
    500          ! 
    501          psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & 
    502             &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     & 
    503             &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 ) 
    504          ! 
     484      ! 
     485      zta = pzeta(ji,jj) 
     486      ! 
     487      zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable) 
     488      ! 
     489      zpsi_k = 2.*LOG((1. + zphi_h)/2.) 
     490      ! 
     491      zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective 
     492      ! 
     493      zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 
     494         &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 
     495      ! 
     496      zf = zta*zta 
     497      zf = zf/(1. + zf) 
     498      zc = MIN(50._wp,0.35_wp*zta) 
     499      zstab = 0.5 + SIGN(0.5_wp, zta) 
     500      ! 
     501      psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & 
     502         &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     & 
     503         &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 ) 
     504      ! 
    505505      END_2D 
    506506      ! 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk_algo_ecmwf.F90

    r13159 r13197  
    411411      !!---------------------------------------------------------------------------------- 
    412412      DO_2D_11_11 
    413          ! 
    414          zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!): 
    415          ! 
    416          ! Unstable (Paulson 1970): 
    417          !   eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 
    418          zx = SQRT(ABS(1._wp - 16._wp*zzeta)) 
    419          ztmp = 1._wp + SQRT(zx) 
    420          ztmp = ztmp*ztmp 
    421          psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) )   & 
    422             &       -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi 
    423          ! 
    424          ! Unstable: 
    425          ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 
    426          psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & 
    427             &       - zzeta - 2._wp/3._wp*5._wp/0.35_wp 
    428          ! 
    429          ! Combining: 
    430          stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 
    431          ! 
    432          psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 
    433             &                +      stab  * psi_stab      ! (zzeta > 0) Stable 
    434          ! 
     413      ! 
     414      zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!): 
     415      ! 
     416      ! Unstable (Paulson 1970): 
     417      !   eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 
     418      zx = SQRT(ABS(1._wp - 16._wp*zzeta)) 
     419      ztmp = 1._wp + SQRT(zx) 
     420      ztmp = ztmp*ztmp 
     421      psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) )   & 
     422         &       -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi 
     423      ! 
     424      ! Unstable: 
     425      ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 
     426      psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & 
     427         &       - zzeta - 2._wp/3._wp*5._wp/0.35_wp 
     428      ! 
     429      ! Combining: 
     430      stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 
     431      ! 
     432      psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 
     433         &                +      stab  * psi_stab      ! (zzeta > 0) Stable 
     434      ! 
    435435      END_2D 
    436436   END FUNCTION psi_m_ecmwf 
     
    456456      ! 
    457457      DO_2D_11_11 
    458          ! 
    459          zzeta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!): 
    460          ! 
    461          zx  = ABS(1._wp - 16._wp*zzeta)**.25        ! this is actually (1/phi_m)**2  !!! 
    462          !                                     ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1 
    463          ! Unstable (Paulson 1970) : 
    464          psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx))   ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 
    465          ! 
    466          ! Stable: 
    467          psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 
    468             &       - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp 
    469          ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution... 
    470          ! 
    471          stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 
    472          ! 
    473          ! 
    474          psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst &   ! (zzeta < 0) Unstable 
    475             &                +    stab    * psi_stab        ! (zzeta > 0) Stable 
    476          ! 
     458      ! 
     459      zzeta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!): 
     460      ! 
     461      zx  = ABS(1._wp - 16._wp*zzeta)**.25        ! this is actually (1/phi_m)**2  !!! 
     462      !                                     ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1 
     463      ! Unstable (Paulson 1970) : 
     464      psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx))   ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 
     465      ! 
     466      ! Stable: 
     467      psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 
     468         &       - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp 
     469      ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution... 
     470      ! 
     471      stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 
     472      ! 
     473      ! 
     474      psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst &   ! (zzeta < 0) Unstable 
     475         &                +    stab    * psi_stab        ! (zzeta > 0) Stable 
     476      ! 
    477477      END_2D 
    478478   END FUNCTION psi_h_ecmwf 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk_phy.F90

    r13159 r13197  
    524524            &                    pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj),                & 
    525525            &                    pTau(ji,jj), zQsen, zQlat ) 
    526  
     526          
    527527         zTs2  = pTs(ji,jj)*pTs(ji,jj) 
    528528         zQlw  = emiss_w*(prlw(ji,jj) - stefan*zTs2*zTs2) ! Net longwave flux 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/BENCH/MY_SRC/usrdef_hgr.F90

    r9762 r13197  
    2424   PUBLIC   usr_def_hgr   ! called by domhgr.F90 
    2525 
     26   !! * Substitutions 
     27#  include "do_loop_substitute.h90" 
    2628   !!---------------------------------------------------------------------- 
    2729   !! NEMO/OPA 4.0, NEMO Consortium (2016) 
     
    7274      !                           
    7375      ! Position coordinates (in grid points) 
    74       !                          ==========          
    75       DO jj = 1, jpj 
    76          DO ji = 1, jpi 
    77              
    78             zti = REAL( ji - 1 + nimpp - 1, wp )          ;  ztj = REAL( jj - 1 + njmpp - 1, wp ) 
    79             zui = REAL( ji - 1 + nimpp - 1, wp ) + 0.5_wp ;  zvj = REAL( jj - 1 + njmpp - 1, wp ) + 0.5_wp 
     76      !                          ========== 
     77      DO_2D_11_11 
     78          
     79         zti = REAL( ji - 1 + nimpp - 1, wp )          ;  ztj = REAL( jj - 1 + njmpp - 1, wp ) 
     80         zui = REAL( ji - 1 + nimpp - 1, wp ) + 0.5_wp ;  zvj = REAL( jj - 1 + njmpp - 1, wp ) + 0.5_wp 
     81          
     82         plamt(ji,jj) = zti 
     83         plamu(ji,jj) = zui 
     84         plamv(ji,jj) = zti 
     85         plamf(ji,jj) = zui 
     86          
     87         pphit(ji,jj) = ztj 
     88         pphiv(ji,jj) = zvj 
     89         pphiu(ji,jj) = ztj 
     90         pphif(ji,jj) = zvj 
    8091 
    81             plamt(ji,jj) = zti 
    82             plamu(ji,jj) = zui 
    83             plamv(ji,jj) = zti 
    84             plamf(ji,jj) = zui 
    85     
    86             pphit(ji,jj) = ztj 
    87             pphiv(ji,jj) = zvj 
    88             pphiu(ji,jj) = ztj 
    89             pphif(ji,jj) = zvj 
    90              
    91          END DO 
    92       END DO 
     92      END_2D 
    9393      !      
    9494      ! Horizontal scale factors (in meters) 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/BENCH/MY_SRC/usrdef_istate.F90

    r11536 r13197  
    2828   PUBLIC   usr_def_istate   ! called by istate.F90 
    2929 
     30   !! * Substitutions 
     31#  include "do_loop_substitute.h90" 
    3032   !!---------------------------------------------------------------------- 
    3133   !! NEMO/OPA 4.0 , NEMO Consortium (2016) 
     
    6264      ! 
    6365      ! define unique value on each point. z2d ranging from 0.05 to -0.05 
    64       DO jj = 1, jpj 
    65          DO ji = 1, jpi 
    66             z2d(ji,jj) = 0.1 * ( 0.5 - REAL( mig(ji) + mjg(jj) * jpiglo, wp ) / REAL( jpiglo * jpjglo, wp ) ) 
    67          ENDDO 
    68       ENDDO 
     66      DO_2D_11_11 
     67         z2d(ji,jj) = 0.1 * ( 0.5 - REAL( mig(ji) + (mjg(jj)-1) * jpiglo, wp ) / REAL( jpiglo * jpjglo, wp ) ) 
     68      END_2D 
    6969      ! 
    7070      ! sea level: 
     
    7878         pts(:,:,jk,jp_sal) = 30._wp + 1._wp * zfact + z2d(:,:)           ! 30 to 31 +/- 0.05 psu 
    7979         ! velocities: 
    80          pu(:,:,jk) = z2d(:,:) * 0.1_wp                                   ! +/- 0.005  m/s 
    81          pv(:,:,jk) = z2d(:,:) * 0.01_wp                                  ! +/- 0.0005 m/s 
     80         pu(:,:,jk) = z2d(:,:) *  0.1_wp * umask(:,:,jk)                  ! +/- 0.005  m/s 
     81         pv(:,:,jk) = z2d(:,:) * 0.01_wp * vmask(:,:,jk)                  ! +/- 0.0005 m/s 
    8282      ENDDO 
    8383      ! 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/BENCH/MY_SRC/usrdef_sbc.F90

    r12377 r13197  
    3434   PUBLIC   usrdef_sbc_ice_flx  ! routine called by sbcice_lim.F90 for ice thermo 
    3535 
     36   !! * Substitutions 
     37#  include "do_loop_substitute.h90" 
    3638   !!---------------------------------------------------------------------- 
    3739   !! NEMO/OPA 4.0 , NEMO Consortium (2016) 
     
    102104      ! 
    103105      ! define unique value on each point. z2d ranging from 0.05 to -0.05 
    104       DO jj = 1, jpj 
    105          DO ji = 1, jpi 
    106             z2d(ji,jj) = 0.1 * ( 0.5 - REAL( nimpp + ji - 1 + ( njmpp + jj - 2 ) * jpiglo, wp ) / REAL( jpiglo * jpjglo, wp ) ) 
    107          ENDDO 
    108       ENDDO 
     106      DO_2D_11_11 
     107         z2d(ji,jj) = 0.1 * ( 0.5 - REAL( nimpp + ji - 1 + ( njmpp + jj - 2 ) * jpiglo, wp ) / REAL( jpiglo * jpjglo, wp ) ) 
     108      END_2D 
    109109      utau_ice(:,:) = 0.1_wp +  z2d(:,:) 
    110110      vtau_ice(:,:) = 0.1_wp +  z2d(:,:) 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/CANAL/MY_SRC/diawri.F90

    r12489 r13197  
    2626   !!---------------------------------------------------------------------- 
    2727   USE oce            ! ocean dynamics and tracers  
     28   USE isf_oce 
     29   USE isfcpl 
     30   USE abl            ! abl variables in case ln_abl = .true. 
    2831   USE dom_oce        ! ocean space and time domain 
    2932   USE phycst         ! physical constants 
     
    6568   PUBLIC   dia_wri_state 
    6669   PUBLIC   dia_wri_alloc           ! Called by nemogcm module 
    67  
     70#if ! defined key_iomput    
     71   PUBLIC   dia_wri_alloc_abl       ! Called by sbcabl  module (if ln_abl = .true.) 
     72#endif 
    6873   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file 
    6974   INTEGER ::          nb_T              , ndim_bT   ! grid_T file 
     
    7176   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file 
    7277   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file 
     78   INTEGER ::   nid_A, nz_A, nh_A, ndim_A, ndim_hA   ! grid_ABL file    
    7379   INTEGER ::   ndex(1)                              ! ??? 
    7480   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV 
     81   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hA, ndex_A ! ABL 
    7582   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V 
    7683   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT 
    7784 
     85   !! * Substitutions 
     86#  include "do_loop_substitute.h90" 
    7887   !!---------------------------------------------------------------------- 
    7988   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    147156      CALL iom_put(  "sst", ts(:,:,1,jp_tem,Kmm) )    ! surface temperature 
    148157      IF ( iom_use("sbt") ) THEN 
    149          DO jj = 1, jpj 
    150             DO ji = 1, jpi 
    151                ikbot = mbkt(ji,jj) 
    152                z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm) 
    153             END DO 
    154          END DO 
     158         DO_2D_11_11 
     159            ikbot = mbkt(ji,jj) 
     160            z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm) 
     161         END_2D 
    155162         CALL iom_put( "sbt", z2d )                ! bottom temperature 
    156163      ENDIF 
     
    159166      CALL iom_put(  "sss", ts(:,:,1,jp_sal,Kmm) )    ! surface salinity 
    160167      IF ( iom_use("sbs") ) THEN 
    161          DO jj = 1, jpj 
    162             DO ji = 1, jpi 
    163                ikbot = mbkt(ji,jj) 
    164                z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm) 
    165             END DO 
    166          END DO 
     168         DO_2D_11_11 
     169            ikbot = mbkt(ji,jj) 
     170            z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm) 
     171         END_2D 
    167172         CALL iom_put( "sbs", z2d )                ! bottom salinity 
    168173      ENDIF 
     
    171176         zztmp = rho0 * 0.25 
    172177         z2d(:,:) = 0._wp 
    173          DO jj = 2, jpjm1 
    174             DO ji = fs_2, fs_jpim1   ! vector opt. 
    175                zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * uu(ji  ,jj,mbku(ji  ,jj),Kmm)  )**2   & 
    176                   &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm)  )**2   & 
    177                   &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vv(ji,jj  ,mbkv(ji,jj  ),Kmm)  )**2   & 
    178                   &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vv(ji,jj-1,mbkv(ji,jj-1),Kmm)  )**2 
    179                z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)  
    180                ! 
    181             END DO 
    182          END DO 
     178         DO_2D_00_00 
     179            zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * uu(ji  ,jj,mbku(ji  ,jj),Kmm)  )**2   & 
     180               &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm)  )**2   & 
     181               &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vv(ji,jj  ,mbkv(ji,jj  ),Kmm)  )**2   & 
     182               &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vv(ji,jj-1,mbkv(ji,jj-1),Kmm)  )**2 
     183            z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)  
     184            ! 
     185         END_2D 
    183186         CALL lbc_lnk( 'diawri', z2d, 'T', 1. ) 
    184187         CALL iom_put( "taubot", z2d )            
     
    188191      CALL iom_put(  "ssu", uu(:,:,1,Kmm) )            ! surface i-current 
    189192      IF ( iom_use("sbu") ) THEN 
    190          DO jj = 1, jpj 
    191             DO ji = 1, jpi 
    192                ikbot = mbku(ji,jj) 
    193                z2d(ji,jj) = uu(ji,jj,ikbot,Kmm) 
    194             END DO 
    195          END DO 
     193         DO_2D_11_11 
     194            ikbot = mbku(ji,jj) 
     195            z2d(ji,jj) = uu(ji,jj,ikbot,Kmm) 
     196         END_2D 
    196197         CALL iom_put( "sbu", z2d )                ! bottom i-current 
    197198      ENDIF 
     
    200201      CALL iom_put(  "ssv", vv(:,:,1,Kmm) )            ! surface j-current 
    201202      IF ( iom_use("sbv") ) THEN 
    202          DO jj = 1, jpj 
    203             DO ji = 1, jpi 
    204                ikbot = mbkv(ji,jj) 
    205                z2d(ji,jj) = vv(ji,jj,ikbot,Kmm) 
    206             END DO 
    207          END DO 
     203         DO_2D_11_11 
     204            ikbot = mbkv(ji,jj) 
     205            z2d(ji,jj) = vv(ji,jj,ikbot,Kmm) 
     206         END_2D 
    208207         CALL iom_put( "sbv", z2d )                ! bottom j-current 
    209208      ENDIF 
    210209 
     210      IF( ln_zad_Aimp ) ww = ww + wi               ! Recombine explicit and implicit parts of vertical velocity for diagnostic output 
     211      ! 
    211212      CALL iom_put( "woce", ww )                   ! vertical velocity 
    212213      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value 
     
    219220         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 
    220221      ENDIF 
     222      ! 
     223      IF( ln_zad_Aimp ) ww = ww - wi               ! Remove implicit part of vertical velocity that was added for diagnostic output 
    221224 
    222225      CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef. 
     
    227230      IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) ) 
    228231 
    229       IF ( iom_use("salgrad") .OR. iom_use("salgrad2") ) THEN 
    230          z3d(:,:,jpk) = 0. 
    231          DO jk = 1, jpkm1 
    232             DO jj = 2, jpjm1                                    ! sal gradient 
    233                DO ji = fs_2, fs_jpim1   ! vector opt. 
    234                   zztmp  = ts(ji,jj,jk,jp_sal,Kmm) 
    235                   zztmpx = ( ts(ji+1,jj,jk,jp_sal,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj  ,jk,jp_sal,Kmm) ) * r1_e1u(ji-1,jj) 
    236                   zztmpy = ( ts(ji,jj+1,jk,jp_sal,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji  ,jj-1,jk,jp_sal,Kmm) ) * r1_e2v(ji,jj-1) 
    237                   z3d(ji,jj,jk) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   & 
    238                      &                 * umask(ji,jj,jk) * umask(ji-1,jj,jk) * vmask(ji,jj,jk) * umask(ji,jj-1,jk) 
    239                END DO 
    240             END DO 
    241          END DO 
    242          CALL lbc_lnk( 'diawri', z3d, 'T', 1. ) 
    243          CALL iom_put( "salgrad2",  z3d )          ! square of module of sal gradient 
    244          z3d(:,:,:) = SQRT( z3d(:,:,:) ) 
    245          CALL iom_put( "salgrad" ,  z3d )          ! module of sal gradient 
    246       ENDIF 
    247           
    248232      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN 
    249          DO jj = 2, jpjm1                                    ! sst gradient 
    250             DO ji = fs_2, fs_jpim1   ! vector opt. 
    251                zztmp  = ts(ji,jj,1,jp_tem,Kmm) 
    252                zztmpx = ( ts(ji+1,jj,1,jp_tem,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj  ,1,jp_tem,Kmm) ) * r1_e1u(ji-1,jj) 
    253                zztmpy = ( ts(ji,jj+1,1,jp_tem,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji  ,jj-1,1,jp_tem,Kmm) ) * r1_e2v(ji,jj-1) 
    254                z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   & 
    255                   &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1) 
    256             END DO 
    257          END DO 
     233         DO_2D_00_00 
     234            zztmp  = ts(ji,jj,1,jp_tem,Kmm) 
     235            zztmpx = ( ts(ji+1,jj,1,jp_tem,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj  ,1,jp_tem,Kmm) ) * r1_e1u(ji-1,jj) 
     236            zztmpy = ( ts(ji,jj+1,1,jp_tem,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji  ,jj-1,1,jp_tem,Kmm) ) * r1_e2v(ji,jj-1) 
     237            z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   & 
     238               &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1) 
     239         END_2D 
    258240         CALL lbc_lnk( 'diawri', z2d, 'T', 1. ) 
    259241         CALL iom_put( "sstgrad2",  z2d )          ! square of module of sst gradient 
     
    265247      IF( iom_use("heatc") ) THEN 
    266248         z2d(:,:)  = 0._wp  
    267          DO jk = 1, jpkm1 
    268             DO jj = 1, jpj 
    269                DO ji = 1, jpi 
    270                   z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk) 
    271                END DO 
    272             END DO 
    273          END DO 
     249         DO_3D_11_11( 1, jpkm1 ) 
     250            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk) 
     251         END_3D 
    274252         CALL iom_put( "heatc", rho0_rcp * z2d )   ! vertically integrated heat content (J/m2) 
    275253      ENDIF 
     
    277255      IF( iom_use("saltc") ) THEN 
    278256         z2d(:,:)  = 0._wp  
    279          DO jk = 1, jpkm1 
    280             DO jj = 1, jpj 
    281                DO ji = 1, jpi 
    282                   z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 
    283                END DO 
    284             END DO 
    285          END DO 
     257         DO_3D_11_11( 1, jpkm1 ) 
     258            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 
     259         END_3D 
    286260         CALL iom_put( "saltc", rho0 * z2d )          ! vertically integrated salt content (PSU*kg/m2) 
    287261      ENDIF 
     
    289263      IF( iom_use("salt2c") ) THEN 
    290264         z2d(:,:)  = 0._wp  
    291          DO jk = 1, jpkm1 
    292             DO jj = 1, jpj 
    293                DO ji = 1, jpi 
    294                   z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 
    295                END DO 
    296             END DO 
    297          END DO 
     265         DO_3D_11_11( 1, jpkm1 ) 
     266            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 
     267         END_3D 
    298268         CALL iom_put( "salt2c", rho0 * z2d )          ! vertically integrated salt content (PSU*kg/m2) 
    299269      ENDIF 
     
    301271      IF ( iom_use("eken") ) THEN 
    302272         z3d(:,:,jpk) = 0._wp  
    303          DO jk = 1, jpkm1 
    304             DO jj = 2, jpjm1 
    305                DO ji = fs_2, fs_jpim1   ! vector opt. 
    306                   zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    307                   z3d(ji,jj,jk) = zztmp * (  uu(ji-1,jj,jk,Kmm)**2 * e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)   & 
    308                      &                     + uu(ji  ,jj,jk,Kmm)**2 * e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)   & 
    309                      &                     + vv(ji,jj-1,jk,Kmm)**2 * e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)   & 
    310                      &                     + vv(ji,jj  ,jk,Kmm)**2 * e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)   ) 
    311                END DO 
    312             END DO 
    313          END DO 
     273         DO_3D_00_00( 1, jpkm1 ) 
     274            zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     275            z3d(ji,jj,jk) = zztmp * (  uu(ji-1,jj,jk,Kmm)**2 * e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)   & 
     276               &                     + uu(ji  ,jj,jk,Kmm)**2 * e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)   & 
     277               &                     + vv(ji,jj-1,jk,Kmm)**2 * e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)   & 
     278               &                     + vv(ji,jj  ,jk,Kmm)**2 * e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)   ) 
     279         END_3D 
    314280         CALL lbc_lnk( 'diawri', z3d, 'T', 1. ) 
    315281         CALL iom_put( "eken", z3d )                 ! kinetic energy 
     
    321287         z3d(1,:, : ) = 0._wp 
    322288         z3d(:,1, : ) = 0._wp 
    323          DO jk = 1, jpkm1 
    324             DO jj = 2, jpj 
    325                DO ji = 2, jpi 
    326                   z3d(ji,jj,jk) = 0.25_wp * ( uu(ji  ,jj,jk,Kmm) * uu(ji  ,jj,jk,Kmm) * e1e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)  & 
    327                      &                      + uu(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)  & 
    328                      &                      + vv(ji,jj  ,jk,Kmm) * vv(ji,jj  ,jk,Kmm) * e1e2v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)  & 
    329                      &                      + vv(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)  )  & 
    330                      &                    * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) 
    331                END DO 
    332             END DO 
    333          END DO 
    334           
     289         DO_3D_00_00( 1, jpkm1 ) 
     290            z3d(ji,jj,jk) = 0.25_wp * ( uu(ji  ,jj,jk,Kmm) * uu(ji  ,jj,jk,Kmm) * e1e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)  & 
     291               &                      + uu(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)  & 
     292               &                      + vv(ji,jj  ,jk,Kmm) * vv(ji,jj  ,jk,Kmm) * e1e2v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)  & 
     293               &                      + vv(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)  )  & 
     294               &                    * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) 
     295         END_3D 
    335296         CALL lbc_lnk( 'diawri', z3d, 'T', 1. ) 
    336297         CALL iom_put( "ke", z3d ) ! kinetic energy 
    337298 
    338299         z2d(:,:)  = 0._wp  
    339          DO jk = 1, jpkm1 
    340             DO jj = 1, jpj 
    341                DO ji = 1, jpi 
    342                   z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * tmask(ji,jj,jk) 
    343                END DO 
    344             END DO 
    345          END DO 
     300         DO_3D_11_11( 1, jpkm1 ) 
     301            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * tmask(ji,jj,jk) 
     302         END_3D 
    346303         CALL iom_put( "ke_zint", z2d )   ! vertically integrated kinetic energy 
    347304 
     
    353310          
    354311         z3d(:,:,jpk) = 0._wp  
    355          DO jk = 1, jpkm1 
    356             DO jj = 1, jpjm1 
    357                DO ji = 1, fs_jpim1   ! vector opt. 
    358                   z3d(ji,jj,jk) = (   e2v(ji+1,jj  ) * vv(ji+1,jj  ,jk,Kmm) - e2v(ji,jj) * vv(ji,jj,jk,Kmm)    & 
    359                      &              - e1u(ji  ,jj+1) * uu(ji  ,jj+1,jk,Kmm) + e1u(ji,jj) * uu(ji,jj,jk,Kmm)  ) * r1_e1e2f(ji,jj) 
    360                END DO 
    361             END DO 
    362          END DO 
     312         DO_3D_00_00( 1, jpkm1 ) 
     313            z3d(ji,jj,jk) = (   e2v(ji+1,jj  ) * vv(ji+1,jj  ,jk,Kmm) - e2v(ji,jj) * vv(ji,jj,jk,Kmm)    & 
     314               &              - e1u(ji  ,jj+1) * uu(ji  ,jj+1,jk,Kmm) + e1u(ji,jj) * uu(ji,jj,jk,Kmm)  ) * r1_e1e2f(ji,jj) 
     315         END_3D 
    363316         CALL lbc_lnk( 'diawri', z3d, 'F', 1. ) 
    364317         CALL iom_put( "relvor", z3d )                  ! relative vorticity 
    365318 
    366          DO jk = 1, jpkm1 
    367             DO jj = 1, jpj 
    368                DO ji = 1, jpi 
    369                   z3d(ji,jj,jk) = ff_f(ji,jj) + z3d(ji,jj,jk)  
    370                END DO 
    371             END DO 
    372          END DO 
     319         DO_3D_11_11( 1, jpkm1 ) 
     320            z3d(ji,jj,jk) = ff_f(ji,jj) + z3d(ji,jj,jk)  
     321         END_3D 
    373322         CALL iom_put( "absvor", z3d )                  ! absolute vorticity 
    374323 
    375          DO jk = 1, jpkm1 
    376             DO jj = 1, jpjm1 
    377                DO ji = 1, fs_jpim1   ! vector opt. 
    378                   ze3  = (  e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   & 
    379                      &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  ) 
    380                   IF( ze3 /= 0._wp ) THEN   ;   ze3 = 4._wp / ze3 
    381                   ELSE                      ;   ze3 = 0._wp 
    382                   ENDIF 
    383                   z3d(ji,jj,jk) = ze3 * z3d(ji,jj,jk)  
    384                END DO 
    385             END DO 
    386          END DO 
     324         DO_3D_00_00( 1, jpkm1 ) 
     325            ze3  = (  e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   & 
     326               &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  ) 
     327            IF( ze3 /= 0._wp ) THEN   ;   ze3 = 4._wp / ze3 
     328            ELSE                      ;   ze3 = 0._wp 
     329            ENDIF 
     330            z3d(ji,jj,jk) = ze3 * z3d(ji,jj,jk)  
     331         END_3D 
    387332         CALL lbc_lnk( 'diawri', z3d, 'F', 1. ) 
    388333         CALL iom_put( "potvor", z3d )                  ! potential vorticity 
    389334 
    390335      ENDIF 
    391     
    392336      ! 
    393337      IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN 
     
    404348      IF( iom_use("u_heattr") ) THEN 
    405349         z2d(:,:) = 0._wp  
    406          DO jk = 1, jpkm1 
    407             DO jj = 2, jpjm1 
    408                DO ji = fs_2, fs_jpim1   ! vector opt. 
    409                   z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) ) 
    410                END DO 
    411             END DO 
    412          END DO 
     350         DO_3D_00_00( 1, jpkm1 ) 
     351            z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) ) 
     352         END_3D 
    413353         CALL lbc_lnk( 'diawri', z2d, 'U', -1. ) 
    414354         CALL iom_put( "u_heattr", 0.5*rcp * z2d )    ! heat transport in i-direction 
     
    417357      IF( iom_use("u_salttr") ) THEN 
    418358         z2d(:,:) = 0.e0  
    419          DO jk = 1, jpkm1 
    420             DO jj = 2, jpjm1 
    421                DO ji = fs_2, fs_jpim1   ! vector opt. 
    422                   z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) ) 
    423                END DO 
    424             END DO 
    425          END DO 
     359         DO_3D_00_00( 1, jpkm1 ) 
     360            z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) ) 
     361         END_3D 
    426362         CALL lbc_lnk( 'diawri', z2d, 'U', -1. ) 
    427363         CALL iom_put( "u_salttr", 0.5 * z2d )        ! heat transport in i-direction 
     
    439375      IF( iom_use("v_heattr") ) THEN 
    440376         z2d(:,:) = 0.e0  
    441          DO jk = 1, jpkm1 
    442             DO jj = 2, jpjm1 
    443                DO ji = fs_2, fs_jpim1   ! vector opt. 
    444                   z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) ) 
    445                END DO 
    446             END DO 
    447          END DO 
     377         DO_3D_00_00( 1, jpkm1 ) 
     378            z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) ) 
     379         END_3D 
    448380         CALL lbc_lnk( 'diawri', z2d, 'V', -1. ) 
    449381         CALL iom_put( "v_heattr", 0.5*rcp * z2d )    !  heat transport in j-direction 
     
    452384      IF( iom_use("v_salttr") ) THEN 
    453385         z2d(:,:) = 0._wp  
    454          DO jk = 1, jpkm1 
    455             DO jj = 2, jpjm1 
    456                DO ji = fs_2, fs_jpim1   ! vector opt. 
    457                   z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) ) 
    458                END DO 
    459             END DO 
    460          END DO 
     386         DO_3D_00_00( 1, jpkm1 ) 
     387            z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) ) 
     388         END_3D 
    461389         CALL lbc_lnk( 'diawri', z2d, 'V', -1. ) 
    462390         CALL iom_put( "v_salttr", 0.5 * z2d )        !  heat transport in j-direction 
     
    465393      IF( iom_use("tosmint") ) THEN 
    466394         z2d(:,:) = 0._wp 
    467          DO jk = 1, jpkm1 
    468             DO jj = 2, jpjm1 
    469                DO ji = fs_2, fs_jpim1   ! vector opt. 
    470                   z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) *  ts(ji,jj,jk,jp_tem,Kmm) 
    471                END DO 
    472             END DO 
    473          END DO 
     395         DO_3D_00_00( 1, jpkm1 ) 
     396            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) *  ts(ji,jj,jk,jp_tem,Kmm) 
     397         END_3D 
    474398         CALL lbc_lnk( 'diawri', z2d, 'T', -1. ) 
    475399         CALL iom_put( "tosmint", rho0 * z2d )        ! Vertical integral of temperature 
     
    477401      IF( iom_use("somint") ) THEN 
    478402         z2d(:,:)=0._wp 
    479          DO jk = 1, jpkm1 
    480             DO jj = 2, jpjm1 
    481                DO ji = fs_2, fs_jpim1   ! vector opt. 
    482                   z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) 
    483                END DO 
    484             END DO 
    485          END DO 
     403         DO_3D_00_00( 1, jpkm1 ) 
     404            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) 
     405         END_3D 
    486406         CALL lbc_lnk( 'diawri', z2d, 'T', -1. ) 
    487407         CALL iom_put( "somint", rho0 * z2d )         ! Vertical integral of salinity 
     
    490410      CALL iom_put( "bn2", rn2 )                      ! Brunt-Vaisala buoyancy frequency (N^2) 
    491411      ! 
    492            
     412       
    493413      IF (ln_dia25h)   CALL dia_25h( kt, Kmm )        ! 25h averaging 
    494414 
     
    506426      INTEGER, DIMENSION(2) :: ierr 
    507427      !!---------------------------------------------------------------------- 
    508       ierr = 0 
    509       ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     & 
    510          &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     & 
    511          &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) ) 
     428      IF( nn_write == -1 ) THEN 
     429         dia_wri_alloc = 0 
     430      ELSE     
     431         ierr = 0 
     432         ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     & 
     433            &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     & 
     434            &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) ) 
    512435         ! 
    513       dia_wri_alloc = MAXVAL(ierr) 
    514       CALL mpp_sum( 'diawri', dia_wri_alloc ) 
     436         dia_wri_alloc = MAXVAL(ierr) 
     437         CALL mpp_sum( 'diawri', dia_wri_alloc ) 
     438         ! 
     439      ENDIF 
    515440      ! 
    516441   END FUNCTION dia_wri_alloc 
     442  
     443   INTEGER FUNCTION dia_wri_alloc_abl() 
     444      !!---------------------------------------------------------------------- 
     445     ALLOCATE(   ndex_hA(jpi*jpj), ndex_A (jpi*jpj*jpkam1), STAT=dia_wri_alloc_abl) 
     446      CALL mpp_sum( 'diawri', dia_wri_alloc_abl ) 
     447      ! 
     448   END FUNCTION dia_wri_alloc_abl 
    517449 
    518450    
     
    538470      INTEGER  ::   ierr                                     ! error code return from allocation 
    539471      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers 
     472      INTEGER  ::   ipka                                     ! ABL 
    540473      INTEGER  ::   jn, ierror                               ! local integers 
    541474      REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars 
     
    543476      REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace 
    544477      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace 
     478      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl   ! ABL 3D workspace 
    545479      !!---------------------------------------------------------------------- 
    546480      ! 
     
    576510      ijmi = 1      ;      ijma = jpj 
    577511      ipk = jpk 
     512      IF(ln_abl) ipka = jpkam1 
    578513 
    579514      ! define time axis 
     
    678613            &          "m", ipk, gdepw_1d, nz_W, "down" ) 
    679614 
     615         IF( ln_abl ) THEN  
     616         ! Define the ABL grid FILE ( nid_A ) 
     617            CALL dia_nam( clhstnam, nn_write, 'grid_ABL' ) 
     618            IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename 
     619            CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit 
     620               &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       & 
     621               &          nit000-1, zjulian, rn_Dt, nh_A, nid_A, domain_id=nidom, snc4chunks=snc4set ) 
     622            CALL histvert( nid_A, "ght_abl", "Vertical T levels",      &  ! Vertical grid: gdept 
     623               &           "m", ipka, ght_abl(2:jpka), nz_A, "up" ) 
     624            !                                                            ! Index of ocean points 
     625         ALLOCATE( zw3d_abl(jpi,jpj,ipka) )  
     626         zw3d_abl(:,:,:) = 1._wp  
     627         CALL wheneq( jpi*jpj*ipka, zw3d_abl, 1, 1., ndex_A , ndim_A  )      ! volume 
     628            CALL wheneq( jpi*jpj     , zw3d_abl, 1, 1., ndex_hA, ndim_hA )      ! surface 
     629         DEALLOCATE(zw3d_abl) 
     630         ENDIF 
    680631 
    681632         ! Declare all the output fields as NETCDF variables 
     
    727678         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm 
    728679            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    729 ! 
     680         ! 
     681         IF( ln_abl ) THEN 
     682            CALL histdef( nid_A, "t_abl", "Potential Temperature"     , "K"        ,       &  ! t_abl 
     683               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
     684            CALL histdef( nid_A, "q_abl", "Humidity"                  , "kg/kg"    ,       &  ! q_abl 
     685               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )  
     686            CALL histdef( nid_A, "u_abl", "Atmospheric U-wind   "     , "m/s"        ,     &  ! u_abl 
     687               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
     688            CALL histdef( nid_A, "v_abl", "Atmospheric V-wind   "     , "m/s"    ,         &  ! v_abl 
     689               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )  
     690            CALL histdef( nid_A, "tke_abl", "Atmospheric TKE   "     , "m2/s2"    ,        &  ! tke_abl 
     691               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )  
     692            CALL histdef( nid_A, "avm_abl", "Atmospheric turbulent viscosity", "m2/s"   ,  &  ! avm_abl 
     693               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )  
     694            CALL histdef( nid_A, "avt_abl", "Atmospheric turbulent diffusivity", "m2/s2",  &  ! avt_abl 
     695               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )  
     696            CALL histdef( nid_A, "pblh", "Atmospheric boundary layer height "  , "m",      &  ! pblh 
     697               &          jpi, jpj, nh_A,  1  , 1, 1   , -99 , 32, clop, zsto, zout )                  
     698#if defined key_si3 
     699            CALL histdef( nid_A, "oce_frac", "Fraction of open ocean"  , " ",      &  ! ato_i 
     700               &          jpi, jpj, nh_A,  1  , 1, 1   , -99 , 32, clop, zsto, zout ) 
     701#endif 
     702            CALL histend( nid_A, snc4chunks=snc4set ) 
     703         ENDIF 
     704         ! 
    730705         IF( ln_icebergs ) THEN 
    731706            CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , & 
     
    885860      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction    
    886861      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed    
    887 ! 
     862      ! 
     863      IF( ln_abl ) THEN  
     864         ALLOCATE( zw3d_abl(jpi,jpj,jpka) ) 
     865         IF( ln_mskland )   THEN  
     866            DO jk=1,jpka 
     867               zw3d_abl(:,:,jk) = tmask(:,:,1) 
     868            END DO        
     869         ELSE 
     870            zw3d_abl(:,:,:) = 1._wp      
     871         ENDIF        
     872         CALL histwrite( nid_A,  "pblh"   , it, pblh(:,:)                  *zw3d_abl(:,:,1     ), ndim_hA, ndex_hA )   ! pblh  
     873         CALL histwrite( nid_A,  "u_abl"  , it, u_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! u_abl 
     874         CALL histwrite( nid_A,  "v_abl"  , it, v_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! v_abl 
     875         CALL histwrite( nid_A,  "t_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! t_abl 
     876         CALL histwrite( nid_A,  "q_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! q_abl        
     877         CALL histwrite( nid_A,  "tke_abl", it, tke_abl (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! tke_abl 
     878         CALL histwrite( nid_A,  "avm_abl", it, avm_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avm_abl 
     879         CALL histwrite( nid_A,  "avt_abl", it, avt_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avt_abl  
     880#if defined key_si3 
     881         CALL histwrite( nid_A,  "oce_frac"   , it, ato_i(:,:)                                  , ndim_hA, ndex_hA )   ! ato_i 
     882#endif 
     883         DEALLOCATE(zw3d_abl) 
     884      ENDIF 
     885      ! 
    888886      IF( ln_icebergs ) THEN 
    889887         ! 
     
    931929      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress 
    932930 
    933       CALL histwrite( nid_W, "vovecrtz", it, ww             , ndim_T, ndex_T )    ! vert. current 
     931      IF( ln_zad_Aimp ) THEN 
     932         CALL histwrite( nid_W, "vovecrtz", it, ww + wi     , ndim_T, ndex_T )    ! vert. current 
     933      ELSE 
     934         CALL histwrite( nid_W, "vovecrtz", it, ww          , ndim_T, ndex_T )    ! vert. current 
     935      ENDIF 
    934936      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef. 
    935937      CALL histwrite( nid_W, "votkeavm", it, avm            , ndim_T, ndex_T )    ! T vert. eddy visc. coef. 
     
    951953         CALL histclo( nid_V ) 
    952954         CALL histclo( nid_W ) 
     955         IF(ln_abl) CALL histclo( nid_A ) 
    953956      ENDIF 
    954957      ! 
     
    974977      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created 
    975978      !! 
    976       INTEGER :: inum 
     979      INTEGER :: inum, jk 
    977980      !!---------------------------------------------------------------------- 
    978981      !  
     
    981984      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created ' 
    982985      IF(lwp) WRITE(numout,*) '                and named :', cdfile_name, '...nc' 
    983  
    984 #if defined key_si3 
    985      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl ) 
    986 #else 
    987      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. ) 
    988 #endif 
    989  
     986      ! 
     987      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. ) 
     988      ! 
    990989      CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) )    ! now temperature 
    991990      CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) )    ! now salinity 
     
    993992      CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm)                )    ! now i-velocity 
    994993      CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm)                )    ! now j-velocity 
    995       CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww                )    ! now k-velocity 
     994      IF( ln_zad_Aimp ) THEN 
     995         CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww + wi        )    ! now k-velocity 
     996      ELSE 
     997         CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww             )    ! now k-velocity 
     998      ENDIF 
     999      CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep            )    ! now k-velocity 
     1000      CALL iom_rstput( 0, 0, inum, 'ht'     , ht                 )    ! now water column height 
     1001      ! 
     1002      IF ( ln_isf ) THEN 
     1003         IF (ln_isfcav_mlt) THEN 
     1004            CALL iom_rstput( 0, 0, inum, 'fwfisf_cav', fwfisf_cav          )    ! now k-velocity 
     1005            CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav    )    ! now k-velocity 
     1006            CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav    )    ! now k-velocity 
     1007            CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,wp) )    ! now k-velocity 
     1008            CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,wp) )    ! now k-velocity 
     1009            CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,wp), ktype = jp_i1 ) 
     1010         END IF 
     1011         IF (ln_isfpar_mlt) THEN 
     1012            CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,wp) )    ! now k-velocity 
     1013            CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par          )    ! now k-velocity 
     1014            CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par    )    ! now k-velocity 
     1015            CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par    )    ! now k-velocity 
     1016            CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,wp) )    ! now k-velocity 
     1017            CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,wp) )    ! now k-velocity 
     1018            CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,wp), ktype = jp_i1 ) 
     1019         END IF 
     1020      END IF 
     1021      ! 
    9961022      IF( ALLOCATED(ahtu) ) THEN 
    9971023         CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point 
     
    10171043         CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd            )    ! now StokesDrift k-velocity 
    10181044      ENDIF 
    1019   
     1045      IF ( ln_abl ) THEN 
     1046         CALL iom_rstput ( 0, 0, inum, "uz1_abl",   u_abl(:,:,2,nt_a  ) )   ! now first level i-wind 
     1047         CALL iom_rstput ( 0, 0, inum, "vz1_abl",   v_abl(:,:,2,nt_a  ) )   ! now first level j-wind 
     1048         CALL iom_rstput ( 0, 0, inum, "tz1_abl",  tq_abl(:,:,2,nt_a,1) )   ! now first level temperature 
     1049         CALL iom_rstput ( 0, 0, inum, "qz1_abl",  tq_abl(:,:,2,nt_a,2) )   ! now first level humidity 
     1050      ENDIF 
     1051      ! 
     1052      CALL iom_close( inum ) 
     1053      !  
    10201054#if defined key_si3 
    10211055      IF( nn_ice == 2 ) THEN   ! condition needed in case agrif + ice-model but no-ice in child grid 
     1056         CALL iom_open( TRIM(cdfile_name)//'_ice', inum, ldwrt = .TRUE., kdlev = jpl, cdcomp = 'ICE' ) 
    10221057         CALL ice_wri_state( inum ) 
     1058         CALL iom_close( inum ) 
    10231059      ENDIF 
    10241060#endif 
    1025       ! 
    1026       CALL iom_close( inum ) 
    1027       !  
     1061 
    10281062   END SUBROUTINE dia_wri_state 
    10291063 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/CANAL/MY_SRC/domvvl.F90

    r12489 r13197  
    3737 
    3838   PUBLIC  dom_vvl_init       ! called by domain.F90 
     39   PUBLIC  dom_vvl_zgr        ! called by isfcpl.F90 
    3940   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90 
    4041   PUBLIC  dom_vvl_sf_update  ! called by step.F90 
     
    6263   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! retoring period for low freq. divergence 
    6364 
     65   !! * Substitutions 
     66#  include "do_loop_substitute.h90" 
    6467   !!---------------------------------------------------------------------- 
    6568   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    116119      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 
    117120      ! 
     121      IF(lwp) WRITE(numout,*) 
     122      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 
     123      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
     124      ! 
     125      CALL dom_vvl_ctl     ! choose vertical coordinate (z_star, z_tilde or layer) 
     126      ! 
     127      !                    ! Allocate module arrays 
     128      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 
     129      ! 
     130      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 
     131      CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' ) 
     132      e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all 
     133      ! 
     134      CALL dom_vvl_zgr(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column 
     135      ! 
     136   END SUBROUTINE dom_vvl_init 
     137   ! 
     138   SUBROUTINE dom_vvl_zgr(Kbb, Kmm, Kaa) 
     139      !!---------------------------------------------------------------------- 
     140      !!                ***  ROUTINE dom_vvl_init  *** 
     141      !!                    
     142      !! ** Purpose :  Interpolation of all scale factors,  
     143      !!               depths and water column heights 
     144      !! 
     145      !! ** Method  :  - interpolate scale factors 
     146      !! 
     147      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b) 
     148      !!              - Regrid: e3(u/v)_n 
     149      !!                        e3(u/v)_b        
     150      !!                        e3w_n            
     151      !!                        e3(u/v)w_b       
     152      !!                        e3(u/v)w_n       
     153      !!                        gdept_n, gdepw_n and gde3w_n 
     154      !!              - h(t/u/v)_0 
     155      !!              - frq_rst_e3t and frq_rst_hdv 
     156      !! 
     157      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling. 
     158      !!---------------------------------------------------------------------- 
     159      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 
     160      !!---------------------------------------------------------------------- 
    118161      INTEGER ::   ji, jj, jk 
    119162      INTEGER ::   ii0, ii1, ij0, ij1 
    120163      REAL(wp)::   zcoef 
    121164      !!---------------------------------------------------------------------- 
    122       ! 
    123       IF(lwp) WRITE(numout,*) 
    124       IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 
    125       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    126       ! 
    127       CALL dom_vvl_ctl     ! choose vertical coordinate (z_star, z_tilde or layer) 
    128       ! 
    129       !                    ! Allocate module arrays 
    130       IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 
    131       ! 
    132       !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 
    133       CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' ) 
    134       e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all 
    135165      ! 
    136166      !                    !== Set of all other vertical scale factors  ==!  (now and before) 
     
    160190      gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb) 
    161191      gdepw(:,:,1,Kbb) = 0.0_wp 
    162       DO jk = 2, jpk                               ! vertical sum 
    163          DO jj = 1,jpj 
    164             DO ji = 1,jpi 
    165                !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
    166                !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
    167                !                             ! 0.5 where jk = mikt      
     192      DO_3D_11_11( 2, jpk ) 
     193         !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
     194         !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
     195         !                             ! 0.5 where jk = mikt      
    168196!!gm ???????   BUG ?  gdept(:,:,:,Kmm) as well as gde3w  does not include the thickness of ISF ?? 
    169                zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 
    170                gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
    171                gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  & 
    172                   &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm))  
    173                gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
    174                gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 
    175                gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  & 
    176                   &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb))  
    177             END DO 
    178          END DO 
    179       END DO 
     197         zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 
     198         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
     199         gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  & 
     200            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm))  
     201         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
     202         gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 
     203         gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  & 
     204            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb))  
     205      END_3D 
    180206      ! 
    181207      !                    !==  thickness of the water column  !!   (ocean portion only) 
     
    212238         ENDIF 
    213239         IF ( ln_vvl_zstar_at_eqtor ) THEN   ! use z-star in vicinity of the Equator 
    214             DO jj = 1, jpj 
    215                DO ji = 1, jpi 
     240            DO_2D_11_11 
    216241!!gm  case |gphi| >= 6 degrees is useless   initialized just above by default 
    217                   IF( ABS(gphit(ji,jj)) >= 6.) THEN 
    218                      ! values outside the equatorial band and transition zone (ztilde) 
    219                      frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp ) 
    220                      frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 
    221                   ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star 
    222                      ! values inside the equatorial band (ztilde as zstar) 
    223                      frq_rst_e3t(ji,jj) =  0.0_wp 
    224                      frq_rst_hdv(ji,jj) =  1.0_wp / rn_Dt 
    225                   ELSE                                      ! transition band (2.5 to 6 degrees N/S) 
    226                      !                                      ! (linearly transition from z-tilde to z-star) 
    227                      frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   & 
    228                         &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
    229                         &                                          * 180._wp / 3.5_wp ) ) 
    230                      frq_rst_hdv(ji,jj) = (1.0_wp / rn_Dt)                                & 
    231                         &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rn_Dt) )*0.5_wp   & 
    232                         &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
    233                         &                                          * 180._wp / 3.5_wp ) ) 
    234                   ENDIF 
    235                END DO 
    236             END DO 
     242               IF( ABS(gphit(ji,jj)) >= 6.) THEN 
     243                  ! values outside the equatorial band and transition zone (ztilde) 
     244                  frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp ) 
     245                  frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 
     246               ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star 
     247                  ! values inside the equatorial band (ztilde as zstar) 
     248                  frq_rst_e3t(ji,jj) =  0.0_wp 
     249                  frq_rst_hdv(ji,jj) =  1.0_wp / rn_Dt 
     250               ELSE                                      ! transition band (2.5 to 6 degrees N/S) 
     251                  !                                      ! (linearly transition from z-tilde to z-star) 
     252                  frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   & 
     253                     &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
     254                     &                                          * 180._wp / 3.5_wp ) ) 
     255                  frq_rst_hdv(ji,jj) = (1.0_wp / rn_Dt)                                & 
     256                     &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rn_Dt) )*0.5_wp   & 
     257                     &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
     258                     &                                          * 180._wp / 3.5_wp ) ) 
     259               ENDIF 
     260            END_2D 
    237261            IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN 
    238262               IF( nn_cfg == 3 ) THEN   ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 
     
    264288      ENDIF 
    265289      ! 
    266    END SUBROUTINE dom_vvl_init 
     290   END SUBROUTINE dom_vvl_zgr 
    267291 
    268292 
     
    329353      END DO 
    330354      ! 
    331       IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate ! 
    332          !                                                            ! ------baroclinic part------ ! 
     355      IF( (ln_vvl_ztilde .OR. ln_vvl_layer) .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate ! 
     356         !                                                               ! ------baroclinic part------ ! 
    333357         ! I - initialization 
    334358         ! ================== 
     
    383407         zwu(:,:) = 0._wp 
    384408         zwv(:,:) = 0._wp 
    385          DO jk = 1, jpkm1        ! a - first derivative: diffusive fluxes 
    386             DO jj = 1, jpjm1 
    387                DO ji = 1, fs_jpim1   ! vector opt. 
    388                   un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           & 
    389                      &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
    390                   vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           &  
    391                      &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
    392                   zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
    393                   zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 
    394                END DO 
    395             END DO 
    396          END DO 
    397          DO jj = 1, jpj          ! b - correction for last oceanic u-v points 
    398             DO ji = 1, jpi 
    399                un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 
    400                vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 
    401             END DO 
    402          END DO 
    403          DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes 
    404             DO jj = 2, jpjm1 
    405                DO ji = fs_2, fs_jpim1   ! vector opt. 
    406                   tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
    407                      &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    & 
    408                      &                                            ) * r1_e1e2t(ji,jj) 
    409                END DO 
    410             END DO 
    411          END DO 
     409         DO_3D_10_10( 1, jpkm1 ) 
     410            un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           & 
     411               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
     412            vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           &  
     413               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
     414            zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
     415            zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 
     416         END_3D 
     417         DO_2D_11_11 
     418            un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 
     419            vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 
     420         END_2D 
     421         DO_3D_00_00( 1, jpkm1 ) 
     422            tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
     423               &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    & 
     424               &                                            ) * r1_e1e2t(ji,jj) 
     425         END_3D 
    412426         !                       ! d - thickness diffusion transport: boundary conditions 
    413427         !                             (stored for tracer advction and continuity equation) 
     
    416430         ! 4 - Time stepping of baroclinic scale factors 
    417431         ! --------------------------------------------- 
    418          ! Leapfrog time stepping 
    419          ! ~~~~~~~~~~~~~~~~~~~~~~ 
    420432         CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp ) 
    421433         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + rDt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 
     
    613625         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 
    614626      ENDIF 
    615       gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm) 
    616       gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm) 
    617  
    618       e3t(:,:,:,Kmm) = e3t(:,:,:,Kaa) 
    619       e3u(:,:,:,Kmm) = e3u(:,:,:,Kaa) 
    620       e3v(:,:,:,Kmm) = e3v(:,:,:,Kaa) 
    621627 
    622628      ! Compute all missing vertical scale factor and depths 
     
    641647      gdepw(:,:,1,Kmm) = 0.0_wp 
    642648      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 
    643       DO jk = 2, jpk 
    644          DO jj = 1,jpj 
    645             DO ji = 1,jpi 
    646               !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
    647                                                                  ! 1 for jk = mikt 
    648                zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    649                gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
    650                gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  & 
    651                    &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) )  
    652                gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
    653             END DO 
    654          END DO 
    655       END DO 
     649      DO_3D_11_11( 2, jpk ) 
     650        !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
     651                                                           ! 1 for jk = mikt 
     652         zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
     653         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
     654         gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  & 
     655             &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) )  
     656         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
     657      END_3D 
    656658 
    657659      ! Local depth and Inverse of the local depth of the water 
     
    700702         ! 
    701703      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean 
    702          DO jk = 1, jpk 
    703             DO jj = 1, jpjm1 
    704                DO ji = 1, fs_jpim1   ! vector opt. 
    705                   pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   & 
    706                      &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
    707                      &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
    708                END DO 
    709             END DO 
    710          END DO 
     704         DO_3D_10_10( 1, jpk ) 
     705            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   & 
     706               &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
     707               &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
     708         END_3D 
    711709         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp ) 
    712710         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 
    713711         ! 
    714712      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean 
    715          DO jk = 1, jpk 
    716             DO jj = 1, jpjm1 
    717                DO ji = 1, fs_jpim1   ! vector opt. 
    718                   pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   & 
    719                      &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
    720                      &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
    721                END DO 
    722             END DO 
    723          END DO 
     713         DO_3D_10_10( 1, jpk ) 
     714            pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   & 
     715               &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
     716               &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
     717         END_3D 
    724718         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp ) 
    725719         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 
    726720         ! 
    727721      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean 
    728          DO jk = 1, jpk 
    729             DO jj = 1, jpjm1 
    730                DO ji = 1, fs_jpim1   ! vector opt. 
    731                   pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
    732                      &                       *    r1_e1e2f(ji,jj)                                                  & 
    733                      &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
    734                      &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
    735                END DO 
    736             END DO 
    737          END DO 
     722         DO_3D_10_10( 1, jpk ) 
     723            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
     724               &                       *    r1_e1e2f(ji,jj)                                                  & 
     725               &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
     726               &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
     727         END_3D 
    738728         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp ) 
    739729         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 
     
    810800            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 
    811801            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 
     802            ! 
    812803            !                             ! --------- ! 
    813804            !                             ! all cases ! 
    814805            !                             ! --------- ! 
     806            ! 
    815807            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist 
    816808               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 
     
    828820               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files' 
    829821               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 
    830                IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.' 
     822               IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    831823               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 
    832824               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 
     
    835827               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files' 
    836828               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 
    837                IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.' 
     829               IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    838830               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 
    839831               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
     
    842834               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file' 
    843835               IF(lwp) write(numout,*) 'Compute scale factor from sshn' 
    844                IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.' 
     836               IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    845837               DO jk = 1, jpk 
    846838                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 
     
    895887                  ssh(:,:,Kbb) = -ssh_ref 
    896888 
    897                   DO jj = 1, jpj 
    898                      DO ji = 1, jpi 
    899                         IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth 
    900                            ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 
    901                            ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 
    902                         ENDIF 
    903                      ENDDO 
    904                   ENDDO 
     889                  DO_2D_11_11 
     890                     IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth 
     891                        ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 
     892                        ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 
     893                     ENDIF 
     894                  END_2D 
    905895               ENDIF !If test case else 
    906896 
     
    913903               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    914904 
    915                DO ji = 1, jpi 
    916                   DO jj = 1, jpj 
    917                      IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 
    918                        CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 
    919                      ENDIF 
    920                   END DO  
    921                END DO  
     905               DO_2D_11_11 
     906                  IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 
     907                     CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 
     908                  ENDIF 
     909               END_2D 
    922910               ! 
    923911            ELSE 
    924912               ! 
    925                ! usr_def_istate called here only to get sshb, that is needed to initialize e3t(Kbb) and e3t(Kmm) 
    926                CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
    927                ! usr_def_istate will be called again in istate_init to initialize ts(bn), ssh(bn), u(bn) and v(bn) 
     913               ! usr_def_istate called here only to get ssh(Kbb) needed to initialize e3t(Kbb) and e3t(Kmm) 
     914               ! 
     915               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )   
     916               ! 
     917               ! usr_def_istate will be called again in istate_init to initialize ts, ssh, u and v 
    928918               ! 
    929919               DO jk=1,jpk 
    930                   e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) & 
    931                     &                              / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
    932                     &              + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )   ! make sure e3t_b != 0 on land points 
     920                  e3t(:,:,jk,Kbb) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) & 
     921                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
     922                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )   ! make sure e3t(:,:,:,Kbb) != 0 on land points 
    933923               END DO 
    934924               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 
    935                ssh(:,:  ,Kmm) = ssh(:,:  ,Kbb)   ! needed later for gde3w 
     925               ssh(:,:,Kmm) = ssh(:,:,Kbb)                                     ! needed later for gde3w 
    936926               ! 
    937927            END IF           ! end of ll_wd edits 
     
    10251015      ! 
    10261016      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 
    1027       IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 
    10281017      ! 
    10291018      IF(lwp) THEN                   ! Print the choice 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/CANAL/MY_SRC/stpctl.F90

    r13159 r13197  
    1919   USE dom_oce         ! ocean space and time domain variables  
    2020   USE c1d             ! 1D vertical configuration 
    21    USE zdf_oce ,  ONLY : ln_zad_Aimp       ! ocean vertical physics variables 
    22    USE wet_dry,   ONLY : ll_wd, ssh_ref    ! reference depth for negative bathy 
    23    !   
    2421   USE diawri          ! Standard run outputs       (dia_wri_state routine) 
     22   ! 
    2523   USE in_out_manager  ! I/O manager 
    2624   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2725   USE lib_mpp         ! distributed memory computing 
    28    ! 
     26   USE zdf_oce ,  ONLY : ln_zad_Aimp       ! ocean vertical physics variables 
     27   USE wet_dry,   ONLY : ll_wd, ssh_ref    ! reference depth for negative bathy 
     28 
    2929   USE netcdf          ! NetCDF library 
    3030   IMPLICIT NONE 
     
    3333   PUBLIC stp_ctl           ! routine called by step.F90 
    3434 
    35    INTEGER                ::   nrunid   ! netcdf file id 
    36    INTEGER, DIMENSION(8)  ::   nvarid   ! netcdf variable id 
     35   INTEGER  ::   idrun, idtime, idssh, idu, ids1, ids2, idt1, idt2, idc1, idw1, istatus 
     36   LOGICAL  ::   lsomeoce 
    3737   !!---------------------------------------------------------------------- 
    3838   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    4242CONTAINS 
    4343 
    44    SUBROUTINE stp_ctl( kt, Kmm ) 
     44   SUBROUTINE stp_ctl( kt, Kbb, Kmm, kindic ) 
    4545      !!---------------------------------------------------------------------- 
    4646      !!                    ***  ROUTINE stp_ctl  *** 
     
    5050      !! ** Method  : - Save the time step in numstp 
    5151      !!              - Print it each 50 time steps 
    52       !!              - Stop the run IF problem encountered by setting nstop > 0 
     52      !!              - Stop the run IF problem encountered by setting indic=-3 
    5353      !!                Problems checked: |ssh| maximum larger than 10 m 
    5454      !!                                  |U|   maximum larger than 10 m/s  
     
    5757      !! ** Actions :   "time.step" file = last ocean time-step 
    5858      !!                "run.stat"  file = run statistics 
    59       !!                 nstop indicator sheared among all local domain 
     59      !!                nstop indicator sheared among all local domain (lk_mpp=T) 
    6060      !!---------------------------------------------------------------------- 
    6161      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index 
    62       INTEGER, INTENT(in   ) ::   Kmm      ! ocean time level index 
     62      INTEGER, INTENT(in   ) ::   Kbb, Kmm      ! ocean time level index 
     63      INTEGER, INTENT(inout) ::   kindic   ! error indicator 
    6364      !! 
    64       INTEGER                         ::   ji                                    ! dummy loop indices 
    65       INTEGER                         ::   idtime, istatus 
    66       INTEGER , DIMENSION(9)          ::   iareasum, iareamin, iareamax 
    67       INTEGER , DIMENSION(3,4)        ::   iloc                                  ! min/max loc indices 
    68       REAL(wp)                        ::   zzz                                   ! local real  
    69       REAL(wp), DIMENSION(9)          ::   zmax, zmaxlocal 
    70       LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns 
    71       LOGICAL, DIMENSION(jpi,jpj,jpk) ::   llmsk 
    72       CHARACTER(len=20)               ::   clname 
     65      INTEGER                ::   ji, jj, jk          ! dummy loop indices 
     66      INTEGER, DIMENSION(2)  ::   ih                  ! min/max loc indices 
     67      INTEGER, DIMENSION(3)  ::   iu, is1, is2        ! min/max loc indices 
     68      REAL(wp)               ::   zzz                 ! local real  
     69      REAL(wp), DIMENSION(9) ::   zmax 
     70      LOGICAL                ::   ll_wrtstp, ll_colruns, ll_wrtruns 
     71      CHARACTER(len=20) :: clname 
    7372      !!---------------------------------------------------------------------- 
    74       IF( nstop > 0 .AND. ngrdstop > -1 )   RETURN   !   stpctl was already called by a child grid 
    75       ! 
    76       ll_wrtstp  = ( MOD( kt-nit000, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend ) 
    77       ll_colruns = ll_wrtstp .AND. sn_cfctl%l_runstat .AND. jpnij > 1  
    78       ll_wrtruns = ( ll_colruns .OR. jpnij == 1 ) .AND. lwm 
    79       ! 
    80       IF( kt == nit000 ) THEN 
    81          ! 
    82          IF( lwp ) THEN 
    83             WRITE(numout,*) 
    84             WRITE(numout,*) 'stp_ctl : time-stepping control' 
    85             WRITE(numout,*) '~~~~~~~' 
    86          ENDIF 
    87          !                                ! open time.step    ascii file, done only by 1st subdomain 
    88          IF( lwm )   CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 
    89          ! 
    90          IF( ll_wrtruns ) THEN 
    91             !                             ! open run.stat     ascii file, done only by 1st subdomain 
     73      ! 
     74      ll_wrtstp  = ( MOD( kt, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend ) 
     75      ll_colruns = ll_wrtstp .AND. ( sn_cfctl%l_runstat ) 
     76      ll_wrtruns = ll_colruns .AND. lwm 
     77      IF( kt == nit000 .AND. lwp ) THEN 
     78         WRITE(numout,*) 
     79         WRITE(numout,*) 'stp_ctl : time-stepping control' 
     80         WRITE(numout,*) '~~~~~~~' 
     81         !                                ! open time.step file 
     82         IF( lwm ) CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 
     83         !                                ! open run.stat file(s) at start whatever 
     84         !                                ! the value of sn_cfctl%ptimincr 
     85         IF( lwm .AND. ( sn_cfctl%l_runstat ) ) THEN 
    9286            CALL ctl_opn( numrun, 'run.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 
    93             !                             ! open run.stat.nc netcdf file, done only by 1st subdomain 
    9487            clname = 'run.stat.nc' 
    9588            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) 
    96             istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, nrunid ) 
    97             istatus = NF90_DEF_DIM( nrunid, 'time', NF90_UNLIMITED, idtime ) 
    98             istatus = NF90_DEF_VAR( nrunid, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), nvarid(1) ) 
    99             istatus = NF90_DEF_VAR( nrunid,   'abs_u_max', NF90_DOUBLE, (/ idtime /), nvarid(2) ) 
    100             istatus = NF90_DEF_VAR( nrunid,       's_min', NF90_DOUBLE, (/ idtime /), nvarid(3) ) 
    101             istatus = NF90_DEF_VAR( nrunid,       's_max', NF90_DOUBLE, (/ idtime /), nvarid(4) ) 
    102             istatus = NF90_DEF_VAR( nrunid,       't_min', NF90_DOUBLE, (/ idtime /), nvarid(5) ) 
    103             istatus = NF90_DEF_VAR( nrunid,       't_max', NF90_DOUBLE, (/ idtime /), nvarid(6) ) 
     89            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, idrun ) 
     90            istatus = NF90_DEF_DIM( idrun, 'time', NF90_UNLIMITED, idtime ) 
     91            istatus = NF90_DEF_VAR( idrun, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), idssh ) 
     92            istatus = NF90_DEF_VAR( idrun,   'abs_u_max', NF90_DOUBLE, (/ idtime /), idu  ) 
     93            istatus = NF90_DEF_VAR( idrun,       's_min', NF90_DOUBLE, (/ idtime /), ids1 ) 
     94            istatus = NF90_DEF_VAR( idrun,       's_max', NF90_DOUBLE, (/ idtime /), ids2 ) 
     95            istatus = NF90_DEF_VAR( idrun,       't_min', NF90_DOUBLE, (/ idtime /), idt1 ) 
     96            istatus = NF90_DEF_VAR( idrun,       't_max', NF90_DOUBLE, (/ idtime /), idt2 ) 
    10497            IF( ln_zad_Aimp ) THEN 
    105                istatus = NF90_DEF_VAR( nrunid,   'Cf_max', NF90_DOUBLE, (/ idtime /), nvarid(7) ) 
    106                istatus = NF90_DEF_VAR( nrunid,'abs_wi_max',NF90_DOUBLE, (/ idtime /), nvarid(8) ) 
     98               istatus = NF90_DEF_VAR( idrun,   'abs_wi_max', NF90_DOUBLE, (/ idtime /), idw1 ) 
     99               istatus = NF90_DEF_VAR( idrun,       'Cf_max', NF90_DOUBLE, (/ idtime /), idc1 ) 
    107100            ENDIF 
    108             istatus = NF90_ENDDEF(nrunid) 
    109          ENDIF 
    110          !     
    111       ENDIF 
    112       ! 
    113       !                                   !==              write current time step              ==! 
    114       !                                   !==  done only by 1st subdomain at writting timestep  ==! 
    115       IF( lwm .AND. ll_wrtstp ) THEN 
     101            istatus = NF90_ENDDEF(idrun) 
     102            zmax(8:9) = 0._wp    ! initialise to zero in case ln_zad_Aimp option is not in use 
     103         ENDIF 
     104      ENDIF 
     105      IF( kt == nit000 )   lsomeoce = COUNT( ssmask(:,:) == 1._wp ) > 0 
     106      ! 
     107      IF(lwm .AND. ll_wrtstp) THEN        !==  current time step  ==!   ("time.step" file) 
    116108         WRITE ( numstp, '(1x, i8)' )   kt 
    117109         REWIND( numstp ) 
    118110      ENDIF 
    119       !                                   !==            test of local extrema           ==! 
    120       !                                   !==  done by all processes at every time step  ==! 
    121       ! 
    122       ! define zmax default value. needed for land processors 
    123       IF( ll_colruns ) THEN    ! default value: must not be kept when calling mpp_max -> must be as small as possible 
    124          zmax(:) = -HUGE(1._wp) 
    125       ELSE                     ! default value: must not give true for any of the tests bellow (-> avoid manipulating HUGE...) 
    126          zmax(:) =  0._wp 
    127          zmax(3) = -1._wp      ! avoid salinity minimum at 0. 
    128       ENDIF 
    129       ! 
    130       llmsk(:,:,1) = ssmask(:,:) == 1._wp 
    131       IF( COUNT( llmsk(:,:,1) ) > 0 ) THEN   ! avoid huge values sent back for land processors... 
    132          IF( ll_wd ) THEN 
    133             zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref ), mask = llmsk(:,:,1) )   ! ssh max 
    134          ELSE 
    135             zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm)           ), mask = llmsk(:,:,1) )   ! ssh max 
    136          ENDIF 
    137       ENDIF 
    138       zmax(2) = MAXVAL( ABS( uu(:,:,:,Kmm) ) )                                       ! velocity max (zonal only) 
    139       llmsk(:,:,:) = tmask(:,:,:) == 1._wp 
    140       IF( COUNT( llmsk(:,:,:) ) > 0 ) THEN   ! avoid huge values sent back for land processors... 
    141          zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     ! minus salinity max 
    142          zmax(4) = MAXVAL(  ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     !       salinity max 
    143          IF( ll_colruns .OR. jpnij == 1 ) THEN     ! following variables are used only in the netcdf file 
    144             zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  ! minus temperature max 
    145             zmax(6) = MAXVAL(  ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  !       temperature max 
    146             IF( ln_zad_Aimp ) THEN 
    147                zmax(7) = MAXVAL(   Cu_adv(:,:,:)   , mask = llmsk )                  ! partitioning coeff. max 
    148                llmsk(:,:,:) = wmask(:,:,:) == 1._wp 
    149                IF( COUNT( llmsk(:,:,:) ) > 0 ) THEN   ! avoid huge values sent back for land processors... 
    150                   zmax(8) = MAXVAL(ABS( wi(:,:,:) ), mask = llmsk )                  ! implicit vertical vel. max 
    151                ENDIF 
    152             ENDIF 
    153          ENDIF 
    154       ENDIF 
    155       zmax(9) = REAL( nstop, wp )                                              ! stop indicator 
    156       !                                   !==               get global extrema             ==! 
    157       !                                   !==  done by all processes if writting run.stat  ==! 
     111      ! 
     112      !                                   !==  test of extrema  ==! 
     113      IF( ll_wd ) THEN 
     114         zmax(1) = MAXVAL(  ABS( ssh(:,:,Kmm) + ssh_ref*tmask(:,:,1) )  )        ! ssh max  
     115      ELSE 
     116         zmax(1) = MAXVAL(  ABS( ssh(:,:,Kmm) )  )                               ! ssh max 
     117      ENDIF 
     118      zmax(2) = MAXVAL(  ABS( uu(:,:,:,Kmm) )  )                                  ! velocity max (zonal only) 
     119      zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp )   ! minus salinity max 
     120      zmax(4) = MAXVAL(  ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp )   !       salinity max 
     121      zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp )   ! minus temperature max 
     122      zmax(6) = MAXVAL(  ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp )   !       temperature max 
     123      zmax(7) = REAL( nstop , wp )                                            ! stop indicator 
     124      IF( ln_zad_Aimp ) THEN 
     125         zmax(8) = MAXVAL(  ABS( wi(:,:,:) ) , mask = wmask(:,:,:) == 1._wp ) ! implicit vertical vel. max 
     126         zmax(9) = MAXVAL(   Cu_adv(:,:,:)   , mask = tmask(:,:,:) == 1._wp ) ! partitioning coeff. max 
     127      ENDIF 
     128      ! 
    158129      IF( ll_colruns ) THEN 
    159          zmaxlocal(:) = zmax(:) 
    160130         CALL mpp_max( "stpctl", zmax )          ! max over the global domain 
    161          nstop = NINT( zmax(9) )                 ! update nstop indicator (now sheared among all local domains) 
    162       ENDIF 
    163       !                                   !==              write "run.stat" files              ==! 
    164       !                                   !==  done only by 1st subdomain at writting timestep  ==! 
     131         nstop = NINT( zmax(7) )                 ! nstop indicator sheared among all local domains 
     132      ENDIF 
     133      !                                   !==  run statistics  ==!   ("run.stat" files) 
    165134      IF( ll_wrtruns ) THEN 
    166135         WRITE(numrun,9500) kt, zmax(1), zmax(2), -zmax(3), zmax(4) 
    167          istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ zmax(1)/), (/kt/), (/1/) ) 
    168          istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/ zmax(2)/), (/kt/), (/1/) ) 
    169          istatus = NF90_PUT_VAR( nrunid, nvarid(3), (/-zmax(3)/), (/kt/), (/1/) ) 
    170          istatus = NF90_PUT_VAR( nrunid, nvarid(4), (/ zmax(4)/), (/kt/), (/1/) ) 
    171          istatus = NF90_PUT_VAR( nrunid, nvarid(5), (/-zmax(5)/), (/kt/), (/1/) ) 
    172          istatus = NF90_PUT_VAR( nrunid, nvarid(6), (/ zmax(6)/), (/kt/), (/1/) ) 
     136         istatus = NF90_PUT_VAR( idrun, idssh, (/ zmax(1)/), (/kt/), (/1/) ) 
     137         istatus = NF90_PUT_VAR( idrun,   idu, (/ zmax(2)/), (/kt/), (/1/) ) 
     138         istatus = NF90_PUT_VAR( idrun,  ids1, (/-zmax(3)/), (/kt/), (/1/) ) 
     139         istatus = NF90_PUT_VAR( idrun,  ids2, (/ zmax(4)/), (/kt/), (/1/) ) 
     140         istatus = NF90_PUT_VAR( idrun,  idt1, (/-zmax(5)/), (/kt/), (/1/) ) 
     141         istatus = NF90_PUT_VAR( idrun,  idt2, (/ zmax(6)/), (/kt/), (/1/) ) 
    173142         IF( ln_zad_Aimp ) THEN 
    174             istatus = NF90_PUT_VAR( nrunid, nvarid(7), (/ zmax(7)/), (/kt/), (/1/) ) 
    175             istatus = NF90_PUT_VAR( nrunid, nvarid(8), (/ zmax(8)/), (/kt/), (/1/) ) 
    176          ENDIF 
    177          IF( kt == nitend )   istatus = NF90_CLOSE(nrunid) 
     143            istatus = NF90_PUT_VAR( idrun,  idw1, (/ zmax(8)/), (/kt/), (/1/) ) 
     144            istatus = NF90_PUT_VAR( idrun,  idc1, (/ zmax(9)/), (/kt/), (/1/) ) 
     145         ENDIF 
     146         IF( MOD( kt , 100 ) == 0 ) istatus = NF90_SYNC(idrun) 
     147         IF( kt == nitend         ) istatus = NF90_CLOSE(idrun) 
    178148      END IF 
    179       !                                   !==               error handling               ==! 
    180       !                                   !==  done by all processes at every time step  ==! 
    181       ! 
    182       IF(   zmax(1) >   20._wp .OR.   &                   ! too large sea surface height ( > 20 m ) 
    183          &  zmax(2) >   10._wp .OR.   &                   ! too large velocity ( > 10 m/s) 
    184 !!$         &  zmax(3) >=   0._wp .OR.   &                   ! negative or zero sea surface salinity 
    185 !!$         &  zmax(4) >= 100._wp .OR.   &                   ! too large sea surface salinity ( > 100 ) 
    186 !!$         &  zmax(4) <    0._wp .OR.   &                   ! too large sea surface salinity (keep this line for sea-ice) 
    187          &  ISNAN( zmax(1) + zmax(2) + zmax(3) ) .OR.   &               ! NaN encounter in the tests 
    188          &  ABS(   zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN   ! Infinity encounter in the tests 
     149      !                                   !==  error handling  ==! 
     150      IF( ( sn_cfctl%l_glochk .OR. lsomeoce ) .AND. (   &  ! domain contains some ocean points, check for sensible ranges 
     151         &  zmax(1) >   20._wp .OR.   &                    ! too large sea surface height ( > 20 m ) 
     152         &  zmax(2) >   10._wp .OR.   &                    ! too large velocity ( > 10 m/s) 
     153!!$         &  zmax(3) >=   0._wp .OR.   &                    ! negative or zero sea surface salinity 
     154!!$         &  zmax(4) >= 100._wp .OR.   &                    ! too large sea surface salinity ( > 100 ) 
     155!!$         &  zmax(4) <    0._wp .OR.   &                    ! too large sea surface salinity (keep this line for sea-ice) 
     156         &  ISNAN( zmax(1) + zmax(2) + zmax(3) ) ) ) THEN   ! NaN encounter in the tests 
     157         IF( lk_mpp .AND. sn_cfctl%l_glochk ) THEN 
     158            ! have use mpp_max (because sn_cfctl%l_glochk=.T. and distributed) 
     159            CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:,Kmm))        , ssmask(:,:)  , zzz, ih  ) 
     160            CALL mpp_maxloc( 'stpctl', ABS(uu(:,:,:,Kmm))          , umask (:,:,:), zzz, iu  ) 
     161            CALL mpp_minloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is1 ) 
     162            CALL mpp_maxloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is2 ) 
     163         ELSE 
     164            ! find local min and max locations 
     165            ih(:)  = MAXLOC( ABS( ssh(:,:,Kmm)   )                              ) + (/ nimpp - 1, njmpp - 1    /) 
     166            iu(:)  = MAXLOC( ABS( uu  (:,:,:,Kmm) )                              ) + (/ nimpp - 1, njmpp - 1, 0 /) 
     167            is1(:) = MINLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
     168            is2(:) = MAXLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
     169         ENDIF 
     170          
     171         WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests' 
     172         WRITE(ctmp2,9100) kt,   zmax(1), ih(1) , ih(2) 
     173         WRITE(ctmp3,9200) kt,   zmax(2), iu(1) , iu(2) , iu(3) 
     174         WRITE(ctmp4,9300) kt, - zmax(3), is1(1), is1(2), is1(3) 
     175         WRITE(ctmp5,9400) kt,   zmax(4), is2(1), is2(2), is2(3) 
     176         WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort.nc file' 
     177          
     178         CALL dia_wri_state( Kmm, 'output.abort' )     ! create an output.abort file 
     179          
     180         IF( .NOT. sn_cfctl%l_glochk ) THEN 
     181            WRITE(ctmp8,*) 'E R R O R message from sub-domain: ', narea 
     182            CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp8, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ctmp6 ) 
     183         ELSE 
     184            CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6, ' ' ) 
     185         ENDIF 
     186 
     187         kindic = -3 
    189188         ! 
    190          iloc(:,:) = 0 
    191          IF( ll_colruns ) THEN   ! zmax is global, so it is the same on all subdomains -> no dead lock with mpp_maxloc 
    192             ! first: close the netcdf file, so we can read it 
    193             IF( lwm .AND. kt /= nitend )   istatus = NF90_CLOSE(nrunid) 
    194             ! get global loc on the min/max 
    195             CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:,         Kmm)), ssmask(:,:  ), zzz, iloc(1:2,1) )   ! mpp_maxloc ok if mask = F  
    196             CALL mpp_maxloc( 'stpctl', ABS( uu(:,:,:,       Kmm)),  umask(:,:,:), zzz, iloc(1:3,2) ) 
    197             CALL mpp_minloc( 'stpctl',      ts(:,:,:,jp_sal,Kmm) ,  tmask(:,:,:), zzz, iloc(1:3,3) ) 
    198             CALL mpp_maxloc( 'stpctl',      ts(:,:,:,jp_sal,Kmm) ,  tmask(:,:,:), zzz, iloc(1:3,4) ) 
    199             ! find which subdomain has the max. 
    200             iareamin(:) = jpnij+1   ;   iareamax(:) = 0   ;   iareasum(:) = 0 
    201             DO ji = 1, 9 
    202                IF( zmaxlocal(ji) == zmax(ji) ) THEN 
    203                   iareamin(ji) = narea   ;   iareamax(ji) = narea   ;   iareasum(ji) = 1 
    204                ENDIF 
    205             END DO 
    206             CALL mpp_min( "stpctl", iareamin )         ! min over the global domain 
    207             CALL mpp_max( "stpctl", iareamax )         ! max over the global domain 
    208             CALL mpp_sum( "stpctl", iareasum )         ! sum over the global domain 
    209          ELSE                    ! find local min and max locations: 
    210             ! if we are here, this means that the subdomain contains some oce points -> no need to test the mask used in maxloc 
    211             iloc(1:2,1) = MAXLOC( ABS( ssh(:,:,         Kmm)), mask = ssmask(:,:  ) == 1._wp ) + (/ nimpp - 1, njmpp - 1    /) 
    212             iloc(1:3,2) = MAXLOC( ABS(  uu(:,:,:,       Kmm)), mask =  umask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
    213             iloc(1:3,3) = MINLOC(       ts(:,:,:,jp_sal,Kmm) , mask =  tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
    214             iloc(1:3,4) = MAXLOC(       ts(:,:,:,jp_sal,Kmm) , mask =  tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
    215             iareamin(:) = narea   ;   iareamax(:) = narea   ;   iareasum(:) = 1         ! this is local information 
    216          ENDIF 
    217          ! 
    218          WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests' 
    219          CALL wrt_line( ctmp2, kt, '|ssh| max',  zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) ) 
    220          CALL wrt_line( ctmp3, kt, '|U|   max',  zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) ) 
    221          CALL wrt_line( ctmp4, kt, 'Sal   min', -zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) ) 
    222          CALL wrt_line( ctmp5, kt, 'Sal   max',  zmax(4), iloc(:,4), iareasum(4), iareamin(4), iareamax(4) ) 
    223          IF( Agrif_Root() ) THEN 
    224             WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort* files' 
    225          ELSE 
    226             WRITE(ctmp6,*) '      ===> output of last computed fields in '//TRIM(Agrif_CFixed())//'_output.abort* files' 
    227          ENDIF 
    228          ! 
    229          CALL dia_wri_state( Kmm, 'output.abort' )     ! create an output.abort file 
    230          ! 
    231          IF( ll_colruns .or. jpnij == 1 ) THEN   ! all processes synchronized -> use lwp to print in opened ocean.output files 
    232             IF(lwp) THEN   ;   CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 ) 
    233             ELSE           ;   nstop = MAX(1, nstop)   ! make sure nstop > 0 (automatically done when calling ctl_stop) 
    234             ENDIF 
    235          ELSE                                    ! only mpi subdomains with errors are here -> STOP now 
    236             CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 ) 
    237          ENDIF 
    238          ! 
    239       ENDIF 
    240       ! 
    241       IF( nstop > 0 ) THEN                                                  ! an error was detected and we did not abort yet... 
    242          ngrdstop = Agrif_Fixed()                                           ! store which grid got this error 
    243          IF( .NOT. ll_colruns .AND. jpnij > 1 )   CALL ctl_stop( 'STOP' )   ! we must abort here to avoid MPI deadlock 
    244       ENDIF 
    245       ! 
     189      ENDIF 
     190      ! 
     1919100  FORMAT (' kt=',i8,'   |ssh| max: ',1pg11.4,', at  i j  : ',2i5) 
     1929200  FORMAT (' kt=',i8,'   |U|   max: ',1pg11.4,', at  i j k: ',3i5) 
     1939300  FORMAT (' kt=',i8,'   S     min: ',1pg11.4,', at  i j k: ',3i5) 
     1949400  FORMAT (' kt=',i8,'   S     max: ',1pg11.4,', at  i j k: ',3i5) 
    2461959500  FORMAT(' it :', i8, '    |ssh|_max: ', D23.16, ' |U|_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16) 
    247196      ! 
    248197   END SUBROUTINE stp_ctl 
    249  
    250  
    251    SUBROUTINE wrt_line( cdline, kt, cdprefix, pval, kloc, ksum, kmin, kmax ) 
    252       !!---------------------------------------------------------------------- 
    253       !!                     ***  ROUTINE wrt_line  *** 
    254       !! 
    255       !! ** Purpose :   write information line 
    256       !! 
    257       !!---------------------------------------------------------------------- 
    258       CHARACTER(len=*),      INTENT(  out) ::   cdline 
    259       CHARACTER(len=*),      INTENT(in   ) ::   cdprefix 
    260       REAL(wp),              INTENT(in   ) ::   pval 
    261       INTEGER, DIMENSION(3), INTENT(in   ) ::   kloc 
    262       INTEGER,               INTENT(in   ) ::   kt, ksum, kmin, kmax 
    263       ! 
    264       CHARACTER(len=80) ::   clsuff 
    265       CHARACTER(len=9 ) ::   clkt, clsum, clmin, clmax 
    266       CHARACTER(len=9 ) ::   cli, clj, clk 
    267       CHARACTER(len=1 ) ::   clfmt 
    268       CHARACTER(len=4 ) ::   cl4   ! needed to be able to compile with Agrif, I don't know why 
    269       INTEGER           ::   ifmtk 
    270       !!---------------------------------------------------------------------- 
    271       WRITE(clkt , '(i9)') kt 
    272        
    273       WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpnij  ,wp))) + 1     ! how many digits to we need to write ? (we decide max = 9) 
    274       !!! WRITE(clsum, '(i'//clfmt//')') ksum                   ! this is creating a compilation error with AGRIF 
    275       cl4 = '(i'//clfmt//')'   ;   WRITE(clsum, cl4) ksum 
    276       WRITE(clfmt, '(i1)') INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1    ! how many digits to we need to write ? (we decide max = 9) 
    277       cl4 = '(i'//clfmt//')'   ;   WRITE(clmin, cl4) kmin-1 
    278                                    WRITE(clmax, cl4) kmax-1 
    279       ! 
    280       WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpiglo,wp))) + 1      ! how many digits to we need to write jpiglo? (we decide max = 9) 
    281       cl4 = '(i'//clfmt//')'   ;   WRITE(cli, cl4) kloc(1)      ! this is ok with AGRIF 
    282       WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpjglo,wp))) + 1      ! how many digits to we need to write jpjglo? (we decide max = 9) 
    283       cl4 = '(i'//clfmt//')'   ;   WRITE(clj, cl4) kloc(2)      ! this is ok with AGRIF 
    284       ! 
    285       IF( ksum == 1 ) THEN   ;   WRITE(clsuff,9100) TRIM(clmin) 
    286       ELSE                   ;   WRITE(clsuff,9200) TRIM(clsum), TRIM(clmin), TRIM(clmax) 
    287       ENDIF 
    288       IF(kloc(3) == 0) THEN 
    289          ifmtk = INT(LOG10(REAL(jpk,wp))) + 1                   ! how many digits to we need to write jpk? (we decide max = 9) 
    290          clk = REPEAT(' ', ifmtk)                               ! create the equivalent in blank string 
    291          WRITE(cdline,9300) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj), clk(1:ifmtk), TRIM(clsuff) 
    292       ELSE 
    293          WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpk,wp))) + 1      ! how many digits to we need to write jpk? (we decide max = 9) 
    294          !!! WRITE(clk, '(i'//clfmt//')') kloc(3)               ! this is creating a compilation error with AGRIF 
    295          cl4 = '(i'//clfmt//')'   ;   WRITE(clk, cl4) kloc(3)   ! this is ok with AGRIF 
    296          WRITE(cdline,9400) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj),    TRIM(clk), TRIM(clsuff) 
    297       ENDIF 
    298       ! 
    299 9100  FORMAT('MPI rank ', a) 
    300 9200  FORMAT('found in ', a, ' MPI tasks, spread out among ranks ', a, ' to ', a) 
    301 9300  FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j   ', a, ' ', a, ' ', a, ' ', a) 
    302 9400  FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j k ', a, ' ', a, ' ', a, ' ', a) 
    303       ! 
    304    END SUBROUTINE wrt_line 
    305  
    306198 
    307199   !!====================================================================== 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/CANAL/MY_SRC/trazdf.F90

    r12489 r13197  
    3535   PUBLIC   tra_zdf_imp   ! called by trczdf.F90 
    3636 
     37   !! * Substitutions 
     38#  include "do_loop_substitute.h90" 
    3739   !!---------------------------------------------------------------------- 
    3840   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    7779      ! JMM avoid negative salinities near river outlet ! Ugly fix 
    7880      ! JMM : restore negative salinities to small salinities: 
    79 !!$   WHERE( pts(:,:,:,jp_sal,Kaa) < 0._wp )   pts(:,:,:,jp_sal,Kaa) = 0.1_wp 
     81!!$      WHERE( pts(:,:,:,jp_sal,Kaa) < 0._wp )   pts(:,:,:,jp_sal,Kaa) = 0.1_wp 
    8082!!gm 
    8183 
     
    9597      ENDIF 
    9698      !                                          ! print mean trends (used for debugging) 
    97       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kaa), clinfo1=' zdf  - Ta: ', mask1=tmask,               & 
    98          &                       tab3d_2=pts(:,:,:,jp_sal,Kaa), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     99      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kaa), clinfo1=' zdf  - Ta: ', mask1=tmask,               & 
     100         &                                  tab3d_2=pts(:,:,:,jp_sal,Kaa), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    99101      ! 
    100102      IF( ln_timing )   CALL timing_stop('tra_zdf') 
     
    154156            IF( l_ldfslp ) THEN            ! isoneutral diffusion: add the contribution  
    155157               IF( ln_traldf_msc  ) THEN     ! MSC iso-neutral operator  
    156                   DO jk = 2, jpkm1 
    157                      DO jj = 2, jpjm1 
    158                         DO ji = fs_2, fs_jpim1   ! vector opt. 
    159                            zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)   
    160                         END DO 
    161                      END DO 
    162                   END DO 
     158                  DO_3D_00_00( 2, jpkm1 ) 
     159                     zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)   
     160                  END_3D 
    163161               ELSE                          ! standard or triad iso-neutral operator 
    164                   DO jk = 2, jpkm1 
    165                      DO jj = 2, jpjm1 
    166                         DO ji = fs_2, fs_jpim1   ! vector opt. 
    167                            zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk) 
    168                         END DO 
    169                      END DO 
    170                   END DO 
     162                  DO_3D_00_00( 2, jpkm1 ) 
     163                     zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk) 
     164                  END_3D 
    171165               ENDIF 
    172166            ENDIF 
     
    174168            ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked) 
    175169            IF( ln_zad_Aimp ) THEN         ! Adaptive implicit vertical advection 
    176                DO jk = 1, jpkm1 
    177                   DO jj = 2, jpjm1 
    178                      DO ji = fs_2, fs_jpim1   ! vector opt. (ensure same order of calculation as below if wi=0.) 
    179                         zzwi = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk  ,Kmm) 
    180                         zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 
    181                         zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws   & 
    182                            &                 + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) 
    183                         zwi(ji,jj,jk) = zzwi + p2dt *   MIN( wi(ji,jj,jk  ) , 0._wp ) 
    184                         zws(ji,jj,jk) = zzws - p2dt *   MAX( wi(ji,jj,jk+1) , 0._wp ) 
    185                     END DO 
    186                   END DO 
    187                END DO 
     170               DO_3D_00_00( 1, jpkm1 ) 
     171                  zzwi = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk  ,Kmm) 
     172                  zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 
     173                  zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws   & 
     174                     &                 + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) 
     175                  zwi(ji,jj,jk) = zzwi + p2dt *   MIN( wi(ji,jj,jk  ) , 0._wp ) 
     176                  zws(ji,jj,jk) = zzws - p2dt *   MAX( wi(ji,jj,jk+1) , 0._wp ) 
     177               END_3D 
    188178            ELSE 
    189                DO jk = 1, jpkm1 
    190                   DO jj = 2, jpjm1 
    191                      DO ji = fs_2, fs_jpim1   ! vector opt. 
    192                         zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk,Kmm) 
    193                         zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 
    194                         zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk) 
    195                     END DO 
    196                   END DO 
    197                END DO 
     179               DO_3D_00_00( 1, jpkm1 ) 
     180                  zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk,Kmm) 
     181                  zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 
     182                  zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk) 
     183               END_3D 
    198184            ENDIF 
    199185            ! 
     
    217203            !   used as a work space array: its value is modified. 
    218204            ! 
    219             DO jj = 2, jpjm1        !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k) 
    220                DO ji = fs_2, fs_jpim1            ! done one for all passive tracers (so included in the IF instruction) 
    221                   zwt(ji,jj,1) = zwd(ji,jj,1) 
    222                END DO 
    223             END DO 
    224             DO jk = 2, jpkm1 
    225                DO jj = 2, jpjm1 
    226                   DO ji = fs_2, fs_jpim1 
    227                      zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 
    228                   END DO 
    229                END DO 
    230             END DO 
     205            DO_2D_00_00 
     206               zwt(ji,jj,1) = zwd(ji,jj,1) 
     207            END_2D 
     208            DO_3D_00_00( 2, jpkm1 ) 
     209               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 
     210            END_3D 
    231211            ! 
    232212         ENDIF  
    233213         !          
    234          DO jj = 2, jpjm1           !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
    235             DO ji = fs_2, fs_jpim1 
    236                pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs) 
    237             END DO 
    238          END DO 
    239          DO jk = 2, jpkm1 
    240             DO jj = 2, jpjm1 
    241                DO ji = fs_2, fs_jpim1 
    242                   zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs)   ! zrhs=right hand side 
    243                   pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa) 
    244                END DO 
    245             END DO 
    246          END DO 
     214         DO_2D_00_00 
     215            pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs) 
     216         END_2D 
     217         DO_3D_00_00( 2, jpkm1 ) 
     218            zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs)   ! zrhs=right hand side 
     219            pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa) 
     220         END_3D 
    247221         ! 
    248          DO jj = 2, jpjm1           !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer) 
    249             DO ji = fs_2, fs_jpim1 
    250                pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 
    251             END DO 
    252          END DO 
    253          DO jk = jpk-2, 1, -1 
    254             DO jj = 2, jpjm1 
    255                DO ji = fs_2, fs_jpim1 
    256                   pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) )   & 
    257                      &             / zwt(ji,jj,jk) * tmask(ji,jj,jk) 
    258                END DO 
    259             END DO 
    260          END DO 
     222         DO_2D_00_00 
     223            pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 
     224         END_2D 
     225         DO_3DS_00_00( jpk-2, 1, -1 ) 
     226            pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) )   & 
     227               &             / zwt(ji,jj,jk) * tmask(ji,jj,jk) 
     228         END_3D 
    261229         !                                            ! ================= ! 
    262230      END DO                                          !  end tracer loop  ! 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/CANAL/MY_SRC/usrdef_hgr.F90

    r10074 r13197  
    2626   PUBLIC   usr_def_hgr   ! called by domhgr.F90 
    2727 
     28   !! * Substitutions 
     29#  include "do_loop_substitute.h90" 
    2830   !!---------------------------------------------------------------------- 
    2931   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    8890#endif 
    8991          
    90       DO jj = 1, jpj 
    91          DO ji = 1, jpi 
    92             zti = FLOAT( ji - 1 + nimpp - 1 )          ;  ztj = FLOAT( jj - 1 + njmpp - 1 ) 
    93             zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5_wp ;  zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5_wp 
    94  
    95             plamt(ji,jj) = zlam0 + rn_dx * zti 
    96             plamu(ji,jj) = zlam0 + rn_dx * zui 
    97             plamv(ji,jj) = plamt(ji,jj)  
    98             plamf(ji,jj) = plamu(ji,jj)  
    99     
    100             pphit(ji,jj) = zphi0 + rn_dy * ztj 
    101             pphiv(ji,jj) = zphi0 + rn_dy * zvj 
    102             pphiu(ji,jj) = pphit(ji,jj)  
    103             pphif(ji,jj) = pphiv(ji,jj)  
    104          END DO 
    105       END DO 
     92      DO_2D_11_11 
     93         zti = FLOAT( ji - 1 + nimpp - 1 )          ;  ztj = FLOAT( jj - 1 + njmpp - 1 ) 
     94         zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5_wp ;  zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5_wp 
     95          
     96         plamt(ji,jj) = zlam0 + rn_dx * zti 
     97         plamu(ji,jj) = zlam0 + rn_dx * zui 
     98         plamv(ji,jj) = plamt(ji,jj)  
     99         plamf(ji,jj) = plamu(ji,jj)  
     100          
     101         pphit(ji,jj) = zphi0 + rn_dy * ztj 
     102         pphiv(ji,jj) = zphi0 + rn_dy * zvj 
     103         pphiu(ji,jj) = pphit(ji,jj)  
     104         pphif(ji,jj) = pphiv(ji,jj)  
     105      END_2D 
    106106      !      
    107107      ! Horizontal scale factors (in meters) 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/CANAL/MY_SRC/usrdef_istate.F90

    r12489 r13197  
    2828   PUBLIC   usr_def_istate   ! called by istate.F90 
    2929 
     30   !! * Substitutions 
     31#  include "do_loop_substitute.h90" 
    3032   !!---------------------------------------------------------------------- 
    3133   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    164166         pssh(:,1) = - ff_t(:,1) / grav * pu(:,1,1) * e2t(:,1) 
    165167         DO jl=1, jpnj 
    166             DO jj=nldj, nlej 
    167                DO ji=nldi, nlei 
    168                   pssh(ji,jj) = pssh(ji,jj-1) - ff_t(ji,jj) / grav * pu(ji,jj,1) * e2t(ji,jj) 
    169                END DO 
    170             END DO 
     168            DO_2D_00_00 
     169               pssh(ji,jj) = pssh(ji,jj-1) - ff_t(ji,jj) / grav * pu(ji,jj,1) * e2t(ji,jj) 
     170            END_2D 
    171171            CALL lbc_lnk( 'usrdef_istate', pssh, 'T',  1. ) 
    172172         END DO 
     
    183183      CASE(4)    ! geostrophic zonal pulse 
    184184    
    185          DO jj=1, jpj 
    186             DO ji=1, jpi 
    187                IF ( ABS(glamt(ji,jj)) <= zjetx ) THEN 
    188                   zdu = rn_uzonal 
    189                ELSEIF ( ABS(glamt(ji,jj)) <= zjetx + 100. ) THEN 
    190                   zdu = rn_uzonal * ( ( zjetx-ABS(glamt(ji,jj)) )/100. + 1. ) 
    191                ELSE 
    192                   zdu = 0. 
    193                END IF 
    194                IF ( ABS(gphit(ji,jj)) <= zjety ) THEN 
    195                   pssh(ji,jj) = - ff_t(ji,jj) * zdu * gphit(ji,jj) * 1.e3 / grav 
    196                   pu(ji,jj,:) = zdu 
    197                   pts(ji,jj,:,jp_sal) = zdu / rn_uzonal + 1. 
    198                ELSE 
    199                   pssh(ji,jj) = - ff_t(ji,jj) * zdu * SIGN(zjety,gphit(ji,jj)) * 1.e3 / grav  
    200                   pu(ji,jj,:) = 0. 
    201                   pts(ji,jj,:,jp_sal) = 1. 
    202                END IF 
    203             END DO 
    204          END DO 
     185         DO_2D_11_11 
     186            IF ( ABS(glamt(ji,jj)) <= zjetx ) THEN 
     187               zdu = rn_uzonal 
     188            ELSEIF ( ABS(glamt(ji,jj)) <= zjetx + 100. ) THEN 
     189               zdu = rn_uzonal * ( ( zjetx-ABS(glamt(ji,jj)) )/100. + 1. ) 
     190            ELSE 
     191               zdu = 0. 
     192            END IF 
     193            IF ( ABS(gphit(ji,jj)) <= zjety ) THEN 
     194               pssh(ji,jj) = - ff_t(ji,jj) * zdu * gphit(ji,jj) * 1.e3 / grav 
     195               pu(ji,jj,:) = zdu 
     196               pts(ji,jj,:,jp_sal) = zdu / rn_uzonal + 1. 
     197            ELSE 
     198               pssh(ji,jj) = - ff_t(ji,jj) * zdu * SIGN(zjety,gphit(ji,jj)) * 1.e3 / grav  
     199               pu(ji,jj,:) = 0. 
     200               pts(ji,jj,:,jp_sal) = 1. 
     201            END IF 
     202         END_2D 
    205203          
    206204         ! temperature: 
    207205         pts(:,:,:,jp_tem) = 10._wp * ptmask(:,:,:)         
    208206         pv(:,:,:) = 0. 
    209           
    210207          
    211208       CASE(5)    ! vortex 
     
    220217         zP0 = rho0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp) 
    221218         ! 
    222          DO jj=1, jpj 
    223             DO ji=1, jpi 
    224                zx = glamt(ji,jj) * 1.e3 
    225                zy = gphit(ji,jj) * 1.e3 
    226                ! Surface pressure: P(x,y,z) = F(z) * Psurf(x,y) 
    227                zpsurf = zP0 * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal * zy 
    228                ! Sea level: 
    229                pssh(ji,jj) = 0. 
    230                DO jl=1,5 
    231                   zdt = pssh(ji,jj) 
    232                   zdzF = (1._wp - EXP(zdt-zH)) / (zH - 1._wp + EXP(-zH))   ! F'(z) 
    233                   zrho1 = rho0 * (1._wp + zn2*zdt/grav) - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 
    234                   pssh(ji,jj) = zpsurf / (zrho1*grav) * ptmask(ji,jj,1)   ! ssh = Psurf / (Rho*g) 
    235                END DO 
    236                ! temperature: 
    237                DO jk=1,jpk 
    238                   zdt =  pdept(ji,jj,jk)  
    239                   zrho1 = rho0 * (1._wp + zn2*zdt/grav) 
    240                   IF (zdt < zH) THEN 
    241                      zdzF = (1._wp-EXP(zdt-zH)) / (zH-1._wp + EXP(-zH))   ! F'(z) 
    242                      zrho1 = zrho1 - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 
    243                   ENDIF 
    244                   !               pts(ji,jj,jk,jp_tem) = (20._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 
    245                   pts(ji,jj,jk,jp_tem) = (10._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 
    246                END DO 
    247             END DO 
    248          END DO 
     219         DO_2D_11_11 
     220            zx = glamt(ji,jj) * 1.e3 
     221            zy = gphit(ji,jj) * 1.e3 
     222            ! Surface pressure: P(x,y,z) = F(z) * Psurf(x,y) 
     223            zpsurf = zP0 * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal * zy 
     224            ! Sea level: 
     225            pssh(ji,jj) = 0. 
     226            DO jl=1,5 
     227               zdt = pssh(ji,jj) 
     228               zdzF = (1._wp - EXP(zdt-zH)) / (zH - 1._wp + EXP(-zH))   ! F'(z) 
     229               zrho1 = rho0 * (1._wp + zn2*zdt/grav) - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 
     230               pssh(ji,jj) = zpsurf / (zrho1*grav) * ptmask(ji,jj,1)   ! ssh = Psurf / (Rho*g) 
     231            END DO 
     232            ! temperature: 
     233            DO jk=1,jpk 
     234               zdt =  pdept(ji,jj,jk)  
     235               zrho1 = rho0 * (1._wp + zn2*zdt/grav) 
     236               IF (zdt < zH) THEN 
     237                  zdzF = (1._wp-EXP(zdt-zH)) / (zH-1._wp + EXP(-zH))   ! F'(z) 
     238                  zrho1 = zrho1 - zdzF * zpsurf / grav    ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 
     239               ENDIF 
     240               !               pts(ji,jj,jk,jp_tem) = (20._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 
     241               pts(ji,jj,jk,jp_tem) = (10._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 
     242            END DO 
     243         END_2D 
    249244         ! 
    250245         ! salinity:   
     
    253248         ! velocities: 
    254249         za = 2._wp * zP0 / zlambda**2 
    255          DO jj=1, jpj 
    256             DO ji=1, jpim1 
    257                zx = glamu(ji,jj) * 1.e3 
    258                zy = gphiu(ji,jj) * 1.e3 
    259                DO jk=1, jpk 
    260                   zdu = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji+1,jj,jk)) 
    261                   IF (zdu < zH) THEN 
    262                      zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH)) 
    263                      zdyPs = - za * zy * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal 
    264                      pu(ji,jj,jk) = - zf / ( rho0 * ff_t(ji,jj) ) * zdyPs * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk) 
    265                   ELSE 
    266                      pu(ji,jj,jk) = 0._wp 
    267                   ENDIF 
    268                END DO 
    269             END DO 
    270          END DO 
    271          ! 
    272          DO jj=1, jpjm1 
    273             DO ji=1, jpi 
    274                zx = glamv(ji,jj) * 1.e3 
    275                zy = gphiv(ji,jj) * 1.e3 
    276                DO jk=1, jpk 
    277                   zdv = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji,jj+1,jk)) 
    278                   IF (zdv < zH) THEN 
    279                      zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH)) 
    280                      zdxPs = - za * zx * EXP(-(zx**2+zy**2)*zr_lambda2) 
    281                      pv(ji,jj,jk) = zf / ( rho0 * ff_f(ji,jj) ) * zdxPs * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk) 
    282                   ELSE 
    283                      pv(ji,jj,jk) = 0._wp 
    284                   ENDIF 
    285                END DO 
    286             END DO 
    287          END DO 
     250         DO_2D_00_00 
     251            zx = glamu(ji,jj) * 1.e3 
     252            zy = gphiu(ji,jj) * 1.e3 
     253            DO jk=1, jpk 
     254               zdu = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji+1,jj,jk)) 
     255               IF (zdu < zH) THEN 
     256                  zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH)) 
     257                  zdyPs = - za * zy * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal 
     258                  pu(ji,jj,jk) = - zf / ( rho0 * ff_t(ji,jj) ) * zdyPs * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk) 
     259               ELSE 
     260                  pu(ji,jj,jk) = 0._wp 
     261               ENDIF 
     262            END DO 
     263         END_2D 
     264         ! 
     265         DO_2D_00_00 
     266            zx = glamv(ji,jj) * 1.e3 
     267            zy = gphiv(ji,jj) * 1.e3 
     268            DO jk=1, jpk 
     269               zdv = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji,jj+1,jk)) 
     270               IF (zdv < zH) THEN 
     271                  zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH)) 
     272                  zdxPs = - za * zx * EXP(-(zx**2+zy**2)*zr_lambda2) 
     273                  pv(ji,jj,jk) = zf / ( rho0 * ff_f(ji,jj) ) * zdxPs * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk) 
     274               ELSE 
     275                  pv(ji,jj,jk) = 0._wp 
     276               ENDIF 
     277            END DO 
     278         END_2D 
    288279         !             
    289280      END SELECT 
    290  
     281       
    291282      IF (ln_sshnoise) THEN 
    292283         CALL RANDOM_NUMBER(zrandom) 
     
    294285      END IF 
    295286      CALL lbc_lnk( 'usrdef_istate', pssh, 'T',  1. ) 
    296       CALL lbc_lnk(  'usrdef_istate', pts, 'T',  1. ) 
    297       CALL lbc_lnk(   'usrdef_istate', pu, 'U', -1. ) 
    298       CALL lbc_lnk(   'usrdef_istate', pv, 'V', -1. ) 
     287      CALL lbc_lnk( 'usrdef_istate', pts , 'T',  1. ) 
     288      CALL lbc_lnk_multi( 'usrdef_istate', pu, 'U', -1., pv, 'V', -1. ) 
    299289 
    300290   END SUBROUTINE usr_def_istate 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/CANAL/MY_SRC/usrdef_sbc.F90

    r12377 r13197  
    3838CONTAINS 
    3939 
    40    SUBROUTINE usrdef_sbc_oce( kt, Kmm, Kbb ) 
     40   SUBROUTINE usrdef_sbc_oce( kt, Kbb ) 
    4141      !!--------------------------------------------------------------------- 
    4242      !!                    ***  ROUTINE usr_def_sbc  *** 
     
    5353      !!---------------------------------------------------------------------- 
    5454      INTEGER, INTENT(in) ::   kt        ! ocean time step 
    55       INTEGER, INTENT(in) ::   Kbb, Kmm  ! ocean time index 
     55      INTEGER, INTENT(in) ::   Kbb       ! ocean time index 
    5656      INTEGER  ::   ji, jj               ! dummy loop indices 
    5757      REAL(wp) :: zrhoair = 1.22     ! approximate air density [Kg/m3] 
     
    8686          
    8787         WHERE( ABS(gphit) <= rn_windszy/2. ) 
    88             zwndrel(:,:) = rn_u10 - rn_uofac * uu(:,:,1,Kmm) 
     88            zwndrel(:,:) = rn_u10 - rn_uofac * uu(:,:,1,Kbb) 
    8989         ELSEWHERE 
    90             zwndrel(:,:) =        - rn_uofac * uu(:,:,1,Kmm) 
     90            zwndrel(:,:) =        - rn_uofac * uu(:,:,1,Kbb) 
    9191         END WHERE 
    9292         utau(:,:) = zrhocd * zwndrel(:,:) * zwndrel(:,:) 
    9393 
    94          zwndrel(:,:) = - rn_uofac * vv(:,:,1,Kmm) 
     94         zwndrel(:,:) = - rn_uofac * vv(:,:,1,Kbb) 
    9595         vtau(:,:) = zrhocd * zwndrel(:,:) * zwndrel(:,:) 
    9696 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/CANAL/MY_SRC/usrdef_zgr.F90

    r12377 r13197  
    204204      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. )           ! set surrounding land to zero (here jperio=0 ==>> closed) 
    205205      ! 
    206       k_bot(:,:) = INT( z2d(:,:) )           ! =jpkm1 over the ocean point, =0 elsewhere 
     206      k_bot(:,:) = NINT( z2d(:,:) )          ! =jpkm1 over the ocean point, =0 elsewhere 
    207207      ! 
    208208      k_top(:,:) = MIN( 1 , k_bot(:,:) )     ! = 1    over the ocean point, =0 elsewhere 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/ICE_ADV1D/MY_SRC/usrdef_hgr.F90

    r10513 r13197  
    2626   PUBLIC   usr_def_hgr   ! called by domhgr.F90 
    2727 
     28   !! * Substitutions 
     29#  include "do_loop_substitute.h90" 
    2830   !!---------------------------------------------------------------------- 
    2931   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    7678      zphi0 = -(jpjglo-1)/2 * 1.e-3 * rn_dy 
    7779 
    78       DO jj = 1, jpj 
    79          DO ji = 1, jpi 
    80             zti = FLOAT( ji - 1 + nimpp - 1 )          ;  ztj = FLOAT( jj - 1 + njmpp - 1 ) 
    81             zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5_wp ;  zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5_wp 
    82  
    83             plamt(ji,jj) = zlam0 + rn_dx * 1.e-3 * zti 
    84             plamu(ji,jj) = zlam0 + rn_dx * 1.e-3 * zui 
    85             plamv(ji,jj) = plamt(ji,jj)  
    86             plamf(ji,jj) = plamu(ji,jj)  
    87     
    88             pphit(ji,jj) = zphi0 + rn_dy * 1.e-3 * ztj 
    89             pphiv(ji,jj) = zphi0 + rn_dy * 1.e-3 * zvj 
    90             pphiu(ji,jj) = pphit(ji,jj)  
    91             pphif(ji,jj) = pphiv(ji,jj)  
    92          END DO 
    93       END DO 
     80      DO_2D_11_11 
     81         zti = FLOAT( ji - 1 + nimpp - 1 )          ;  ztj = FLOAT( jj - 1 + njmpp - 1 ) 
     82         zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5_wp ;  zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5_wp 
     83          
     84         plamt(ji,jj) = zlam0 + rn_dx * 1.e-3 * zti 
     85         plamu(ji,jj) = zlam0 + rn_dx * 1.e-3 * zui 
     86         plamv(ji,jj) = plamt(ji,jj)  
     87         plamf(ji,jj) = plamu(ji,jj)  
     88          
     89         pphit(ji,jj) = zphi0 + rn_dy * 1.e-3 * ztj 
     90         pphiv(ji,jj) = zphi0 + rn_dy * 1.e-3 * zvj 
     91         pphiu(ji,jj) = pphit(ji,jj)  
     92         pphif(ji,jj) = pphiv(ji,jj)  
     93      END_2D 
    9494          
    9595      ! constant scale factors 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/ICE_ADV2D/MY_SRC/usrdef_hgr.F90

    r10515 r13197  
    2626   PUBLIC   usr_def_hgr   ! called by domhgr.F90 
    2727 
     28   !! * Substitutions 
     29#  include "do_loop_substitute.h90" 
    2830   !!---------------------------------------------------------------------- 
    2931   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    8890#endif          
    8991 
    90       DO jj = 1, jpj 
    91          DO ji = 1, jpi 
    92             zti = FLOAT( ji - 1 + nimpp - 1 )          ;  ztj = FLOAT( jj - 1 + njmpp - 1 ) 
    93             zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5_wp ;  zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5_wp 
    94  
    95             plamt(ji,jj) = zlam0 + rn_dx * 1.e-3 * zti 
    96             plamu(ji,jj) = zlam0 + rn_dx * 1.e-3 * zui 
    97             plamv(ji,jj) = plamt(ji,jj)  
    98             plamf(ji,jj) = plamu(ji,jj)  
    99     
    100             pphit(ji,jj) = zphi0 + rn_dy * 1.e-3 * ztj 
    101             pphiv(ji,jj) = zphi0 + rn_dy * 1.e-3 * zvj 
    102             pphiu(ji,jj) = pphit(ji,jj)  
    103             pphif(ji,jj) = pphiv(ji,jj)  
    104          END DO 
    105       END DO 
     92      DO_2D_11_11 
     93         zti = FLOAT( ji - 1 + nimpp - 1 )          ;  ztj = FLOAT( jj - 1 + njmpp - 1 ) 
     94         zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5_wp ;  zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5_wp 
     95          
     96         plamt(ji,jj) = zlam0 + rn_dx * 1.e-3 * zti 
     97         plamu(ji,jj) = zlam0 + rn_dx * 1.e-3 * zui 
     98         plamv(ji,jj) = plamt(ji,jj)  
     99         plamf(ji,jj) = plamu(ji,jj)  
     100          
     101         pphit(ji,jj) = zphi0 + rn_dy * 1.e-3 * ztj 
     102         pphiv(ji,jj) = zphi0 + rn_dy * 1.e-3 * zvj 
     103         pphiu(ji,jj) = pphit(ji,jj)  
     104         pphif(ji,jj) = pphiv(ji,jj)  
     105      END_2D 
    106106          
    107107         ! Horizontal scale factors (in meters) 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/ICE_ADV2D/MY_SRC/usrdef_nam.F90

    r12377 r13197  
    1414   !!   usr_def_hgr   : initialize the horizontal mesh  
    1515   !!---------------------------------------------------------------------- 
    16    USE dom_oce  , ONLY: nimpp , njmpp            ! i- & j-indices of the local domain 
     16   USE dom_oce  , ONLY: nimpp , njmpp, Agrif_Root            ! i- & j-indices of the local domain 
    1717   USE par_oce        ! ocean space and time domain 
    1818   USE phycst         ! physical constants 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/ICE_AGRIF/MY_SRC/usrdef_hgr.F90

    r10516 r13197  
    2626   PUBLIC   usr_def_hgr   ! called by domhgr.F90 
    2727 
     28   !! * Substitutions 
     29#  include "do_loop_substitute.h90" 
    2830   !!---------------------------------------------------------------------- 
    2931   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    8890#endif          
    8991 
    90       DO jj = 1, jpj 
    91          DO ji = 1, jpi 
    92             zti = FLOAT( ji - 1 + nimpp - 1 )          ;  ztj = FLOAT( jj - 1 + njmpp - 1 ) 
    93             zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5_wp ;  zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5_wp 
    94  
    95             plamt(ji,jj) = zlam0 + rn_dx * 1.e-3 * zti 
    96             plamu(ji,jj) = zlam0 + rn_dx * 1.e-3 * zui 
    97             plamv(ji,jj) = plamt(ji,jj)  
    98             plamf(ji,jj) = plamu(ji,jj)  
    99     
    100             pphit(ji,jj) = zphi0 + rn_dy * 1.e-3 * ztj 
    101             pphiv(ji,jj) = zphi0 + rn_dy * 1.e-3 * zvj 
    102             pphiu(ji,jj) = pphit(ji,jj)  
    103             pphif(ji,jj) = pphiv(ji,jj)  
    104          END DO 
    105       END DO 
     92      DO_2D_11_11 
     93         zti = FLOAT( ji - 1 + nimpp - 1 )          ;  ztj = FLOAT( jj - 1 + njmpp - 1 ) 
     94         zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5_wp ;  zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5_wp 
     95          
     96         plamt(ji,jj) = zlam0 + rn_dx * 1.e-3 * zti 
     97         plamu(ji,jj) = zlam0 + rn_dx * 1.e-3 * zui 
     98         plamv(ji,jj) = plamt(ji,jj)  
     99         plamf(ji,jj) = plamu(ji,jj)  
     100          
     101         pphit(ji,jj) = zphi0 + rn_dy * 1.e-3 * ztj 
     102         pphiv(ji,jj) = zphi0 + rn_dy * 1.e-3 * zvj 
     103         pphiu(ji,jj) = pphit(ji,jj)  
     104         pphif(ji,jj) = pphiv(ji,jj)  
     105      END_2D 
    106106          
    107107         ! Horizontal scale factors (in meters) 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/ICE_AGRIF/MY_SRC/usrdef_nam.F90

    r12377 r13197  
    8989         kpj = nbcellsy + 2 + 2*nbghostcells 
    9090      ENDIF 
    91       kpk = 1 
     91      kpk = 2 
    9292      ! 
    9393!!      zlx = (kpi-2)*rn_dx*1.e-3 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/ICE_AGRIF/MY_SRC/usrdef_zgr.F90

    r12377 r13197  
    8989      !                       !==  z-coordinate  ==!   (step-like topography) 
    9090      !                                !* bottom ocean compute from the depth of grid-points 
    91       jpkm1 = jpk 
     91      jpkm1 = jpk-1 
    9292      k_bot(:,:) = 1    ! here use k_top as a land mask 
    9393      !                                !* horizontally uniform coordinate (reference z-co everywhere) 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/ISOMIP/MY_SRC/usrdef_hgr.F90

    r10074 r13197  
    2727   PUBLIC   usr_def_hgr   ! called by domhgr.F90 
    2828 
     29   !! * Substitutions 
     30#  include "do_loop_substitute.h90" 
    2931   !!---------------------------------------------------------------------- 
    3032   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    7577      ! 
    7678      !                       !==  grid point position  ==!   (in degrees) 
    77       DO jj = 1, jpj 
    78          DO ji = 1, jpi             ! longitude   (west coast at lon=0°) 
    79             plamt(ji,jj) = rn_e1deg * (  - 0.5 + REAL( ji-1 + nimpp-1 , wp )  )   
    80             plamu(ji,jj) = rn_e1deg * (          REAL( ji-1 + nimpp-1 , wp )  ) 
    81             plamv(ji,jj) = plamt(ji,jj) 
    82             plamf(ji,jj) = plamu(ji,jj) 
    83             !                       ! latitude   (south coast at lat= 81°) 
    84             pphit(ji,jj) = rn_e2deg * (  - 0.5 + REAL( jj-1 + njmpp-1 , wp )  ) - 80._wp 
    85             pphiu(ji,jj) = pphit(ji,jj) 
    86             pphiv(ji,jj) = rn_e2deg * (          REAL( jj-1 + njmpp-1 , wp )  ) - 80_wp 
    87             pphif(ji,jj) = pphiv(ji,jj) 
    88          END DO 
    89       END DO 
     79      DO_2D_11_11 
     80         !                       ! longitude   (west coast at lon=0°) 
     81         plamt(ji,jj) = rn_e1deg * (  - 0.5 + REAL( ji-1 + nimpp-1 , wp )  )   
     82         plamu(ji,jj) = rn_e1deg * (          REAL( ji-1 + nimpp-1 , wp )  ) 
     83         plamv(ji,jj) = plamt(ji,jj) 
     84         plamf(ji,jj) = plamu(ji,jj) 
     85         !                       ! latitude   (south coast at lat= 81°) 
     86         pphit(ji,jj) = rn_e2deg * (  - 0.5 + REAL( jj-1 + njmpp-1 , wp )  ) - 80._wp 
     87         pphiu(ji,jj) = pphit(ji,jj) 
     88         pphiv(ji,jj) = rn_e2deg * (          REAL( jj-1 + njmpp-1 , wp )  ) - 80_wp 
     89         pphif(ji,jj) = pphiv(ji,jj) 
     90      END_2D 
    9091      ! 
    9192      !                       !==  Horizontal scale factors  ==!   (in meters) 
    92       DO jj = 1, jpj 
    93          DO ji = 1, jpi 
    94             !                       ! e1   (zonal) 
    95             pe1t(ji,jj) = ra * rad * COS( rad * pphit(ji,jj) ) * rn_e1deg 
    96             pe1u(ji,jj) = ra * rad * COS( rad * pphiu(ji,jj) ) * rn_e1deg 
    97             pe1v(ji,jj) = ra * rad * COS( rad * pphiv(ji,jj) ) * rn_e1deg 
    98             pe1f(ji,jj) = ra * rad * COS( rad * pphif(ji,jj) ) * rn_e1deg 
    99             !                       ! e2   (meridional) 
    100             pe2t(ji,jj) = ra * rad * rn_e2deg 
    101             pe2u(ji,jj) = ra * rad * rn_e2deg 
    102             pe2v(ji,jj) = ra * rad * rn_e2deg 
    103             pe2f(ji,jj) = ra * rad * rn_e2deg 
    104          END DO 
    105       END DO 
     93      DO_2D_11_11 
     94         !                       ! e1   (zonal) 
     95         pe1t(ji,jj) = ra * rad * COS( rad * pphit(ji,jj) ) * rn_e1deg 
     96         pe1u(ji,jj) = ra * rad * COS( rad * pphiu(ji,jj) ) * rn_e1deg 
     97         pe1v(ji,jj) = ra * rad * COS( rad * pphiv(ji,jj) ) * rn_e1deg 
     98         pe1f(ji,jj) = ra * rad * COS( rad * pphif(ji,jj) ) * rn_e1deg 
     99         !                       ! e2   (meridional) 
     100         pe2t(ji,jj) = ra * rad * rn_e2deg 
     101         pe2u(ji,jj) = ra * rad * rn_e2deg 
     102         pe2v(ji,jj) = ra * rad * rn_e2deg 
     103         pe2f(ji,jj) = ra * rad * rn_e2deg 
     104      END_2D 
    106105      !                             ! NO reduction of grid size in some straits  
    107106      ke1e2u_v    = 0               !    ==>> u_ & v_surfaces will be computed in dom_ghr routine 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/ISOMIP/MY_SRC/usrdef_zgr.F90

    r12377 r13197  
    3030   PUBLIC   usr_def_zgr   ! called by domzgr.F90 
    3131 
     32   !! * Substitutions 
     33#  include "do_loop_substitute.h90" 
    3234   !!---------------------------------------------------------------------- 
    3335   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    132134            pe3vw(:,:,jk) = pe3w_1d (jk) 
    133135         END DO 
    134          DO jj = 1, jpj                      ! top scale factors and depth at T- and W-points 
    135             DO ji = 1, jpi 
    136                ik = k_top(ji,jj) 
    137                IF ( ik > 2 ) THEN 
    138                   ! pdeptw at the interface 
    139                   pdepw(ji,jj,ik  ) = MAX( zhisf(ji,jj) , pdepw(ji,jj,ik) ) 
    140                   ! e3t in both side of the interface 
    141                   pe3t (ji,jj,ik  ) = pdepw(ji,jj,ik+1) - pdepw(ji,jj,ik) 
    142                   ! pdept in both side of the interface (from previous e3t) 
    143                   pdept(ji,jj,ik  ) = pdepw(ji,jj,ik  ) + pe3t (ji,jj,ik  ) * 0.5_wp 
    144                   pdept(ji,jj,ik-1) = pdepw(ji,jj,ik  ) - pe3t (ji,jj,ik  ) * 0.5_wp 
    145                   ! pe3w on both side of the interface 
    146                   pe3w (ji,jj,ik+1) = pdept(ji,jj,ik+1) - pdept(ji,jj,ik  ) 
    147                   pe3w (ji,jj,ik  ) = pdept(ji,jj,ik  ) - pdept(ji,jj,ik-1) 
    148                   ! e3t into the ice shelf 
    149                   pe3t (ji,jj,ik-1) = pdepw(ji,jj,ik  ) - pdepw(ji,jj,ik-1) 
    150                   pe3w (ji,jj,ik-1) = pdept(ji,jj,ik-1) - pdept(ji,jj,ik-2) 
    151                END IF 
    152             END DO 
    153          END DO          
    154          DO jj = 1, jpj                      ! bottom scale factors and depth at T- and W-points 
    155             DO ji = 1, jpi 
    156                ik = k_bot(ji,jj) 
    157                pdepw(ji,jj,ik+1) = MIN( zht(ji,jj) , pdepw_1d(ik+1) ) 
     136         ! top scale factors and depth at T- and W-points 
     137         DO_2D_11_11 
     138            ik = k_top(ji,jj) 
     139            IF ( ik > 2 ) THEN 
     140               ! pdeptw at the interface 
     141               pdepw(ji,jj,ik  ) = MAX( zhisf(ji,jj) , pdepw(ji,jj,ik) ) 
     142               ! e3t in both side of the interface 
    158143               pe3t (ji,jj,ik  ) = pdepw(ji,jj,ik+1) - pdepw(ji,jj,ik) 
    159                pe3t (ji,jj,ik+1) = pe3t (ji,jj,ik  )  
    160                ! 
     144               ! pdept in both side of the interface (from previous e3t) 
    161145               pdept(ji,jj,ik  ) = pdepw(ji,jj,ik  ) + pe3t (ji,jj,ik  ) * 0.5_wp 
    162                pdept(ji,jj,ik+1) = pdepw(ji,jj,ik+1) + pe3t (ji,jj,ik+1) * 0.5_wp 
    163                pe3w (ji,jj,ik+1) = pdept(ji,jj,ik+1) - pdept(ji,jj,ik) 
    164             END DO 
    165          END DO          
     146               pdept(ji,jj,ik-1) = pdepw(ji,jj,ik  ) - pe3t (ji,jj,ik  ) * 0.5_wp 
     147               ! pe3w on both side of the interface 
     148               pe3w (ji,jj,ik+1) = pdept(ji,jj,ik+1) - pdept(ji,jj,ik  ) 
     149               pe3w (ji,jj,ik  ) = pdept(ji,jj,ik  ) - pdept(ji,jj,ik-1) 
     150               ! e3t into the ice shelf 
     151               pe3t (ji,jj,ik-1) = pdepw(ji,jj,ik  ) - pdepw(ji,jj,ik-1) 
     152               pe3w (ji,jj,ik-1) = pdept(ji,jj,ik-1) - pdept(ji,jj,ik-2) 
     153            END IF 
     154         END_2D 
     155         ! bottom scale factors and depth at T- and W-points 
     156         DO_2D_11_11 
     157            ik = k_bot(ji,jj) 
     158            pdepw(ji,jj,ik+1) = MIN( zht(ji,jj) , pdepw_1d(ik+1) ) 
     159            pe3t (ji,jj,ik  ) = pdepw(ji,jj,ik+1) - pdepw(ji,jj,ik) 
     160            pe3t (ji,jj,ik+1) = pe3t (ji,jj,ik  )  
     161            ! 
     162            pdept(ji,jj,ik  ) = pdepw(ji,jj,ik  ) + pe3t (ji,jj,ik  ) * 0.5_wp 
     163            pdept(ji,jj,ik+1) = pdepw(ji,jj,ik+1) + pe3t (ji,jj,ik+1) * 0.5_wp 
     164            pe3w (ji,jj,ik+1) = pdept(ji,jj,ik+1) - pdept(ji,jj,ik) 
     165         END_2D        
    166166         !                                   ! bottom scale factors and depth at  U-, V-, UW and VW-points 
    167167         pe3u (:,:,:) = pe3t(:,:,:) 
    168168         pe3uw(:,:,:) = pe3w(:,:,:) 
    169          DO jk = 1, jpk                      ! Computed as the minimum of neighbooring scale factors 
    170             DO jj = 1, jpjm1 
    171                DO ji = 1, jpi 
    172                   pe3v (ji,jj,jk) = MIN( pe3t(ji,jj,jk), pe3t(ji,jj+1,jk) ) 
    173                   pe3vw(ji,jj,jk) = MIN( pe3w(ji,jj,jk), pe3w(ji,jj+1,jk) ) 
    174                   pe3f (ji,jj,jk) = pe3v(ji,jj,jk) 
    175                END DO 
    176             END DO 
    177          END DO 
     169         DO_3D_00_00( 1, jpk ) 
     170         !                                   ! Computed as the minimum of neighbooring scale factors 
     171            pe3v (ji,jj,jk) = MIN( pe3t(ji,jj,jk), pe3t(ji,jj+1,jk) ) 
     172            pe3vw(ji,jj,jk) = MIN( pe3w(ji,jj,jk), pe3w(ji,jj+1,jk) ) 
     173            pe3f (ji,jj,jk) = pe3v(ji,jj,jk) 
     174         END_3D 
    178175         CALL lbc_lnk( 'usrdef_zgr', pe3v , 'V', 1._wp )   ;   CALL lbc_lnk( 'usrdef_zgr', pe3vw, 'V', 1._wp ) 
    179176         CALL lbc_lnk( 'usrdef_zgr', pe3f , 'F', 1._wp ) 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/LOCK_EXCHANGE/MY_SRC/usrdef_hgr.F90

    r10074 r13197  
    2626   PUBLIC   usr_def_hgr   ! called by domhgr.F90 
    2727 
     28   !! * Substitutions 
     29#  include "do_loop_substitute.h90" 
    2830   !!---------------------------------------------------------------------- 
    2931   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    7274      !                       !==  grid point position  ==!   (in kilometers) 
    7375      zfact = rn_dx * 1.e-3         ! conversion in km 
    74       DO jj = 1, jpj 
    75          DO ji = 1, jpi             ! longitude 
    76             plamt(ji,jj) = zfact * (  - 0.5 + REAL( ji-1 + nimpp-1 , wp )  )   
    77             plamu(ji,jj) = zfact * (          REAL( ji-1 + nimpp-1 , wp )  ) 
    78             plamv(ji,jj) = plamt(ji,jj) 
    79             plamf(ji,jj) = plamu(ji,jj) 
    80             !                       ! latitude 
    81             pphit(ji,jj) = zfact * (  - 0.5 + REAL( jj-1 + njmpp-1 , wp )  ) 
    82             pphiu(ji,jj) = pphit(ji,jj) 
    83             pphiv(ji,jj) = zfact * (          REAL( jj-1 + njmpp-1 , wp )  ) 
    84             pphif(ji,jj) = pphiv(ji,jj) 
    85          END DO 
    86       END DO 
     76      DO_2D_11_11 
     77         !                       ! longitude 
     78         plamt(ji,jj) = zfact * (  - 0.5 + REAL( ji-1 + nimpp-1 , wp )  )   
     79         plamu(ji,jj) = zfact * (          REAL( ji-1 + nimpp-1 , wp )  ) 
     80         plamv(ji,jj) = plamt(ji,jj) 
     81         plamf(ji,jj) = plamu(ji,jj) 
     82         !                       ! latitude 
     83         pphit(ji,jj) = zfact * (  - 0.5 + REAL( jj-1 + njmpp-1 , wp )  ) 
     84         pphiu(ji,jj) = pphit(ji,jj) 
     85         pphiv(ji,jj) = zfact * (          REAL( jj-1 + njmpp-1 , wp )  ) 
     86         pphif(ji,jj) = pphiv(ji,jj) 
     87      END_2D 
    8788      ! 
    8889      !                       !==  Horizontal scale factors  ==!   (in meters)  
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/OVERFLOW/MY_SRC/usrdef_hgr.F90

    r10074 r13197  
    2626   PUBLIC   usr_def_hgr   ! called by domhgr.F90 
    2727 
     28   !! * Substitutions 
     29#  include "do_loop_substitute.h90" 
    2830   !!---------------------------------------------------------------------- 
    2931   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    7274      !                       !==  grid point position  ==!   (in kilometers) 
    7375      zfact = rn_dx * 1.e-3         ! conversion in km 
    74       DO jj = 1, jpj 
    75          DO ji = 1, jpi             ! longitude 
    76             plamt(ji,jj) = zfact * (  - 0.5 + REAL( ji-1 + nimpp-1 , wp )  )   
    77             plamu(ji,jj) = zfact * (          REAL( ji-1 + nimpp-1 , wp )  ) 
    78             plamv(ji,jj) = plamt(ji,jj) 
    79             plamf(ji,jj) = plamu(ji,jj) 
    80             !                       ! latitude 
    81             pphit(ji,jj) = zfact * (  - 0.5 + REAL( jj-1 + njmpp-1 , wp )  ) 
    82             pphiu(ji,jj) = pphit(ji,jj) 
    83             pphiv(ji,jj) = zfact * (          REAL( jj-1 + njmpp-1 , wp )  ) 
    84             pphif(ji,jj) = pphiv(ji,jj) 
    85          END DO 
    86       END DO 
     76      DO_2D_11_11 
     77         !                       ! longitude 
     78         plamt(ji,jj) = zfact * (  - 0.5 + REAL( ji-1 + nimpp-1 , wp )  )   
     79         plamu(ji,jj) = zfact * (          REAL( ji-1 + nimpp-1 , wp )  ) 
     80         plamv(ji,jj) = plamt(ji,jj) 
     81         plamf(ji,jj) = plamu(ji,jj) 
     82         !                       ! latitude 
     83         pphit(ji,jj) = zfact * (  - 0.5 + REAL( jj-1 + njmpp-1 , wp )  ) 
     84         pphiu(ji,jj) = pphit(ji,jj) 
     85         pphiv(ji,jj) = zfact * (          REAL( jj-1 + njmpp-1 , wp )  ) 
     86         pphif(ji,jj) = pphiv(ji,jj) 
     87      END_2D 
    8788      ! 
    8889      !                       !==  Horizontal scale factors  ==!   (in meters)  
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/OVERFLOW/MY_SRC/usrdef_zgr.F90

    r12377 r13197  
    2929   PUBLIC   usr_def_zgr   ! called by domzgr.F90 
    3030 
     31   !! * Substitutions 
     32#  include "do_loop_substitute.h90" 
    3133   !!---------------------------------------------------------------------- 
    3234   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    182184            pe3vw(:,:,jk) = pe3w_1d (jk) 
    183185         END DO 
    184          DO jj = 1, jpj                      ! bottom scale factors and depth at T- and W-points 
    185             DO ji = 1, jpi 
    186                ik = k_bot(ji,jj) 
    187                   pdepw(ji,jj,ik+1) = MIN( zht(ji,jj) , pdepw_1d(ik+1) ) 
    188                   pe3t (ji,jj,ik  ) = pdepw(ji,jj,ik+1) - pdepw(ji,jj,ik) 
    189                   pe3t (ji,jj,ik+1) = pe3t (ji,jj,ik  )  
    190                   ! 
    191                   pdept(ji,jj,ik  ) = pdepw(ji,jj,ik  ) + pe3t (ji,jj,ik  ) * 0.5_wp 
    192                   pdept(ji,jj,ik+1) = pdepw(ji,jj,ik+1) + pe3t (ji,jj,ik+1) * 0.5_wp 
    193                   pe3w (ji,jj,ik+1) = pdept(ji,jj,ik+1) - pdept(ji,jj,ik)              ! = pe3t (ji,jj,ik  ) 
    194             END DO 
    195          END DO          
     186         DO_2D_11_11 
     187            ik = k_bot(ji,jj) 
     188            pdepw(ji,jj,ik+1) = MIN( zht(ji,jj) , pdepw_1d(ik+1) ) 
     189            pe3t (ji,jj,ik  ) = pdepw(ji,jj,ik+1) - pdepw(ji,jj,ik) 
     190            pe3t (ji,jj,ik+1) = pe3t (ji,jj,ik  )  
     191            ! 
     192            pdept(ji,jj,ik  ) = pdepw(ji,jj,ik  ) + pe3t (ji,jj,ik  ) * 0.5_wp 
     193            pdept(ji,jj,ik+1) = pdepw(ji,jj,ik+1) + pe3t (ji,jj,ik+1) * 0.5_wp 
     194            pe3w (ji,jj,ik+1) = pdept(ji,jj,ik+1) - pdept(ji,jj,ik)              ! = pe3t (ji,jj,ik  ) 
     195         END_2D          
    196196         !                                   ! bottom scale factors and depth at  U-, V-, UW and VW-points 
    197197         !                                   ! usually Computed as the minimum of neighbooring scale factors 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/STATION_ASF/EXPREF/launch_sasf.sh

    r13159 r13197  
    11#!/bin/bash 
    22 
    3 ################################################################ 
    4 # 
    5 # Script to launch a set of STATION_ASF simulations 
    6 # 
    7 # L. Brodeau, 2020 
    8 # 
    9 ################################################################ 
     3# NEMO directory where to fetch compiled STATION_ASF nemo.exe + setup: 
     4NEMO_DIR=`pwd | sed -e "s|/tests/STATION_ASF/EXPREF||g"` 
    105 
    11 # What directory inside "tests" actually contains the compiled "nemo.exe" for STATION_ASF ? 
     6echo "Using NEMO_DIR=${NEMO_DIR}" 
     7 
     8# what directory inside "tests" actually contains the compiled test-case? 
    129TC_DIR="STATION_ASF2" 
    1310 
    14 # DATA_IN_DIR => Directory containing sea-surface + atmospheric forcings 
     11# => so the executable to use is: 
     12NEMO_EXE="${NEMO_DIR}/tests/${TC_DIR}/BLD/bin/nemo.exe" 
     13 
     14# Directory where to run the simulation: 
     15WORK_DIR="${HOME}/tmp/STATION_ASF" 
     16 
     17 
     18# FORC_DIR => Directory containing sea-surface + atmospheric forcings 
    1519#             (get it there https://drive.google.com/file/d/1MxNvjhRHmMrL54y6RX7WIaM9-LGl--ZP/): 
    1620if [ `hostname` = "merlat"        ]; then 
    17     DATA_IN_DIR="/MEDIA/data/STATION_ASF/input_data_STATION_ASF_2016-2018" 
     21    FORC_DIR="/MEDIA/data/STATION_ASF/input_data_STATION_ASF_2016-2018" 
    1822elif [ `hostname` = "luitel"        ]; then 
    19     DATA_IN_DIR="/data/gcm_setup/STATION_ASF/input_data_STATION_ASF_2016-2018" 
     23    FORC_DIR="/data/gcm_setup/STATION_ASF/input_data_STATION_ASF_2016-2018" 
    2024elif [ `hostname` = "ige-meom-cal1" ]; then 
    21     DATA_IN_DIR="/mnt/meom/workdir/brodeau/STATION_ASF/input_data_STATION_ASF_2016-2018" 
     25    FORC_DIR="/mnt/meom/workdir/brodeau/STATION_ASF/input_data_STATION_ASF_2016-2018" 
    2226elif [ `hostname` = "salvelinus" ]; then 
    23     DATA_IN_DIR="/opt/data/STATION_ASF/input_data_STATION_ASF_2016-2018" 
     27    FORC_DIR="/opt/data/STATION_ASF/input_data_STATION_ASF_2016-2018" 
    2428else 
    25     echo "Oops! We don't know `hostname` yet! Define 'DATA_IN_DIR' in the script!"; exit  
     29    echo "Boo!"; exit 
    2630fi 
    27  
    28 expdir=`basename ${PWD}`; # we expect "EXPREF" or "EXP00" normally... 
    29  
    30 # NEMOGCM root directory where to fetch compiled STATION_ASF nemo.exe + setup: 
    31 NEMO_WRK_DIR=`pwd | sed -e "s|/tests/STATION_ASF/${expdir}||g"` 
    32  
    33 # Directory where to run the simulation: 
    34 PROD_DIR="${HOME}/tmp/STATION_ASF" 
     31#====================== 
     32mkdir -p ${WORK_DIR} 
    3533 
    3634 
    37 ####### End of normal user configurable section ####### 
     35if [ ! -f ${NEMO_EXE} ]; then echo " Mhhh, no compiled nemo.exe found into ${NEMO_DIR}/tests/STATION_ASF/BLD/bin !"; exit; fi 
    3836 
    39 #================================================================================ 
    40  
    41 # NEMO executable to use is: 
    42 NEMO_EXE="${NEMO_WRK_DIR}/tests/${TC_DIR}/BLD/bin/nemo.exe" 
    43  
    44  
    45 echo "###########################################################" 
    46 echo "#        S T A T I O N   A i r  -  S e a   F l u x        #" 
    47 echo "###########################################################" 
    48 echo 
    49 echo " We shall work in here: ${STATION_ASF_DIR}/" 
    50 echo " NEMOGCM   work    depository is: ${NEMO_WRK_DIR}/" 
    51 echo "   ==> NEMO EXE to use: ${NEMO_EXE}" 
    52 echo " Input forcing data into: ${DATA_IN_DIR}/" 
    53 echo " Production will be done into: ${PROD_DIR}/" 
    54 echo 
    55  
    56 mkdir -p ${PROD_DIR} 
    57  
    58 if [ ! -f ${NEMO_EXE} ]; then echo " Mhhh, no compiled 'nemo.exe' found into `dirname ${NEMO_EXE}` !"; exit; fi 
    59  
    60 echo 
    61 echo " *** Using the following NEMO executable:" 
    62 echo "  ${NEMO_EXE} " 
    63 echo 
    64  
    65 NEMO_EXPREF="${NEMO_WRK_DIR}/tests/STATION_ASF/EXPREF" 
     37NEMO_EXPREF="${NEMO_DIR}/tests/STATION_ASF/EXPREF" 
    6638if [ ! -d ${NEMO_EXPREF} ]; then echo " Mhhh, no EXPREF directory ${NEMO_EXPREF} !"; exit; fi 
    6739 
    68 rsync -avP ${NEMO_EXE}          ${PROD_DIR}/ 
     40rsync -avP ${NEMO_EXE}          ${WORK_DIR}/ 
    6941 
    7042for ff in "context_nemo.xml" "domain_def_nemo.xml" "field_def_nemo-oce.xml" "file_def_nemo-oce.xml" "grid_def_nemo.xml" "iodef.xml" "namelist_ref"; do 
    7143    if [ ! -f ${NEMO_EXPREF}/${ff} ]; then echo " Mhhh, ${ff} not found into ${NEMO_EXPREF} !"; exit; fi 
    72     rsync -avPL ${NEMO_EXPREF}/${ff} ${PROD_DIR}/ 
     44    rsync -avPL ${NEMO_EXPREF}/${ff} ${WORK_DIR}/ 
    7345done 
    7446 
    7547# Copy forcing to work directory: 
    76 rsync -avP ${DATA_IN_DIR}/Station_PAPA_50N-145W*.nc ${PROD_DIR}/ 
     48rsync -avP ${FORC_DIR}/Station_PAPA_50N-145W*.nc ${WORK_DIR}/ 
    7749 
    7850for CASE in "ECMWF" "COARE3p6" "NCAR" "ECMWF-noskin" "COARE3p6-noskin"; do 
     
    8658    scase=`echo "${CASE}" | tr '[:upper:]' '[:lower:]'` 
    8759 
    88     rm -f ${PROD_DIR}/namelist_cfg 
    89     rsync -avPL ${NEMO_EXPREF}/namelist_${scase}_cfg ${PROD_DIR}/namelist_cfg 
     60    rm -f ${WORK_DIR}/namelist_cfg 
     61    rsync -avPL ${NEMO_EXPREF}/namelist_${scase}_cfg ${WORK_DIR}/namelist_cfg 
    9062 
    91     cd ${PROD_DIR}/ 
     63    cd ${WORK_DIR}/ 
    9264    echo 
    9365    echo "Launching NEMO !" 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/STATION_ASF/EXPREF/namelist_coare3p6-noskin_cfg

    r13159 r13197  
    3535   nn_time0    =       0   !  initial time of day in hhmm 
    3636   nn_leapy    =       1   !  Leap year calendar (1) or not (0) 
    37    ln_rstart   =  .false.   !  start from rest (F) or from a restart file (T) 
    38       nn_euler    =    1      !  = 0 : start with forward time step if ln_rstart=T 
     37   ln_rstart   = .false.   !  start from rest (F) or from a restart file (T) 
     38      ln_1st_euler = .false.  !  =T force a start with forward time step (ln_rstart=T) 
    3939      nn_rstctl   = 2      !  restart control ==> activated only if ln_rstart=T 
    4040      !                          !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/STATION_ASF/EXPREF/namelist_coare3p6_cfg

    r13159 r13197  
    3535   nn_time0    =       0   !  initial time of day in hhmm 
    3636   nn_leapy    =       1   !  Leap year calendar (1) or not (0) 
    37    ln_rstart   =  .false.   !  start from rest (F) or from a restart file (T) 
    38       nn_euler    =    1      !  = 0 : start with forward time step if ln_rstart=T 
     37   ln_rstart   = .false.   !  start from rest (F) or from a restart file (T) 
     38      ln_1st_euler = .false.  !  =T force a start with forward time step (ln_rstart=T) 
    3939      nn_rstctl   = 2      !  restart control ==> activated only if ln_rstart=T 
    4040      !                          !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/STATION_ASF/EXPREF/namelist_ecmwf-noskin_cfg

    r13159 r13197  
    3535   nn_time0    =       0   !  initial time of day in hhmm 
    3636   nn_leapy    =       1   !  Leap year calendar (1) or not (0) 
    37    ln_rstart   =  .false.   !  start from rest (F) or from a restart file (T) 
    38       nn_euler    =    1      !  = 0 : start with forward time step if ln_rstart=T 
     37   ln_rstart   = .false.   !  start from rest (F) or from a restart file (T) 
     38      ln_1st_euler = .false.  !  =T force a start with forward time step (ln_rstart=T) 
    3939      nn_rstctl   = 2      !  restart control ==> activated only if ln_rstart=T 
    4040      !                          !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/STATION_ASF/EXPREF/namelist_ecmwf_cfg

    r13159 r13197  
    3535   nn_time0    =       0   !  initial time of day in hhmm 
    3636   nn_leapy    =       1   !  Leap year calendar (1) or not (0) 
    37    ln_rstart   =  .false.   !  start from rest (F) or from a restart file (T) 
    38       nn_euler    =    1      !  = 0 : start with forward time step if ln_rstart=T 
     37   ln_rstart   = .false.   !  start from rest (F) or from a restart file (T) 
     38      ln_1st_euler = .false.  !  =T force a start with forward time step (ln_rstart=T) 
    3939      nn_rstctl   = 2      !  restart control ==> activated only if ln_rstart=T 
    4040      !                          !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/STATION_ASF/EXPREF/namelist_ncar_cfg

    r13159 r13197  
    3535   nn_time0    =       0   !  initial time of day in hhmm 
    3636   nn_leapy    =       1   !  Leap year calendar (1) or not (0) 
    37    ln_rstart   =  .false.   !  start from rest (F) or from a restart file (T) 
    38       nn_euler    =    1      !  = 0 : start with forward time step if ln_rstart=T 
     37   ln_rstart   = .false.   !  start from rest (F) or from a restart file (T) 
     38      ln_1st_euler = .false.  !  =T force a start with forward time step (ln_rstart=T) 
    3939      nn_rstctl   = 2      !  restart control ==> activated only if ln_rstart=T 
    4040      !                          !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/STATION_ASF/EXPREF/plot_station_asf.py

    r12031 r13197  
    5353L_VARL  = [ r'$Q_{lat}$', r'$Q_{sens}$' , r'$Q_{net}$' , r'$Q_{lw}$' , r'$|\tau|$' , r'$\Delta T_{skin}$' ] ; # name of variable in latex mode 
    5454L_VUNT  = [ r'$W/m^2$'  , r'$W/m^2$'    , r'$W/m^2$'   , r'$W/m^2$'  , r'$N/m^2$'  ,      'K'             ] 
    55 L_VMAX  = [     75.     ,     75.       ,    800.      ,     25.     ,    1.2      ,      -0.7            ] 
    56 L_VMIN  = [   -250.     ,   -125.       ,   -400.      ,   -150.     ,    0.       ,       0.7            ] 
     55L_VMAX  = [     75.     ,     75.       ,    800.      ,     25.     ,    1.2      ,       0.7            ] 
     56L_VMIN  = [   -250.     ,   -125.       ,   -400.      ,   -150.     ,    0.       ,      -0.7            ] 
    5757L_ANOM  = [   True      ,    True       ,    True      ,    True     ,   True      ,      False           ] 
    5858 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/STATION_ASF/MY_SRC/diawri.F90

    r12489 r13197  
    3535   USE iom            ! 
    3636   USE ioipsl         ! 
     37 
    3738#if defined key_si3 
    3839   USE ice 
     
    5657 
    5758   !!---------------------------------------------------------------------- 
    58    !! NEMO/SAS 4.0 , NEMO Consortium (2018) 
    59    !! $Id: diawri.F90 10425 2018-12-19 21:54:16Z smasson $ 
     59   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     60   !! $Id: diawri.F90 12493 2020-03-02 07:56:31Z smasson $ 
    6061   !! Software governed by the CeCILL license (see ./LICENSE) 
    6162   !!---------------------------------------------------------------------- 
     
    114115      INTEGER, DIMENSION(2) :: ierr 
    115116      !!---------------------------------------------------------------------- 
    116       ierr = 0 
    117       ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     & 
    118          &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     & 
    119          &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) ) 
    120       ! 
    121       dia_wri_alloc = MAXVAL(ierr) 
    122       CALL mpp_sum( 'diawri', dia_wri_alloc ) 
     117      IF( nn_write == -1 ) THEN 
     118         dia_wri_alloc = 0 
     119      ELSE 
     120         ierr = 0 
     121         ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     & 
     122            &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     & 
     123            &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) ) 
     124         ! 
     125         dia_wri_alloc = MAXVAL(ierr) 
     126         CALL mpp_sum( 'diawri', dia_wri_alloc ) 
     127         ! 
     128      ENDIF 
    123129      ! 
    124130   END FUNCTION dia_wri_alloc 
     
    374380      CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm)                )    ! now i-velocity 
    375381      CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm)                )    ! now j-velocity 
    376          CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww             )    ! now k-velocity 
     382      CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww             )    ! now k-velocity 
    377383      CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf         )    ! freshwater budget 
    378384      CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns         )    ! total heat flux 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/STATION_ASF/MY_SRC/nemogcm.F90

    r13159 r13197  
    22   !!====================================================================== 
    33   !!                       ***  MODULE nemogcm   *** 
    4    !! StandAlone Surface module : surface fluxes 
     4   !!                      STATION_ASF (SAS meets C1D) 
    55   !!====================================================================== 
    66   !! History :  3.6  ! 2011-11  (S. Alderson, G. Madec) original code 
     
    1919   !!---------------------------------------------------------------------- 
    2020   USE step_oce       ! module used in the ocean time stepping module (step.F90) 
    21    USE sbc_oce        ! surface boundary condition: ocean #LB: rm? 
    2221   USE phycst         ! physical constant                  (par_cst routine) 
    2322   USE domain         ! domain initialization   (dom_init & dom_cfg routines) 
    2423   USE closea         ! treatment of closed seas (for ln_closea) 
    2524   USE usrdef_nam     ! user defined configuration 
     25   USE istate         ! initial state setting          (istate_init routine) 
    2626   USE step, ONLY : Nbb, Nnn, Naa, Nrhs ! time level indices 
    2727   USE daymod         ! calendar 
    2828   USE restart        ! open  restart file 
    29    !LB:USE step           ! NEMO time-stepping                 (stp     routine) 
    3029   USE c1d            ! 1D configuration 
    3130   USE step_c1d       ! Time stepping loop for the 1D configuration 
    32    USE sbcssm         ! 
    3331   ! 
     32   USE in_out_manager ! I/O manager 
    3433   USE lib_mpp        ! distributed memory computing 
    3534   USE mppini         ! shared/distributed memory setting (mpp_init routine) 
     
    4948   !!---------------------------------------------------------------------- 
    5049   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    51    !! $Id: nemogcm.F90 11536 2019-09-11 13:54:18Z smasson $ 
     50   !! $Id: nemogcm.F90 12489 2020-02-28 15:55:11Z davestorkey $ 
    5251   !! Software governed by the CeCILL license (see ./LICENSE) 
    5352   !!---------------------------------------------------------------------- 
     
    8483      !                            !==   time stepping   ==! 
    8584      !                            !-----------------------! 
     85      ! 
     86      !                                               !== set the model time-step  ==! 
     87      ! 
    8688      istp = nit000 
    8789      ! 
     
    107109      ! 
    108110#if defined key_iomput 
    109       CALL xios_finalize  ! end mpp communications with xios 
     111                                    CALL xios_finalize  ! end mpp communications with xios 
    110112#else 
    111       IF( lk_mpp   ) THEN   ;   CALL mppstop      ! end mpp communications 
    112       ENDIF 
     113      IF( lk_mpp )                  CALL mppstop      ! end mpp communications 
    113114#endif 
    114115      ! 
     
    162163      IF( lwm )   CALL ctl_opn(     numout,        'ocean.output', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, -1, .FALSE. ) 
    163164      ! open reference and configuration namelist files 
    164       CALL load_nml( numnam_ref,        'namelist_ref',                                           -1, lwm ) 
    165       CALL load_nml( numnam_cfg,        'namelist_cfg',                                           -1, lwm ) 
     165                  CALL load_nml( numnam_ref,        'namelist_ref',                                           -1, lwm ) 
     166                  CALL load_nml( numnam_cfg,        'namelist_cfg',                                           -1, lwm ) 
    166167      IF( lwm )   CALL ctl_opn(     numond, 'output.namelist.dyn', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, -1, .FALSE. ) 
    167168      ! open /dev/null file to be able to supress output write easily 
    168       CALL ctl_opn(     numnul,           '/dev/null', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, -1, .FALSE. ) 
     169      IF( Agrif_Root() ) THEN 
     170                  CALL ctl_opn(     numnul,           '/dev/null', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, -1, .FALSE. ) 
     171#ifdef key_agrif 
     172      ELSE 
     173                  numnul = Agrif_Parent(numnul)    
     174#endif 
     175      ENDIF 
    169176      ! 
    170177      !                             !--------------------! 
     
    222229903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namcfg in reference namelist' ) 
    223230      READ  ( numnam_cfg, namcfg, IOSTAT = ios, ERR = 904 ) 
    224 904   IF( ios >  0 )   CALL ctl_nam ( ios , 'namcfg in configuration namelist' ) 
     231904   IF( ios >  0 )   CALL ctl_nam ( ios , 'namcfg in configuration namelist' )    
    225232      ! 
    226233      IF( ln_read_cfg ) THEN            ! Read sizes in domain configuration file 
     
    253260      IF( ln_timing    )   CALL timing_start( 'nemo_init') 
    254261      ! 
    255       CALL     phy_cst         ! Physical constants 
    256       CALL     eos_init        ! Equation of state 
     262                           CALL     phy_cst         ! Physical constants 
     263                           CALL     eos_init        ! Equation of state 
    257264      IF( lk_c1d       )   CALL     c1d_init        ! 1D column configuration 
    258       CALL     dom_init( Nbb, Nnn, Naa, "OPA") ! Domain 
     265                           CALL     dom_init( Nbb, Nnn, Naa, "OPA") ! Domain 
    259266      IF( sn_cfctl%l_prtctl )   & 
    260267         &                 CALL prt_ctl_init        ! Print control 
    261  
    262       IF( ln_rstart ) THEN                    ! Restart from a file                                                                                  
    263          !                                    ! -------------------                                                                                  
    264          CALL rst_read( Nbb, Nnn )            ! Read the restart file                                                                                
    265          CALL day_init                        ! model calendar (using both namelist and restart infos)                                               
    266          !                                                                                                                                           
    267       ELSE                                    ! Start from rest                                                                                      
    268          !                                    ! ---------------                                                                                      
    269          numror = 0                           ! define numror = 0 -> no restart file to read                                                         
    270          neuler = 0                           ! Set time-step indicator at nit000 (euler forward)                                                    
    271          CALL day_init                        ! model calendar (using both namelist and restart infos)                                               
    272       ENDIF 
    273       ! 
    274  
    275       !                                      ! external forcing 
    276       CALL     sbc_init( Nbb, Nnn, Naa )    ! surface boundary conditions (including sea-ice) 
     268      ! 
     269       
     270                           CALL  istate_init( Nbb, Nnn, Naa )    ! ocean initial state (Dynamics and tracers) 
     271 
     272      !                                      ! external forcing  
     273                           CALL     sbc_init( Nbb, Nnn, Naa )    ! surface boundary conditions (including sea-ice) 
    277274 
    278275      ! 
     
    305302         WRITE(numout,*) '                              sn_cfctl%l_prttrc  = ', sn_cfctl%l_prttrc 
    306303         WRITE(numout,*) '                              sn_cfctl%l_oasout  = ', sn_cfctl%l_oasout 
    307          WRITE(numout,*) '                              sn_cfctl%procmin   = ', sn_cfctl%procmin 
    308          WRITE(numout,*) '                              sn_cfctl%procmax   = ', sn_cfctl%procmax 
    309          WRITE(numout,*) '                              sn_cfctl%procincr  = ', sn_cfctl%procincr 
    310          WRITE(numout,*) '                              sn_cfctl%ptimincr  = ', sn_cfctl%ptimincr 
     304         WRITE(numout,*) '                              sn_cfctl%procmin   = ', sn_cfctl%procmin   
     305         WRITE(numout,*) '                              sn_cfctl%procmax   = ', sn_cfctl%procmax   
     306         WRITE(numout,*) '                              sn_cfctl%procincr  = ', sn_cfctl%procincr  
     307         WRITE(numout,*) '                              sn_cfctl%ptimincr  = ', sn_cfctl%ptimincr  
    311308         WRITE(numout,*) '      level of print                  nn_print   = ', nn_print 
    312309         WRITE(numout,*) '      Start i indice for SUM control  nn_ictls   = ', nn_ictls 
     
    423420      !!---------------------------------------------------------------------- 
    424421      ! 
    425       ierr =        oce_alloc    ()    ! ocean 
     422      ierr =        oce_alloc    ()    ! ocean  
    426423      ierr = ierr + dia_wri_alloc() 
    427424      ierr = ierr + dom_oce_alloc()    ! ocean domain 
     
    432429   END SUBROUTINE nemo_alloc 
    433430 
    434  
     431    
    435432   SUBROUTINE nemo_set_cfctl(sn_cfctl, setto ) 
    436433      !!---------------------------------------------------------------------- 
     
    456453   !!====================================================================== 
    457454END MODULE nemogcm 
     455 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/STATION_ASF/MY_SRC/sbcssm.F90

    r12249 r13197  
    5454   !!---------------------------------------------------------------------- 
    5555   !! NEMO/SAS 4.0 , NEMO Consortium (2018) 
    56    !! $Id: sbcssm.F90 10068 2018-08-28 14:09:04Z nicolasmartin $ 
     56   !! $Id: sbcssm.F90 12615 2020-03-26 15:18:49Z laurent $ 
    5757   !! Software governed by the CeCILL license (see ./LICENSE) 
    5858   !!---------------------------------------------------------------------- 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/STATION_ASF/MY_SRC/usrdef_hgr.F90

    r11930 r13197  
    1414   !!   usr_def_hgr   : initialize the horizontal mesh  
    1515   !!---------------------------------------------------------------------- 
    16    USE dom_oce  , ONLY: nimpp, njmpp        ! ocean space and time domain 
     16   USE dom_oce  , ONLY: nimpp, njmpp       ! ocean space and time domain 
    1717   USE c1d      ,  ONLY: rn_lon1d, rn_lat1d ! ocean lon/lat define by namelist 
    1818   USE par_oce        ! ocean space and time domain 
     
    3030   !!---------------------------------------------------------------------- 
    3131   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    32    !! $Id: usrdef_hgr.F90 10072 2018-08-28 15:21:50Z nicolasmartin $ 
     32   !! $Id: usrdef_hgr.F90 12489 2020-02-28 15:55:11Z davestorkey $  
    3333   !! Software governed by the CeCILL license (see ./LICENSE) 
    3434   !!---------------------------------------------------------------------- 
     
    5454      !! 
    5555      !! ** Action  : - define longitude & latitude of t-, u-, v- and f-points (in degrees)  
    56       !!              - define coriolis parameter at f-point if the domain in not on the sphere 
     56      !!              - define coriolis parameter at f-point if the domain in not on the sphere (on beta-plane) 
    5757      !!              - define i- & j-scale factors at t-, u-, v- and f-points (in meters) 
    5858      !!              - define u- & v-surfaces (if gridsize reduction is used in some straits) (in m2) 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/STATION_ASF/MY_SRC/usrdef_nam.F90

    r12249 r13197  
    88   !!====================================================================== 
    99   !! History :  4.0  ! 2016-03  (S. Flavoni, G. Madec)  Original code 
    10    !! History :  4.x  ! 2019-10  (L. Brodeau) for STATION_ASF (C1D meets SAS) 
     10   !!            4.x  ! 2019-10  (L. Brodeau) for STATION_ASF (C1D meets SAS) 
    1111   !!---------------------------------------------------------------------- 
    1212 
     
    1515   !!   usr_def_hgr   : initialize the horizontal mesh  
    1616   !!---------------------------------------------------------------------- 
    17    USE dom_oce  , ONLY: nimpp, njmpp             ! ocean space and time domain 
    18    USE dom_oce  , ONLY: ln_zco, ln_zps, ln_sco   ! flag of type of coordinate 
     17   USE dom_oce  , ONLY: nimpp, njmpp       ! ocean space and time domain 
     18!!!   USE dom_oce  , ONLY: ln_zco, ln_zps, ln_sco   ! flag of type of coordinate 
    1919   USE par_oce        ! ocean space and time domain 
    2020   USE phycst         ! physical constants 
     
    3333   !!---------------------------------------------------------------------- 
    3434   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    35    !! $Id: usrdef_nam.F90 11536 2019-09-11 13:54:18Z smasson $  
     35   !! $Id: usrdef_nam.F90 12377 2020-02-12 14:39:06Z acc $  
    3636   !! Software governed by the CeCILL license (see ./LICENSE) 
    3737   !!---------------------------------------------------------------------- 
     
    6868      kk_cfg = 0 
    6969 
    70       ! Global Domain size: STATION_ASF domain is 3 x 3 grid-points x 75 or vertical levels 
     70      ! Global Domain size: STATION_ASF domain is 3 x 3 grid-points x 2 or vertical levels 
    7171      kpi = 3 
    7272      kpj = 3 
    73       kpk = 1 
     73      kpk = 2    ! 2, rather than 1, because 1 would cause some issues... like overflow in array boundary indexes, etc... 
    7474      ! 
    7575      !                             ! Set the lateral boundary condition of the global domain 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/STATION_ASF/MY_SRC/usrdef_zgr.F90

    r12038 r13197  
    11MODULE usrdef_zgr 
    22   !!====================================================================== 
    3    !!                     ***  MODULE usrdef_zgr  *** 
     3   !!                       ***  MODULE usrdef_zgr  *** 
    44   !! 
    55   !!                       ===  STATION_ASF case  === 
    66   !! 
    7    !! user defined : vertical coordinate system of a user configuration 
     7   !! User defined : vertical coordinate system of a user configuration 
    88   !!====================================================================== 
    9    !! History :  4.0  ! 2019-10  (L. Brodeau)  Original code 
     9   !! History :  4.0  ! 2016-06  (G. Madec)  Original code 
     10   !!            4.x  ! 2019-10  (L. Brodeau) Station ASF 
    1011   !!---------------------------------------------------------------------- 
    1112 
    1213   !!---------------------------------------------------------------------- 
    13    !!   usr_def_zgr   : user defined vertical coordinate system (required) 
     14   !!   usr_def_zgr   : user defined vertical coordinate system 
     15   !!      zgr_z      : reference 1D z-coordinate  
     16   !!      zgr_top_bot: ocean top and bottom level indices 
     17   !!      zgr_zco    : 3D verticl coordinate in pure z-coordinate case 
    1418   !!--------------------------------------------------------------------- 
    1519   USE oce            ! ocean variables 
    16    !USE dom_oce        ! ocean domain 
    17    !USE depth_e3       ! depth <=> e3 
    1820   USE usrdef_nam     ! User defined : namelist variables 
    1921   ! 
     
    2123   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    2224   USE lib_mpp        ! distributed memory computing library 
    23    USE timing         ! Timing 
    2425 
    2526   IMPLICIT NONE 
    2627   PRIVATE 
    2728 
    28    PUBLIC   usr_def_zgr   ! called by domzgr.F90 
     29   PUBLIC   usr_def_zgr        ! called by domzgr.F90 
    2930 
    3031   !!---------------------------------------------------------------------- 
    3132   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    32    !! $Id: usrdef_zgr.F90 10072 2018-08-28 15:21:50Z nicolasmartin $ 
     33   !! $Id: usrdef_zgr.F90 12377 2020-02-12 14:39:06Z acc $ 
    3334   !! Software governed by the CeCILL license (see ./LICENSE) 
    3435   !!---------------------------------------------------------------------- 
     
    4748      !! 
    4849      !!---------------------------------------------------------------------- 
    49       LOGICAL                   , INTENT(  out) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags ( read in namusr_def ) 
    50       LOGICAL                   , INTENT(  out) ::   ld_isfcav                   ! under iceshelf cavity flag 
    51       REAL(wp), DIMENSION(:)    , INTENT(  out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth     [m] 
    52       REAL(wp), DIMENSION(:)    , INTENT(  out) ::   pe3t_1d , pe3w_1d           ! 1D grid-point depth     [m] 
    53       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdept, pdepw                ! grid-point depth        [m] 
    54       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors  [m] 
    55       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3w , pe3uw, pe3vw         ! i-scale factors 
    56       INTEGER , DIMENSION(:,:)  , INTENT(  out) ::   k_top, k_bot                ! first & last ocean level 
     50      LOGICAL                   , INTENT(out) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags 
     51      LOGICAL                   , INTENT(out) ::   ld_isfcav                   ! under iceshelf cavity flag 
     52      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth     [m] 
     53      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d           ! 1D grid-point depth     [m] 
     54      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pdept, pdepw                ! grid-point depth        [m] 
     55      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors  [m] 
     56      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3w , pe3uw, pe3vw         ! i-scale factors  
     57      INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_top, k_bot                ! first & last ocean level 
    5758      !!---------------------------------------------------------------------- 
    5859      ! 
     
    6162      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    6263      ! 
    63  
     64      ! 
     65      ! type of vertical coordinate 
     66      ! --------------------------- 
    6467      ld_zco    = .TRUE.         ! z-coordinate without ocean cavities 
    6568      ld_zps    = .FALSE. 
    6669      ld_sco    = .FALSE. 
    6770      ld_isfcav = .FALSE. 
    68        
     71 
     72      !! 1st level (the only one that matters) 
    6973      pdept_1d(1) = rn_dept1 ! depth (m) at which the SST is taken/measured == depth of first T point! 
    7074      pdepw_1d(1) = 0._wp 
     
    7276      pe3w_1d(1)  = rn_dept1 ! LB??? 
    7377 
    74       pdept(:,:,:) = rn_dept1 
    75       pdepw(:,:,:) = 0._wp 
    76       pe3t(:,:,:) = 2._wp*rn_dept1 
    77       pe3u(:,:,:) = 2._wp*rn_dept1 
    78       pe3v(:,:,:) = 2._wp*rn_dept1 
    79       pe3f(:,:,:) = 2._wp*rn_dept1 
    80       pe3w(:,:,:)  = rn_dept1  ! LB??? 
    81       pe3uw(:,:,:) = rn_dept1  ! LB??? 
    82       pe3vw(:,:,:) = rn_dept1  ! LB??? 
     78      pdept(:,:,1) = rn_dept1 
     79      pdepw(:,:,1) = 0._wp 
     80      pe3t(:,:,1) = 2._wp*rn_dept1 
     81      pe3u(:,:,1) = 2._wp*rn_dept1 
     82      pe3v(:,:,1) = 2._wp*rn_dept1 
     83      pe3f(:,:,1) = 2._wp*rn_dept1 
     84      pe3w(:,:,1)  = rn_dept1  ! LB??? 
     85      pe3uw(:,:,1) = rn_dept1  ! LB??? 
     86      pe3vw(:,:,1) = rn_dept1  ! LB??? 
     87       
     88      !! 2nd level, technically useless (only for the sake of code stability) 
     89      pdept_1d(2) = 3._wp*rn_dept1 
     90      pdepw_1d(2) = 2._wp*rn_dept1 
     91      pe3t_1d(2)  = 2._wp*rn_dept1 
     92      pe3w_1d(2)  = 2._wp*rn_dept1 
     93 
     94      pdept(:,:,2) = 3._wp*rn_dept1 
     95      pdepw(:,:,2) = 2._wp*rn_dept1 
     96      pe3t(:,:,2) = 2._wp*rn_dept1 
     97      pe3u(:,:,2) = 2._wp*rn_dept1 
     98      pe3v(:,:,2) = 2._wp*rn_dept1 
     99      pe3f(:,:,2) = 2._wp*rn_dept1 
     100      pe3w(:,:,2)  = 2._wp*rn_dept1 
     101      pe3uw(:,:,2) = 2._wp*rn_dept1 
     102      pe3vw(:,:,2) = 2._wp*rn_dept1 
     103 
    83104      k_top = 1 
    84105      k_bot = 1 
    85       ! 
     106 
    86107   END SUBROUTINE usr_def_zgr 
    87108   !!====================================================================== 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/VORTEX/MY_SRC/domvvl.F90

    r12489 r13197  
    6363   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! retoring period for low freq. divergence 
    6464 
     65   !! * Substitutions 
     66#  include "do_loop_substitute.h90" 
    6567   !!---------------------------------------------------------------------- 
    6668   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    188190      gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb) 
    189191      gdepw(:,:,1,Kbb) = 0.0_wp 
    190       DO jk = 2, jpk                               ! vertical sum 
    191          DO jj = 1,jpj 
    192             DO ji = 1,jpi 
    193                !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
    194                !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
    195                !                             ! 0.5 where jk = mikt      
     192      DO_3D_11_11( 2, jpk ) 
     193         !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
     194         !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
     195         !                             ! 0.5 where jk = mikt      
    196196!!gm ???????   BUG ?  gdept(:,:,:,Kmm) as well as gde3w  does not include the thickness of ISF ?? 
    197                zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 
    198                gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
    199                gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  & 
    200                   &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm))  
    201                gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
    202                gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 
    203                gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  & 
    204                   &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb))  
    205             END DO 
    206          END DO 
    207       END DO 
     197         zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 
     198         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
     199         gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  & 
     200            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm))  
     201         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
     202         gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 
     203         gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  & 
     204            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb))  
     205      END_3D 
    208206      ! 
    209207      !                    !==  thickness of the water column  !!   (ocean portion only) 
     
    240238         ENDIF 
    241239         IF ( ln_vvl_zstar_at_eqtor ) THEN   ! use z-star in vicinity of the Equator 
    242             DO jj = 1, jpj 
    243                DO ji = 1, jpi 
     240            DO_2D_11_11 
    244241!!gm  case |gphi| >= 6 degrees is useless   initialized just above by default 
    245                   IF( ABS(gphit(ji,jj)) >= 6.) THEN 
    246                      ! values outside the equatorial band and transition zone (ztilde) 
    247                      frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp ) 
    248                      frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 
    249                   ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star 
    250                      ! values inside the equatorial band (ztilde as zstar) 
    251                      frq_rst_e3t(ji,jj) =  0.0_wp 
    252                      frq_rst_hdv(ji,jj) =  1.0_wp / rn_Dt 
    253                   ELSE                                      ! transition band (2.5 to 6 degrees N/S) 
    254                      !                                      ! (linearly transition from z-tilde to z-star) 
    255                      frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   & 
    256                         &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
    257                         &                                          * 180._wp / 3.5_wp ) ) 
    258                      frq_rst_hdv(ji,jj) = (1.0_wp / rn_Dt)                                & 
    259                         &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rn_Dt) )*0.5_wp   & 
    260                         &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
    261                         &                                          * 180._wp / 3.5_wp ) ) 
    262                   ENDIF 
    263                END DO 
    264             END DO 
     242               IF( ABS(gphit(ji,jj)) >= 6.) THEN 
     243                  ! values outside the equatorial band and transition zone (ztilde) 
     244                  frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp ) 
     245                  frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 
     246               ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star 
     247                  ! values inside the equatorial band (ztilde as zstar) 
     248                  frq_rst_e3t(ji,jj) =  0.0_wp 
     249                  frq_rst_hdv(ji,jj) =  1.0_wp / rn_Dt 
     250               ELSE                                      ! transition band (2.5 to 6 degrees N/S) 
     251                  !                                      ! (linearly transition from z-tilde to z-star) 
     252                  frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   & 
     253                     &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
     254                     &                                          * 180._wp / 3.5_wp ) ) 
     255                  frq_rst_hdv(ji,jj) = (1.0_wp / rn_Dt)                                & 
     256                     &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rn_Dt) )*0.5_wp   & 
     257                     &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
     258                     &                                          * 180._wp / 3.5_wp ) ) 
     259               ENDIF 
     260            END_2D 
    265261            IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN 
    266262               IF( nn_cfg == 3 ) THEN   ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 
     
    357353      END DO 
    358354      ! 
    359       IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate ! 
    360          !                                                            ! ------baroclinic part------ ! 
     355      IF( (ln_vvl_ztilde .OR. ln_vvl_layer) .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate ! 
     356         !                                                               ! ------baroclinic part------ ! 
    361357         ! I - initialization 
    362358         ! ================== 
     
    411407         zwu(:,:) = 0._wp 
    412408         zwv(:,:) = 0._wp 
    413          DO jk = 1, jpkm1        ! a - first derivative: diffusive fluxes 
    414             DO jj = 1, jpjm1 
    415                DO ji = 1, jpim1   ! vector opt. 
    416                   un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           & 
    417                      &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
    418                   vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           &  
    419                      &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
    420                   zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
    421                   zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 
    422                END DO 
    423             END DO 
    424          END DO 
    425          DO jj = 1, jpj          ! b - correction for last oceanic u-v points 
    426             DO ji = 1, jpi 
    427                un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 
    428                vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 
    429             END DO 
    430          END DO 
    431          DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes 
    432             DO jj = 2, jpjm1 
    433                DO ji = 2, jpim1   ! vector opt. 
    434                   tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
    435                      &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    & 
    436                      &                                            ) * r1_e1e2t(ji,jj) 
    437                END DO 
    438             END DO 
    439          END DO 
     409         DO_3D_10_10( 1, jpkm1 ) 
     410            un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           & 
     411               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
     412            vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           &  
     413               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
     414            zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
     415            zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 
     416         END_3D 
     417         DO_2D_11_11 
     418            un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 
     419            vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 
     420         END_2D 
     421         DO_3D_00_00( 1, jpkm1 ) 
     422            tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
     423               &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    & 
     424               &                                            ) * r1_e1e2t(ji,jj) 
     425         END_3D 
    440426         !                       ! d - thickness diffusion transport: boundary conditions 
    441427         !                             (stored for tracer advction and continuity equation) 
     
    444430         ! 4 - Time stepping of baroclinic scale factors 
    445431         ! --------------------------------------------- 
    446          ! Leapfrog time stepping 
    447          ! ~~~~~~~~~~~~~~~~~~~~~~ 
    448432         CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp ) 
    449433         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + rDt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 
     
    646630      ! Horizontal scale factor interpolations 
    647631      ! -------------------------------------- 
    648       ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are allready computed in dynnxt 
     632      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt 
    649633      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also 
    650634       
     
    663647      gdepw(:,:,1,Kmm) = 0.0_wp 
    664648      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 
    665       DO jk = 2, jpk 
    666          DO jj = 1,jpj 
    667             DO ji = 1,jpi 
    668               !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
    669                                                                  ! 1 for jk = mikt 
    670                zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    671                gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
    672                gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  & 
    673                    &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) )  
    674                gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
    675             END DO 
    676          END DO 
    677       END DO 
     649      DO_3D_11_11( 2, jpk ) 
     650        !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
     651                                                           ! 1 for jk = mikt 
     652         zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
     653         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
     654         gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  & 
     655             &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) )  
     656         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
     657      END_3D 
    678658 
    679659      ! Local depth and Inverse of the local depth of the water 
     
    722702         ! 
    723703      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean 
    724          DO jk = 1, jpk 
    725             DO jj = 1, jpjm1 
    726                DO ji = 1, jpim1   ! vector opt. 
    727                   pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   & 
    728                      &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
    729                      &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
    730                END DO 
    731             END DO 
    732          END DO 
     704         DO_3D_10_10( 1, jpk ) 
     705            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   & 
     706               &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
     707               &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
     708         END_3D 
    733709         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp ) 
    734710         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 
    735711         ! 
    736712      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean 
    737          DO jk = 1, jpk 
    738             DO jj = 1, jpjm1 
    739                DO ji = 1, jpim1   ! vector opt. 
    740                   pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   & 
    741                      &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
    742                      &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
    743                END DO 
    744             END DO 
    745          END DO 
     713         DO_3D_10_10( 1, jpk ) 
     714            pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   & 
     715               &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
     716               &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
     717         END_3D 
    746718         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp ) 
    747719         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 
    748720         ! 
    749721      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean 
    750          DO jk = 1, jpk 
    751             DO jj = 1, jpjm1 
    752                DO ji = 1, jpim1   ! vector opt. 
    753                   pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
    754                      &                       *    r1_e1e2f(ji,jj)                                                  & 
    755                      &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
    756                      &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
    757                END DO 
    758             END DO 
    759          END DO 
     722         DO_3D_10_10( 1, jpk ) 
     723            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
     724               &                       *    r1_e1e2f(ji,jj)                                                  & 
     725               &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
     726               &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
     727         END_3D 
    760728         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp ) 
    761729         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 
     
    832800            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 
    833801            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 
     802            ! 
    834803            !                             ! --------- ! 
    835804            !                             ! all cases ! 
    836805            !                             ! --------- ! 
     806            ! 
    837807            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist 
    838808               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 
     
    850820               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files' 
    851821               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 
    852                IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.' 
     822               IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    853823               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 
    854824               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 
     
    857827               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files' 
    858828               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 
    859                IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.' 
     829               IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    860830               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 
    861831               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
     
    864834               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file' 
    865835               IF(lwp) write(numout,*) 'Compute scale factor from sshn' 
    866                IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.' 
     836               IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    867837               DO jk = 1, jpk 
    868838                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 
     
    917887                  ssh(:,:,Kbb) = -ssh_ref 
    918888 
    919                   DO jj = 1, jpj 
    920                      DO ji = 1, jpi 
    921                         IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth 
    922                            ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 
    923                            ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 
    924                         ENDIF 
    925                      ENDDO 
    926                   ENDDO 
     889                  DO_2D_11_11 
     890                     IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth 
     891                        ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 
     892                        ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 
     893                     ENDIF 
     894                  END_2D 
    927895               ENDIF !If test case else 
    928896 
     
    935903               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    936904 
    937                DO ji = 1, jpi 
    938                   DO jj = 1, jpj 
    939                      IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 
    940                        CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 
    941                      ENDIF 
    942                   END DO  
    943                END DO  
     905               DO_2D_11_11 
     906                  IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 
     907                     CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 
     908                  ENDIF 
     909               END_2D 
    944910               ! 
    945911            ELSE 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/VORTEX/MY_SRC/usrdef_hgr.F90

    r10074 r13197  
    2626   PUBLIC   usr_def_hgr   ! called by domhgr.F90 
    2727 
     28   !! * Substitutions 
     29#  include "do_loop_substitute.h90" 
    2830   !!---------------------------------------------------------------------- 
    2931   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    8890#endif 
    8991          
    90       DO jj = 1, jpj 
    91          DO ji = 1, jpi 
    92             zti = FLOAT( ji - 1 + nimpp - 1 )          ;  ztj = FLOAT( jj - 1 + njmpp - 1 ) 
    93             zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5_wp ;  zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5_wp 
    94  
    95             plamt(ji,jj) = zlam0 + rn_dx * 1.e-3 * zti 
    96             plamu(ji,jj) = zlam0 + rn_dx * 1.e-3 * zui 
    97             plamv(ji,jj) = plamt(ji,jj)  
    98             plamf(ji,jj) = plamu(ji,jj)  
    99     
    100             pphit(ji,jj) = zphi0 + rn_dy * 1.e-3 * ztj 
    101             pphiv(ji,jj) = zphi0 + rn_dy * 1.e-3 * zvj 
    102             pphiu(ji,jj) = pphit(ji,jj)  
    103             pphif(ji,jj) = pphiv(ji,jj)  
    104          END DO 
    105       END DO 
     92      DO_2D_11_11 
     93         zti = FLOAT( ji - 1 + nimpp - 1 )          ;  ztj = FLOAT( jj - 1 + njmpp - 1 ) 
     94         zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5_wp ;  zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5_wp 
     95          
     96         plamt(ji,jj) = zlam0 + rn_dx * 1.e-3 * zti 
     97         plamu(ji,jj) = zlam0 + rn_dx * 1.e-3 * zui 
     98         plamv(ji,jj) = plamt(ji,jj)  
     99         plamf(ji,jj) = plamu(ji,jj)  
     100          
     101         pphit(ji,jj) = zphi0 + rn_dy * 1.e-3 * ztj 
     102         pphiv(ji,jj) = zphi0 + rn_dy * 1.e-3 * zvj 
     103         pphiu(ji,jj) = pphit(ji,jj)  
     104         pphif(ji,jj) = pphiv(ji,jj)  
     105      END_2D 
    106106      !      
    107107      ! Horizontal scale factors (in meters) 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/VORTEX/MY_SRC/usrdef_istate.F90

    r12489 r13197  
    2828   PUBLIC   usr_def_istate   ! called by istate.F90 
    2929 
     30   !! * Substitutions 
     31#  include "do_loop_substitute.h90" 
    3032   !!---------------------------------------------------------------------- 
    3133   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    7375      ! Sea level: 
    7476      za = -zP0 * (1._wp-EXP(-zH)) / (grav*(zH-1._wp + EXP(-zH))) 
    75       DO ji=1, jpi 
    76          DO jj=1, jpj 
    77             zx = glamt(ji,jj) * 1.e3 
    78             zy = gphit(ji,jj) * 1.e3 
    79             zrho1 = rho0 + za * EXP(-(zx**2+zy**2)/zlambda**2) 
    80             pssh(ji,jj) = zP0 * EXP(-(zx**2+zy**2)/zlambda**2)/(zrho1*grav) * ptmask(ji,jj,1) 
    81          END DO 
    82       END DO 
     77      DO_2D_11_11 
     78         zx = glamt(ji,jj) * 1.e3 
     79         zy = gphit(ji,jj) * 1.e3 
     80         zrho1 = rho0 + za * EXP(-(zx**2+zy**2)/zlambda**2) 
     81         pssh(ji,jj) = zP0 * EXP(-(zx**2+zy**2)/zlambda**2)/(zrho1*grav) * ptmask(ji,jj,1) 
     82      END_2D 
    8383      ! 
    8484      ! temperature:          
    85       DO ji=1, jpi 
    86          DO jj=1, jpj 
    87             zx = glamt(ji,jj) * 1.e3 
    88             zy = gphit(ji,jj) * 1.e3 
    89             DO jk=1,jpk 
    90                zdt =  pdept(ji,jj,jk)  
    91                zrho1 = rho0 * (1._wp + zn2*zdt/grav) 
    92                IF (zdt < zH) THEN 
    93                   zrho1 = zrho1 - zP0 * (1._wp-EXP(zdt-zH)) & 
    94                           & * EXP(-(zx**2+zy**2)/zlambda**2) / (grav*(zH -1._wp + exp(-zH))); 
    95                ENDIF 
    96                pts(ji,jj,jk,jp_tem) = (20._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 
    97             END DO 
     85      DO_2D_11_11 
     86         zx = glamt(ji,jj) * 1.e3 
     87         zy = gphit(ji,jj) * 1.e3 
     88         DO jk=1,jpk 
     89            zdt =  pdept(ji,jj,jk)  
     90            zrho1 = rho0 * (1._wp + zn2*zdt/grav) 
     91            IF (zdt < zH) THEN 
     92               zrho1 = zrho1 - zP0 * (1._wp-EXP(zdt-zH)) & 
     93                  & * EXP(-(zx**2+zy**2)/zlambda**2) / (grav*(zH -1._wp + EXP(-zH))); 
     94            ENDIF 
     95            pts(ji,jj,jk,jp_tem) = (20._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 
    9896         END DO 
    99       END DO 
     97      END_2D 
    10098      ! 
    10199      ! salinity:   
     
    104102      ! velocities: 
    105103      za = 2._wp * zP0 / (zf0 * rho0 * zlambda**2) 
    106       DO ji=1, jpim1 
    107          DO jj=1, jpj 
    108             zx = glamu(ji,jj) * 1.e3 
    109             zy = gphiu(ji,jj) * 1.e3 
    110             DO jk=1, jpk 
    111                zdu = 0.5_wp * (pdept(ji  ,jj,jk) + pdept(ji+1,jj,jk)) 
    112                IF (zdu < zH) THEN 
    113                   zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH)) 
    114                   pu(ji,jj,jk) = (za * zf * zy * EXP(-(zx**2+zy**2)/zlambda**2)) * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk) 
    115                ELSE 
    116                   pu(ji,jj,jk) = 0._wp 
    117                ENDIF 
    118             END DO 
     104      DO_2D_00_00 
     105         zx = glamu(ji,jj) * 1.e3 
     106         zy = gphiu(ji,jj) * 1.e3 
     107         DO jk=1, jpk 
     108            zdu = 0.5_wp * (pdept(ji  ,jj,jk) + pdept(ji+1,jj,jk)) 
     109            IF (zdu < zH) THEN 
     110               zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH)) 
     111               pu(ji,jj,jk) = (za * zf * zy * EXP(-(zx**2+zy**2)/zlambda**2)) * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk) 
     112            ELSE 
     113               pu(ji,jj,jk) = 0._wp 
     114            ENDIF 
    119115         END DO 
    120       END DO 
     116      END_2D 
    121117      ! 
    122       DO ji=1, jpi 
    123          DO jj=1, jpjm1 
    124             zx = glamv(ji,jj) * 1.e3 
    125             zy = gphiv(ji,jj) * 1.e3 
    126             DO jk=1, jpk 
    127                zdv = 0.5_wp * (pdept(ji  ,jj,jk) + pdept(ji,jj+1,jk)) 
    128                IF (zdv < zH) THEN 
    129                   zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH)) 
    130                   pv(ji,jj,jk) = -(za * zf * zx * EXP(-(zx**2+zy**2)/zlambda**2)) * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk) 
    131                ELSE 
    132                   pv(ji,jj,jk) = 0._wp 
    133                ENDIF 
    134             END DO 
     118      DO_2D_00_00 
     119         zx = glamv(ji,jj) * 1.e3 
     120         zy = gphiv(ji,jj) * 1.e3 
     121         DO jk=1, jpk 
     122            zdv = 0.5_wp * (pdept(ji  ,jj,jk) + pdept(ji,jj+1,jk)) 
     123            IF (zdv < zH) THEN 
     124               zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH)) 
     125               pv(ji,jj,jk) = -(za * zf * zx * EXP(-(zx**2+zy**2)/zlambda**2)) * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk) 
     126            ELSE 
     127               pv(ji,jj,jk) = 0._wp 
     128            ENDIF 
    135129         END DO 
    136       END DO 
    137  
    138       CALL lbc_lnk( 'usrdef_istate', pu, 'U', -1. ) 
    139       CALL lbc_lnk( 'usrdef_istate', pv, 'V', -1. ) 
     130      END_2D 
     131      ! 
     132      CALL lbc_lnk_multi( 'usrdef_istate', pu, 'U', -1., pv, 'V', -1. ) 
    140133      !    
    141134   END SUBROUTINE usr_def_istate 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/VORTEX/MY_SRC/usrdef_zgr.F90

    r12377 r13197  
    192192      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. )           ! set surrounding land to zero (here jperio=0 ==>> closed) 
    193193      ! 
    194       k_bot(:,:) = INT( z2d(:,:) )           ! =jpkm1 over the ocean point, =0 elsewhere 
     194      k_bot(:,:) = NINT( z2d(:,:) )          ! =jpkm1 over the ocean point, =0 elsewhere 
    195195      ! 
    196196      k_top(:,:) = MIN( 1 , k_bot(:,:) )     ! = 1    over the ocean point, =0 elsewhere 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/WAD/MY_SRC/usrdef_hgr.F90

    r10074 r13197  
    2626   PUBLIC   usr_def_hgr   ! called by domhgr.F90 
    2727 
     28   !! * Substitutions 
     29#  include "do_loop_substitute.h90" 
    2830   !!---------------------------------------------------------------------- 
    2931   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    7274      !                       !==  grid point position  ==!   (in kilometers) 
    7375      zfact = rn_dx * 1.e-3         ! conversion in km 
    74       DO jj = 1, jpj 
    75          DO ji = 1, jpi             ! longitude 
    76             plamt(ji,jj) = zfact * (  - 0.5 + REAL( ji-1 + nimpp-1 , wp )  )   
    77             plamu(ji,jj) = zfact * (          REAL( ji-1 + nimpp-1 , wp )  ) 
    78             plamv(ji,jj) = plamt(ji,jj) 
    79             plamf(ji,jj) = plamu(ji,jj) 
    80             !                       ! latitude 
    81             pphit(ji,jj) = zfact * (  - 0.5 + REAL( jj-1 + njmpp-1 , wp )  ) 
    82             pphiu(ji,jj) = pphit(ji,jj) 
    83             pphiv(ji,jj) = zfact * (          REAL( jj-1 + njmpp-1 , wp )  ) 
    84             pphif(ji,jj) = pphiv(ji,jj) 
    85          END DO 
    86       END DO 
     76      DO_2D_11_11 
     77         !                       ! longitude 
     78         plamt(ji,jj) = zfact * (  - 0.5 + REAL( ji-1 + nimpp-1 , wp )  )   
     79         plamu(ji,jj) = zfact * (          REAL( ji-1 + nimpp-1 , wp )  ) 
     80         plamv(ji,jj) = plamt(ji,jj) 
     81         plamf(ji,jj) = plamu(ji,jj) 
     82         !                       ! latitude 
     83         pphit(ji,jj) = zfact * (  - 0.5 + REAL( jj-1 + njmpp-1 , wp )  ) 
     84         pphiu(ji,jj) = pphit(ji,jj) 
     85         pphiv(ji,jj) = zfact * (          REAL( jj-1 + njmpp-1 , wp )  ) 
     86         pphif(ji,jj) = pphiv(ji,jj) 
     87      END_2D 
    8788      ! 
    8889      !                       !==  Horizontal scale factors  ==!   (in meters)  
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/WAD/MY_SRC/usrdef_istate.F90

    r10074 r13197  
    2626   PUBLIC   usr_def_istate   ! called in istate.F90 
    2727 
     28   !! * Substitutions 
     29#  include "do_loop_substitute.h90" 
    2830   !!---------------------------------------------------------------------- 
    2931   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    174176      ! Apply minimum wetdepth criterion 
    175177      ! 
    176       do jj = 1,jpj 
    177          do ji = 1,jpi 
    178             IF( ht_0(ji,jj) + pssh(ji,jj) < rn_wdmin1 ) THEN 
    179                pssh(ji,jj) = ptmask(ji,jj,1)*( rn_wdmin1 - ht_0(ji,jj) ) 
    180             ENDIF 
    181          end do 
    182       end do 
     178      DO_2D_11_11 
     179         IF( ht_0(ji,jj) + pssh(ji,jj) < rn_wdmin1 ) THEN 
     180            pssh(ji,jj) = ptmask(ji,jj,1)*( rn_wdmin1 - ht_0(ji,jj) ) 
     181         ENDIF 
     182      END_2D 
    183183      ! 
    184184   END SUBROUTINE usr_def_istate 
  • NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/WAD/MY_SRC/usrdef_zgr.F90

    r12377 r13197  
    2929   PUBLIC   usr_def_zgr        ! called by domzgr.F90 
    3030 
     31   !! * Substitutions 
     32#  include "do_loop_substitute.h90" 
    3133   !!---------------------------------------------------------------------- 
    3234   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    242244      ! at v-point: averaging zht 
    243245      zhv = 0._wp 
    244       DO jj = 1, jpjm1 
    245          zhv(:,jj) = 0.5_wp * ( zht(:,jj) + zht(:,jj+1) ) 
    246       END DO 
     246      DO_2D_00_00 
     247         zhv(ji,jj) = 0.5_wp * ( zht(ji,jj) + zht(ji,jj+1) ) 
     248      END_2D 
    247249      CALL lbc_lnk( 'usrdef_zgr', zhv, 'V', 1. )     ! boundary condition: this mask the surrounding grid-points 
    248250      DO jj = mj0(1), mj1(1)   ! first  row of global domain only 
     
    279281         ht_0 = zht 
    280282         k_bot(:,:) = jpkm1 * k_top(:,:)  !* bottom ocean = jpk-1 (here use k_top as a land mask) 
    281          DO jj = 1, jpj 
    282             DO ji = 1, jpi 
    283               IF( zht(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN 
    284                 k_bot(ji,jj) = 0 
    285                 k_top(ji,jj) = 0 
    286               ENDIF 
    287            END DO 
    288          END DO 
     283         DO_2D_11_11 
     284            IF( zht(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN 
     285               k_bot(ji,jj) = 0 
     286               k_top(ji,jj) = 0 
     287            ENDIF 
     288         END_2D 
    289289         ! 
    290290         !                                !* terrain-following coordinate with e3.(k)=cst) 
    291291         !                                !  OVERFLOW case : identical with j-index (T=V, U=F) 
    292          DO jj = 1, jpjm1 
    293             DO ji = 1, jpim1 
    294               z1_jpkm1 = 1._wp / REAL( k_bot(ji,jj) - k_top(ji,jj) + 1 , wp) 
    295               DO jk = 1, jpk 
    296                   zwet = MAX( zht(ji,jj), rn_wdmin1 ) 
    297                   pdept(ji,jj,jk) = zwet * z1_jpkm1 * ( REAL( jk   , wp ) - 0.5_wp ) 
    298                   pdepw(ji,jj,jk) = zwet * z1_jpkm1 * ( REAL( jk-1 , wp )          ) 
    299                   pe3t (ji,jj,jk) = zwet * z1_jpkm1 
    300                   pe3w (ji,jj,jk) = zwet * z1_jpkm1 
    301                   zwet = MAX( zhu(ji,jj), rn_wdmin1 ) 
    302                   pe3u (ji,jj,jk) = zwet * z1_jpkm1 
    303                   pe3uw(ji,jj,jk) = zwet * z1_jpkm1 
    304                   pe3f (ji,jj,jk) = zwet * z1_jpkm1 
    305                   zwet = MAX( zhv(ji,jj), rn_wdmin1 ) 
    306                   pe3v (ji,jj,jk) = zwet * z1_jpkm1 
    307                   pe3vw(ji,jj,jk) = zwet * z1_jpkm1 
    308               END DO       
    309            END DO       
    310          END DO       
     292         DO_2D_00_00 
     293            z1_jpkm1 = 1._wp / REAL( k_bot(ji,jj) - k_top(ji,jj) + 1 , wp) 
     294            DO jk = 1, jpk 
     295               zwet = MAX( zht(ji,jj), rn_wdmin1 ) 
     296               pdept(ji,jj,jk) = zwet * z1_jpkm1 * ( REAL( jk   , wp ) - 0.5_wp ) 
     297               pdepw(ji,jj,jk) = zwet * z1_jpkm1 * ( REAL( jk-1 , wp )          ) 
     298               pe3t (ji,jj,jk) = zwet * z1_jpkm1 
     299               pe3w (ji,jj,jk) = zwet * z1_jpkm1 
     300               zwet = MAX( zhu(ji,jj), rn_wdmin1 ) 
     301               pe3u (ji,jj,jk) = zwet * z1_jpkm1 
     302               pe3uw(ji,jj,jk) = zwet * z1_jpkm1 
     303               pe3f (ji,jj,jk) = zwet * z1_jpkm1 
     304               zwet = MAX( zhv(ji,jj), rn_wdmin1 ) 
     305               pe3v (ji,jj,jk) = zwet * z1_jpkm1 
     306               pe3vw(ji,jj,jk) = zwet * z1_jpkm1 
     307            END DO 
     308         END_2D      
    311309         CALL lbc_lnk( 'usrdef_zgr', pdept, 'T', 1. ) 
    312310         CALL lbc_lnk( 'usrdef_zgr', pdepw, 'T', 1. ) 
Note: See TracChangeset for help on using the changeset viewer.