Changeset 12217


Ignore:
Timestamp:
2019-12-12T17:01:26+01:00 (5 months ago)
Author:
agn
Message:

make FKOSM changes to tramle.F90

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRA/tramle.F90

    r12178 r12217  
    2121   USE lbclnk         ! lateral boundary condition / mpp link 
    2222 
     23   USE zdfosm, ONLY  : ln_osm_mle, hmle, dbdx_mle, dbdy_mle, mld_prof 
     24 
    2325   IMPLICIT NONE 
    2426   PRIVATE 
     
    5658CONTAINS 
    5759 
    58    SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype ) 
    59       !!---------------------------------------------------------------------- 
    60       !!                  ***  ROUTINE tra_mle_trp  *** 
    61       !! 
    62       !! ** Purpose :   Add to the transport the Mixed Layer Eddy induced transport 
    63       !! 
    64       !! ** Method  :   The 3 components of the Mixed Layer Eddy (MLE) induced 
    65       !!              transport are computed as follows : 
    66       !!                zu_mle = dk[ zpsi_uw ] 
    67       !!                zv_mle = dk[ zpsi_vw ] 
    68       !!                zw_mle = - di[ zpsi_uw ] - dj[ zpsi_vw ] 
    69       !!                where zpsi is the MLE streamfunction at uw and vw points (see the doc) 
    70       !!              and added to the input velocity : 
    71       !!                p.n = p.n + z._mle 
    72       !! 
    73       !! ** Action  : - (pun,pvn,pwn) increased by the mle transport 
    74       !!                CAUTION, the transport is not updated at the last line/raw 
    75       !!                         this may be a problem for some advection schemes 
    76       !! 
    77       !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 
    78       !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 
    79       !!---------------------------------------------------------------------- 
    80       INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index 
    81       INTEGER                         , INTENT(in   ) ::   kit000     ! first time step index 
    82       CHARACTER(len=3)                , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
    83       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu         ! in : 3 ocean transport components 
    84       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pv         ! out: same 3  transport components 
    85       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pw         !   increased by the MLE induced transport 
    86       ! 
    87       INTEGER  ::   ji, jj, jk          ! dummy loop indices 
    88       INTEGER  ::   ii, ij, ik, ikmax   ! local integers 
    89       REAL(wp) ::   zcuw, zmuw, zc      ! local scalar 
    90       REAL(wp) ::   zcvw, zmvw          !   -      - 
    91       INTEGER , DIMENSION(jpi,jpj)     :: inml_mle 
    92       REAL(wp), DIMENSION(jpi,jpj)     :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH 
    93       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw 
    94       !!---------------------------------------------------------------------- 
    95       ! 
    96       !                                      !==  MLD used for MLE  ==! 
    97       !                                                ! compute from the 10m density to deal with the diurnal cycle 
    98       inml_mle(:,:) = mbkt(:,:) + 1                    ! init. to number of ocean w-level (T-level + 1) 
    99       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 
    107       ENDIF 
    108       ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 )                  ! max level of the computation 
    109       ! 
    110       ! 
    111       zmld(:,:) = 0._wp                      !==   Horizontal shape of the MLE  ==! 
    112       zbm (:,:) = 0._wp 
    113       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 
    124  
    125       SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts 
    126       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 
    133       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 
    140       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 
    147       END SELECT 
    148       !                                                ! convert density into buoyancy 
    149       zbm(:,:) = + grav * zbm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 
    150       ! 
    151       ! 
    152       !                                      !==  Magnitude of the MLE stream function  ==! 
    153       ! 
    154       !                 di[bm]  Ds 
    155       ! Psi = Ce  H^2 ---------------- e2u  mu(z)   where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 ) 
    156       !                  e1u   Lf fu                      and the e2u for the "transport" 
    157       !                                                      (not *e3u as divided by e3u at the end) 
    158       ! 
    159       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 
    171          ! 
    172       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 
    182       ENDIF 
    183       ! 
    184       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 
    191       ENDIF 
    192       ! 
    193       !                                      !==  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 
    200       ! 
    201       zpsi_uw(:,:,:) = 0._wp 
    202       zpsi_vw(:,:,:) = 0._wp 
    203       ! 
    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 
    219       ! 
    220       !                                      !==  transport increased by the MLE induced transport ==! 
    221       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)   & 
     60  SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype ) 
     61    !!---------------------------------------------------------------------- 
     62    !!                  ***  ROUTINE tra_mle_trp  *** 
     63    !! 
     64    !! ** Purpose :   Add to the transport the Mixed Layer Eddy induced transport 
     65    !! 
     66    !! ** Method  :   The 3 components of the Mixed Layer Eddy (MLE) induced 
     67    !!              transport are computed as follows : 
     68    !!                zu_mle = dk[ zpsi_uw ] 
     69    !!                zv_mle = dk[ zpsi_vw ] 
     70    !!                zw_mle = - di[ zpsi_uw ] - dj[ zpsi_vw ] 
     71    !!                where zpsi is the MLE streamfunction at uw and vw points (see the doc) 
     72    !!              and added to the input velocity : 
     73    !!                p.n = p.n + z._mle 
     74    !! 
     75    !! ** Action  : - (pun,pvn,pwn) increased by the mle transport 
     76    !!                CAUTION, the transport is not updated at the last line/raw 
     77    !!                         this may be a problem for some advection schemes 
     78    !! 
     79    !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008 
     80    !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008 
     81    !!---------------------------------------------------------------------- 
     82    INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index 
     83    INTEGER                         , INTENT(in   ) ::   kit000     ! first time step index 
     84    CHARACTER(len=3)                , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
     85    REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu         ! in : 3 ocean transport components 
     86    REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pv         ! out: same 3  transport components 
     87    REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pw         !   increased by the MLE induced transport 
     88    ! 
     89    INTEGER  ::   ji, jj, jk          ! dummy loop indices 
     90    INTEGER  ::   ii, ij, ik, ikmax   ! local integers 
     91    REAL(wp) ::   zcuw, zmuw, zc      ! local scalar 
     92    REAL(wp) ::   zcvw, zmvw          !   -      - 
     93    INTEGER , DIMENSION(jpi,jpj)     :: inml_mle 
     94    REAL(wp), DIMENSION(jpi,jpj)     :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH 
     95    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw 
     96    !!---------------------------------------------------------------------- 
     97    ! 
     98    ! 
     99    IF(ln_osm_mle) THEN 
     100       ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 )                  ! max level of the computation 
     101       ! 
     102       ! 
     103       SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts 
     104       CASE ( 0 )                                               != min of the 2 neighbour MLDs 
     105          DO jj = 1, jpjm1 
     106             DO ji = 1, fs_jpim1   ! vector opt. 
     107                zhu(ji,jj) = MIN( hmle(ji+1,jj), hmle(ji,jj) ) 
     108                zhv(ji,jj) = MIN( hmle(ji,jj+1), hmle(ji,jj) ) 
     109             END DO 
     110          END DO 
     111       CASE ( 1 )                                               != average of the 2 neighbour MLDs 
     112          DO jj = 1, jpjm1 
     113             DO ji = 1, fs_jpim1   ! vector opt. 
     114                zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 
     115                zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) 
     116             END DO 
     117          END DO 
     118       CASE ( 2 )                                               != max of the 2 neighbour MLDs 
     119          DO jj = 1, jpjm1 
     120             DO ji = 1, fs_jpim1   ! vector opt. 
     121                zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 
     122                zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) 
     123             END DO 
     124          END DO 
     125       END SELECT 
     126       IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation 
     127          DO jj = 1, jpjm1 
     128             DO ji = 1, fs_jpim1   ! vector opt. 
     129                zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2u(ji,jj)                                            & 
     130                     &           * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) )   & 
     131                     &           / (  MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   ) 
     132                ! 
     133                zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1v(ji,jj)                                            & 
     134                     &           * dbdy_mle(ji,jj)  * MIN( 111.e3_wp , e2v(ji,jj) )   & 
     135                     &           / (  MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   ) 
     136             END DO 
     137          END DO 
     138          ! 
     139       ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 
     140          DO jj = 1, jpjm1 
     141             DO ji = 1, fs_jpim1   ! vector opt. 
     142                zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2u(ji,jj)               & 
     143                     &                  * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) 
     144                ! 
     145                zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1v(ji,jj)               & 
     146                     &                  * dbdy_mle(ji,jj) * MIN( 111.e3_wp , e2v(ji,jj) ) 
     147             END DO 
     148          END DO 
     149       ENDIF 
     150 
     151    ELSE !do not use osn_mle 
     152       !                                      !==  MLD used for MLE  ==! 
     153       !                                                ! convert density into buoyancy 
     154       inml_mle(:,:) = mbkt(:,:) + 1                    ! init. to number of ocean w-level (T-level + 1) 
     155       IF ( nla10 > 0 ) THEN                            ! avoid case where first level is thicker than 10m 
     156          DO jk = jpkm1, nlb10, -1                      ! from the bottom to nlb10 (10m) 
     157             DO jj = 1, jpj 
     158                DO ji = 1, jpi                          ! index of the w-level at the ML based 
     159                   IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer 
     160                END DO 
     161             END DO 
     162          END DO 
     163       ENDIF 
     164 
     165       ! 
     166       ! 
     167       zmld(:,:) = 0._wp                      !==   Horizontal shape of the MLE  ==! 
     168       zbm (:,:) = 0._wp 
     169       zn2 (:,:) = 0._wp 
     170       DO jk = 1, ikmax                                 ! MLD and mean buoyancy and N2 over the mixed layer 
     171          DO jj = 1, jpj 
     172             DO ji = 1, jpi 
     173                zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points 
     174                zmld(ji,jj) = zmld(ji,jj) + zc 
     175                zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0 
     176                zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 
     177             END DO 
     178          END DO 
     179       END DO 
     180 
     181       SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts 
     182       CASE ( 0 )                                               != min of the 2 neighbour MLDs 
     183          DO jj = 1, jpjm1 
     184             DO ji = 1, fs_jpim1   ! vector opt. 
     185                zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) ) 
     186                zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) ) 
     187             END DO 
     188          END DO 
     189       CASE ( 1 )                                               != average of the 2 neighbour MLDs 
     190          DO jj = 1, jpjm1 
     191             DO ji = 1, fs_jpim1   ! vector opt. 
     192                zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp 
     193                zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp 
     194             END DO 
     195          END DO 
     196       CASE ( 2 )                                               != max of the 2 neighbour MLDs 
     197          DO jj = 1, jpjm1 
     198             DO ji = 1, fs_jpim1   ! vector opt. 
     199                zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) ) 
     200                zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) ) 
     201             END DO 
     202          END DO 
     203       END SELECT 
     204       !                                                ! convert density into buoyancy 
     205       ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 )                  ! max level of the computation 
     206       zbm(:,:) = + grav * zbm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) ) 
     207       ! 
     208       ! 
     209       !                                      !==  Magnitude of the MLE stream function  ==! 
     210       ! 
     211       !                 di[bm]  Ds 
     212       ! Psi = Ce  H^2 ---------------- e2u  mu(z)   where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 ) 
     213       !                  e1u   Lf fu                      and the e2u for the "transport" 
     214       !                                                      (not *e3u as divided by e3u at the end) 
     215       ! 
     216       IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation 
     217          DO jj = 1, jpjm1 
     218             DO ji = 1, fs_jpim1   ! vector opt. 
     219                zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2_e1u(ji,jj)                                            & 
     220                     &           * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )   & 
     221                     &           / (  MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   ) 
     222                ! 
     223                zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1_e2v(ji,jj)                                            & 
     224                     &           * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )   & 
     225                     &           / (  MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   ) 
     226             END DO 
     227          END DO 
     228          ! 
     229       ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 
     230          DO jj = 1, jpjm1 
     231             DO ji = 1, fs_jpim1   ! vector opt. 
     232                zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2_e1u(ji,jj)               & 
     233                     &                  * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) 
     234                ! 
     235                zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1_e2v(ji,jj)               & 
     236                     &                  * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) ) 
     237             END DO 
     238          END DO 
     239       ENDIF 
     240       ! 
     241       IF( nn_conv == 1 ) THEN              ! No MLE in case of convection 
     242          DO jj = 1, jpjm1 
     243             DO ji = 1, fs_jpim1   ! vector opt. 
     244                IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp )   zpsim_u(ji,jj) = 0._wp 
     245                IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp )   zpsim_v(ji,jj) = 0._wp 
     246             END DO 
     247          END DO 
     248       ENDIF 
     249       ! 
     250    ENDIF  ! end of lm_osm_mle loop 
     251    !                                      !==  structure function value at uw- and vw-points  ==! 
     252    DO jj = 1, jpjm1 
     253       DO ji = 1, fs_jpim1   ! vector opt. 
     254          zhu(ji,jj) = 1._wp / zhu(ji,jj)                   ! hu --> 1/hu 
     255          zhv(ji,jj) = 1._wp / zhv(ji,jj) 
     256       END DO 
     257    END DO 
     258    ! 
     259    zpsi_uw(:,:,:) = 0._wp 
     260    zpsi_vw(:,:,:) = 0._wp 
     261    ! 
     262    DO jk = 2, ikmax                                ! start from 2 : surface value = 0 
     263       DO jj = 1, jpjm1 
     264          DO ji = 1, fs_jpim1   ! vector opt. 
     265             zcuw = 1._wp - ( gdepw_n(ji+1,jj,jk) + gdepw_n(ji,jj,jk) ) * zhu(ji,jj) 
     266             zcvw = 1._wp - ( gdepw_n(ji,jj+1,jk) + gdepw_n(ji,jj,jk) ) * zhv(ji,jj) 
     267             zcuw = zcuw * zcuw 
     268             zcvw = zcvw * zcvw 
     269             zmuw = MAX(  0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw )  ) 
     270             zmvw = MAX(  0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw )  ) 
     271             ! 
     272             zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk) 
     273             zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk) 
     274          END DO 
     275       END DO 
     276    END DO 
     277    ! 
     278    !                                      !==  transport increased by the MLE induced transport ==! 
     279    DO jk = 1, ikmax 
     280       DO jj = 1, jpjm1                          ! CAUTION pu,pv must be defined at row/column i=1 / j=1 
     281          DO ji = 1, fs_jpim1   ! vector opt. 
     282             pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 
     283             pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 
     284          END DO 
     285       END DO 
     286       DO jj = 2, jpjm1 
     287          DO ji = fs_2, fs_jpim1   ! vector opt. 
     288             pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk)   & 
    231289                  &                          + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) ) 
    232             END DO 
    233          END DO 
    234       END DO 
    235  
    236       IF( cdtype == 'TRA') THEN              !==  outputs  ==! 
    237          ! 
    238          zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:)      ! Lf = N H / f 
    239          CALL iom_put( "Lf_NHpf" , zLf_NH  )    ! Lf = N H / f 
    240          ! 
    241          ! divide by cross distance to give streamfunction with dimensions m^2/s 
    242          DO jk = 1, ikmax+1 
    243             zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:) 
    244             zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:) 
    245          END DO 
    246          CALL iom_put( "psiu_mle", zpsi_uw )    ! i-mle streamfunction 
    247          CALL iom_put( "psiv_mle", zpsi_vw )    ! j-mle streamfunction 
    248       ENDIF 
    249       ! 
    250    END SUBROUTINE tra_mle_trp 
     290          END DO 
     291       END DO 
     292    END DO 
     293 
     294    IF( cdtype == 'TRA') THEN              !==  outputs  ==! 
     295       ! 
     296       IF (ln_osm_mle) THEN 
     297          zLf_NH(:,:) = SQRT( rb_c * hmle(:,:) ) * r1_ft(:,:)      ! Lf = N H / f 
     298       ELSE 
     299          zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:)      ! Lf = N H / f 
     300       END IF 
     301       CALL iom_put( "Lf_NHpf" , zLf_NH  )    ! Lf = N H / f 
     302       ! 
     303       ! divide by cross distance to give streamfunction with dimensions m^2/s 
     304       DO jk = 1, ikmax+1 
     305          zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:) 
     306          zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:) 
     307       END DO 
     308       CALL iom_put( "psiu_mle", zpsi_uw )    ! i-mle streamfunction 
     309       CALL iom_put( "psiv_mle", zpsi_vw )    ! j-mle streamfunction 
     310    ENDIF 
     311    ! 
     312  END SUBROUTINE tra_mle_trp 
    251313 
    252314 
     
    290352         WRITE(numout,*) '         Density difference used to define ML for FK              rn_rho_c_mle  = ', rn_rho_c_mle 
    291353      ENDIF 
    292       ! 
    293       IF(lwp) THEN 
     354 
     355 
     356      IF( ln_osm_mle .AND. ln_mle ) THEN 
     357         WRITE(numout,*) 'WARNING: You are running with both OSM-FK and default FK' 
     358         ! CALL ctl_stop('STOP in zdf_osm_init: Cannot run with both OSM-FK and default FK') 
     359     END IF 
     360 
     361     IF(lwp) THEN 
    294362         WRITE(numout,*) 
    295363         IF( ln_mle ) THEN 
Note: See TracChangeset for help on using the changeset viewer.