New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 13463 for NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/TRA/tramle.F90 – NEMO

Ignore:
Timestamp:
2020-09-14T17:40:34+02:00 (4 years ago)
Author:
andmirek
Message:

Ticket #2195:update to trunk 13461

Location:
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11351_fldread_with_XIOS

    • Property svn:externals
      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
         5^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
         8 
         9# SETTE 
         10^/utils/CI/sette@13382        sette 
  • NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/TRA/tramle.F90

    r10425 r13463  
    4141 
    4242   REAL(wp) ::   r5_21 = 5.e0 / 21.e0   ! factor used in mle streamfunction computation 
    43    REAL(wp) ::   rb_c                   ! ML buoyancy criteria = g rho_c /rau0 where rho_c is defined in zdfmld 
     43   REAL(wp) ::   rb_c                   ! ML buoyancy criteria = g rho_c /rho0 where rho_c is defined in zdfmld 
    4444   REAL(wp) ::   rc_f                   ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_mle=1 case 
    4545 
     
    4848 
    4949   !! * Substitutions 
    50 #  include "vectopt_loop_substitute.h90" 
     50#  include "do_loop_substitute.h90" 
     51#  include "domzgr_substitute.h90" 
    5152   !!---------------------------------------------------------------------- 
    5253   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    5657CONTAINS 
    5758 
    58    SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype ) 
     59   SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype, Kmm ) 
    5960      !!---------------------------------------------------------------------- 
    6061      !!                  ***  ROUTINE tra_mle_trp  *** 
     
    7172      !!                p.n = p.n + z._mle 
    7273      !! 
    73       !! ** Action  : - (pun,pvn,pwn) increased by the mle transport 
     74      !! ** Action  : - (pu,pv,pw) increased by the mle transport 
    7475      !!                CAUTION, the transport is not updated at the last line/raw 
    7576      !!                         this may be a problem for some advection schemes 
     
    8081      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index 
    8182      INTEGER                         , INTENT(in   ) ::   kit000     ! first time step index 
     83      INTEGER                         , INTENT(in   ) ::   Kmm        ! ocean time level index 
    8284      CHARACTER(len=3)                , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
    8385      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu         ! in : 3 ocean transport components 
     
    98100      inml_mle(:,:) = mbkt(:,:) + 1                    ! init. to number of ocean w-level (T-level + 1) 
    99101      IF ( nla10 > 0 ) THEN                            ! avoid case where first level is thicker than 10m 
    100          DO jk = jpkm1, nlb10, -1                      ! from the bottom to nlb10 (10m) 
    101             DO jj = 1, jpj 
    102                DO ji = 1, jpi                          ! index of the w-level at the ML based 
    103                   IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer 
    104                END DO 
    105             END DO 
    106          END DO 
     102         DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 ) 
     103            IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer 
     104         END_3D 
    107105      ENDIF 
    108106      ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 )                  ! max level of the computation 
     
    112110      zbm (:,:) = 0._wp 
    113111      zn2 (:,:) = 0._wp 
    114       DO jk = 1, ikmax                                 ! MLD and mean buoyancy and N2 over the mixed layer 
    115          DO jj = 1, jpj 
    116             DO ji = 1, jpi 
    117                zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points 
    118                zmld(ji,jj) = zmld(ji,jj) + zc 
    119                zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0 
    120                zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 
    121             END DO 
    122          END DO 
    123       END DO 
     112      DO_3D( 1, 1, 1, 1, 1, ikmax ) 
     113         zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points 
     114         zmld(ji,jj) = zmld(ji,jj) + zc 
     115         zbm (ji,jj) = zbm (ji,jj) + zc * (rho0 - rhop(ji,jj,jk) ) * r1_rho0 
     116         zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 
     117      END_3D 
    124118 
    125119      SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts 
    126120      CASE ( 0 )                                               != min of the 2 neighbour MLDs 
    127          DO jj = 1, jpjm1 
    128             DO ji = 1, fs_jpim1   ! vector opt. 
    129                zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) ) 
    130                zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) ) 
    131             END DO 
    132          END DO 
     121         DO_2D( 1, 0, 1, 0 ) 
     122            zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) ) 
     123            zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) ) 
     124         END_2D 
    133125      CASE ( 1 )                                               != average of the 2 neighbour MLDs 
    134          DO jj = 1, jpjm1 
    135             DO ji = 1, fs_jpim1   ! vector opt. 
    136                zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp 
    137                zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp 
    138             END DO 
    139          END DO 
     126         DO_2D( 1, 0, 1, 0 ) 
     127            zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp 
     128            zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp 
     129         END_2D 
    140130      CASE ( 2 )                                               != max of the 2 neighbour MLDs 
    141          DO jj = 1, jpjm1 
    142             DO ji = 1, fs_jpim1   ! vector opt. 
    143                zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) ) 
    144                zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) ) 
    145             END DO 
    146          END DO 
     131         DO_2D( 1, 0, 1, 0 ) 
     132            zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) ) 
     133            zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) ) 
     134         END_2D 
    147135      END SELECT 
    148136      !                                                ! convert density into buoyancy 
    149       zbm(:,:) = + grav * zbm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 
     137      zbm(:,:) = + grav * zbm(:,:) / MAX( e3t(:,:,1,Kmm), zmld(:,:) ) 
    150138      ! 
    151139      ! 
     
    158146      ! 
    159147      IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation 
    160          DO jj = 1, jpjm1 
    161             DO ji = 1, fs_jpim1   ! vector opt. 
    162                zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2_e1u(ji,jj)                                            & 
    163                   &           * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )   & 
    164                   &           / (  MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   ) 
    165                   ! 
    166                zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1_e2v(ji,jj)                                            & 
    167                   &           * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )   & 
    168                   &           / (  MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   ) 
    169             END DO 
    170          END DO 
     148         DO_2D( 1, 0, 1, 0 ) 
     149            zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2_e1u(ji,jj)                                            & 
     150               &           * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )   & 
     151               &           / (  MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   ) 
     152               ! 
     153            zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1_e2v(ji,jj)                                            & 
     154               &           * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )   & 
     155               &           / (  MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   ) 
     156         END_2D 
    171157         ! 
    172158      ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 
    173          DO jj = 1, jpjm1 
    174             DO ji = 1, fs_jpim1   ! vector opt. 
    175                zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2_e1u(ji,jj)               & 
    176                   &                  * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) 
    177                   ! 
    178                zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1_e2v(ji,jj)               & 
    179                   &                  * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) 
    180             END DO 
    181          END DO 
     159         DO_2D( 1, 0, 1, 0 ) 
     160            zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2_e1u(ji,jj)               & 
     161               &                  * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) 
     162               ! 
     163            zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1_e2v(ji,jj)               & 
     164               &                  * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) 
     165         END_2D 
    182166      ENDIF 
    183167      ! 
    184168      IF( nn_conv == 1 ) THEN              ! No MLE in case of convection 
    185          DO jj = 1, jpjm1 
    186             DO ji = 1, fs_jpim1   ! vector opt. 
    187                IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp )   zpsim_u(ji,jj) = 0._wp 
    188                IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp )   zpsim_v(ji,jj) = 0._wp 
    189             END DO 
    190          END DO 
     169         DO_2D( 1, 0, 1, 0 ) 
     170            IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp )   zpsim_u(ji,jj) = 0._wp 
     171            IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp )   zpsim_v(ji,jj) = 0._wp 
     172         END_2D 
    191173      ENDIF 
    192174      ! 
    193175      !                                      !==  structure function value at uw- and vw-points  ==! 
    194       DO jj = 1, jpjm1 
    195          DO ji = 1, fs_jpim1   ! vector opt. 
    196             zhu(ji,jj) = 1._wp / zhu(ji,jj)                   ! hu --> 1/hu 
    197             zhv(ji,jj) = 1._wp / zhv(ji,jj) 
    198          END DO 
    199       END DO 
     176      DO_2D( 1, 0, 1, 0 ) 
     177         zhu(ji,jj) = 1._wp / zhu(ji,jj)                   ! hu --> 1/hu 
     178         zhv(ji,jj) = 1._wp / zhv(ji,jj) 
     179      END_2D 
    200180      ! 
    201181      zpsi_uw(:,:,:) = 0._wp 
    202182      zpsi_vw(:,:,:) = 0._wp 
    203183      ! 
    204       DO jk = 2, ikmax                                ! start from 2 : surface value = 0 
    205          DO jj = 1, jpjm1 
    206             DO ji = 1, fs_jpim1   ! vector opt. 
    207                zcuw = 1._wp - ( gdepw_n(ji+1,jj,jk) + gdepw_n(ji,jj,jk) ) * zhu(ji,jj) 
    208                zcvw = 1._wp - ( gdepw_n(ji,jj+1,jk) + gdepw_n(ji,jj,jk) ) * zhv(ji,jj) 
    209                zcuw = zcuw * zcuw 
    210                zcvw = zcvw * zcvw 
    211                zmuw = MAX(  0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw )  ) 
    212                zmvw = MAX(  0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw )  ) 
    213                ! 
    214                zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk) 
    215                zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk) 
    216             END DO 
    217          END DO 
    218       END DO 
     184      DO_3D( 1, 0, 1, 0, 2, ikmax ) 
     185         zcuw = 1._wp - ( gdepw(ji+1,jj,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhu(ji,jj) 
     186         zcvw = 1._wp - ( gdepw(ji,jj+1,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhv(ji,jj) 
     187         zcuw = zcuw * zcuw 
     188         zcvw = zcvw * zcvw 
     189         zmuw = MAX(  0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw )  ) 
     190         zmvw = MAX(  0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw )  ) 
     191         ! 
     192         zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk) 
     193         zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk) 
     194      END_3D 
    219195      ! 
    220196      !                                      !==  transport increased by the MLE induced transport ==! 
    221197      DO jk = 1, ikmax 
    222          DO jj = 1, jpjm1                          ! CAUTION pu,pv must be defined at row/column i=1 / j=1 
    223             DO ji = 1, fs_jpim1   ! vector opt. 
    224                pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 
    225                pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 
    226             END DO 
    227          END DO 
    228          DO jj = 2, jpjm1 
    229             DO ji = fs_2, fs_jpim1   ! vector opt. 
    230                pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk)   & 
    231                   &                          + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) ) 
    232             END DO 
    233          END DO 
     198         DO_2D( 1, 0, 1, 0 ) 
     199            pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 
     200            pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 
     201         END_2D 
     202         DO_2D( 0, 0, 0, 0 ) 
     203            pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk)   & 
     204               &                          + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) ) 
     205         END_2D 
    234206      END DO 
    235207 
     
    266238      !!---------------------------------------------------------------------- 
    267239 
    268       REWIND( numnam_ref )              ! Namelist namtra_mle in reference namelist : Tracer advection scheme 
    269240      READ  ( numnam_ref, namtra_mle, IOSTAT = ios, ERR = 901) 
    270 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_mle in reference namelist', lwp ) 
    271  
    272       REWIND( numnam_cfg )              ! Namelist namtra_mle in configuration namelist : Tracer advection scheme 
     241901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_mle in reference namelist' ) 
     242 
    273243      READ  ( numnam_cfg, namtra_mle, IOSTAT = ios, ERR = 902 ) 
    274 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_mle in configuration namelist', lwp ) 
     244902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_mle in configuration namelist' ) 
    275245      IF(lwm) WRITE ( numond, namtra_mle ) 
    276246 
     
    304274      IF( ln_mle ) THEN                ! MLE initialisation 
    305275         ! 
    306          rb_c = grav * rn_rho_c_mle /rau0        ! Mixed Layer buoyancy criteria 
     276         rb_c = grav * rn_rho_c_mle /rho0        ! Mixed Layer buoyancy criteria 
    307277         IF(lwp) WRITE(numout,*) 
    308278         IF(lwp) WRITE(numout,*) '      ML buoyancy criteria = ', rb_c, ' m/s2 ' 
     
    313283            IF( ierr /= 0 )   CALL ctl_stop( 'tra_adv_mle_init: failed to allocate arrays' ) 
    314284            z1_t2 = 1._wp / ( rn_time * rn_time ) 
    315             DO jj = 2, jpj                           ! "coriolis+ time^-1" at u- & v-points 
    316                DO ji = fs_2, jpi   ! vector opt. 
    317                   zfu = ( ff_f(ji,jj) + ff_f(ji,jj-1) ) * 0.5_wp 
    318                   zfv = ( ff_f(ji,jj) + ff_f(ji-1,jj) ) * 0.5_wp 
    319                   rfu(ji,jj) = SQRT(  zfu * zfu + z1_t2 ) 
    320                   rfv(ji,jj) = SQRT(  zfv * zfv + z1_t2 ) 
    321                END DO 
    322             END DO 
    323             CALL lbc_lnk_multi( 'tramle', rfu, 'U', 1. , rfv, 'V', 1. ) 
     285            DO_2D( 0, 1, 0, 1 ) 
     286               zfu = ( ff_f(ji,jj) + ff_f(ji,jj-1) ) * 0.5_wp 
     287               zfv = ( ff_f(ji,jj) + ff_f(ji-1,jj) ) * 0.5_wp 
     288               rfu(ji,jj) = SQRT(  zfu * zfu + z1_t2 ) 
     289               rfv(ji,jj) = SQRT(  zfv * zfv + z1_t2 ) 
     290            END_2D 
     291            CALL lbc_lnk_multi( 'tramle', rfu, 'U', 1.0_wp , rfv, 'V', 1.0_wp ) 
    324292            ! 
    325293         ELSEIF( nn_mle == 1 ) THEN           ! MLE array allocation & initialisation 
Note: See TracChangeset for help on using the changeset viewer.