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 5260 for branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90 – NEMO

Ignore:
Timestamp:
2015-05-12T12:37:15+02:00 (9 years ago)
Author:
deazer
Message:

Merged branch with Trunk at revision 5253.
Checked with SETTE, passes modified iodef.xml for AMM12 experiment

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90

    r4313 r5260  
    22   !!============================================================================== 
    33   !!                       ***  MODULE  tranpc  *** 
    4    !! Ocean active tracers:  non penetrative convection scheme 
     4   !! Ocean active tracers:  non penetrative convective adjustment scheme 
    55   !!============================================================================== 
    66   !! History :  1.0  ! 1990-09  (G. Madec)  Original code 
     
    99   !!            3.0  ! 2008-06  (G. Madec)  applied on ta, sa and called before tranxt in step.F90 
    1010   !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA 
     11   !!            3.7  ! 2014-06  (L. Brodeau) new algorithm based on local Brunt-Vaisala freq. 
    1112   !!---------------------------------------------------------------------- 
    1213 
     
    1415   !!   tra_npc : apply the non penetrative convection scheme 
    1516   !!---------------------------------------------------------------------- 
    16    USE oce             ! ocean dynamics and active tracers  
     17   USE oce             ! ocean dynamics and active tracers 
    1718   USE dom_oce         ! ocean space and time domain 
     19   USE phycst          ! physical constants 
    1820   USE zdf_oce         ! ocean vertical physics 
    19    USE trdmod_oce      ! ocean active tracer trends 
     21   USE trd_oce         ! ocean active tracer trends 
    2022   USE trdtra          ! ocean active tracer trends 
    21    USE eosbn2          ! equation of state (eos routine)  
     23   USE eosbn2          ! equation of state (eos routine) 
     24   ! 
    2225   USE lbclnk          ! lateral boundary conditions (or mpp link) 
    2326   USE in_out_manager  ! I/O manager 
     
    2932   PRIVATE 
    3033 
    31    PUBLIC   tra_npc       ! routine called by step.F90 
     34   PUBLIC   tra_npc    ! routine called by step.F90 
    3235 
    3336   !! * Substitutions 
    3437#  include "domzgr_substitute.h90" 
    35    !!---------------------------------------------------------------------- 
    36    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    37    !! $Id$  
     38#  include "vectopt_loop_substitute.h90" 
     39   !!---------------------------------------------------------------------- 
     40   !! NEMO/OPA 3.6 , NEMO Consortium (2014) 
     41   !! $Id$ 
    3842   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    3943   !!---------------------------------------------------------------------- 
     
    4448      !!                  ***  ROUTINE tranpc  *** 
    4549      !! 
    46       !! ** Purpose :   Non penetrative convective adjustment scheme. solve  
     50      !! ** Purpose : Non-penetrative convective adjustment scheme. solve 
    4751      !!      the static instability of the water column on after fields 
    4852      !!      while conserving heat and salt contents. 
    4953      !! 
    50       !! ** Method  :   The algorithm used converges in a maximium of jpk  
    51       !!      iterations. instabilities are treated when the vertical density 
    52       !!      gradient is less than 1.e-5. 
    53       !!      l_trdtra=T: the trend associated with this algorithm is saved. 
     54      !! ** Method  : updated algorithm able to deal with non-linear equation of state 
     55      !!              (i.e. static stability computed locally) 
    5456      !! 
    5557      !! ** Action  : - (ta,sa) after the application od the npc scheme 
    56       !!              - save the associated trends (ttrd,strd) ('key_trdtra') 
     58      !!              - send the associated trends for on-line diagnostics (l_trdtra=T) 
    5759      !! 
    58       !! References : Madec, et al., 1991, JPO, 21, 9, 1349-1371. 
     60      !! References :     Madec, et al., 1991, JPO, 21, 9, 1349-1371. 
    5961      !!---------------------------------------------------------------------- 
    60       ! 
    6162      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    6263      ! 
    6364      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    6465      INTEGER  ::   inpcc        ! number of statically instable water column 
    65       INTEGER  ::   inpci        ! number of iteration for npc scheme 
    66       INTEGER  ::   jiter, jkdown, jkp        ! ??? 
    67       INTEGER  ::   ikbot, ik, ikup, ikdown   ! ??? 
    68       REAL(wp) ::   ze3tot, zta, zsa, zraua, ze3dwn 
    69       REAL(wp), POINTER, DIMENSION(:,:  ) :: zwx, zwy, zwz 
    70       REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt, ztrds, zrhop 
     66      INTEGER  ::   jiter, ikbot, ik, ikup, ikdown, ilayer, ikm   ! local integers 
     67      LOGICAL  ::   l_bottom_reached, l_column_treated 
     68      REAL(wp) ::   zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z 
     69      REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_r2dt 
     70      REAL(wp), POINTER, DIMENSION(:)       ::   zvn2   ! vertical profile of N2 at 1 given point... 
     71      REAL(wp), POINTER, DIMENSION(:,:)     ::   zvts   ! vertical profile of T and S at 1 given point... 
     72      REAL(wp), POINTER, DIMENSION(:,:)     ::   zvab   ! vertical profile of alpha and beta 
     73      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zn2    ! N^2  
     74      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zab    ! alpha and beta 
     75      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   ztrdt, ztrds   ! 3D workspace 
     76      ! 
     77      !!LB debug: 
     78      LOGICAL, PARAMETER :: l_LB_debug = .FALSE. 
     79      INTEGER :: ilc1, jlc1, klc1, nncpu 
     80      LOGICAL :: lp_monitor_point = .FALSE. 
     81      !!LB debug. 
    7182      !!---------------------------------------------------------------------- 
    7283      ! 
    7384      IF( nn_timing == 1 )  CALL timing_start('tra_npc') 
    7485      ! 
    75       CALL wrk_alloc(jpi, jpj, jpk, zrhop ) 
    76       CALL wrk_alloc(jpi, jpk, zwx, zwy, zwz ) 
    77       ! 
    7886      IF( MOD( kt, nn_npc ) == 0 ) THEN 
    79  
    80          inpcc = 0 
    81          inpci = 0 
    82  
    83          CALL eos( tsa, rhd, zrhop, fsdept_n(:,:,:) )         ! Potential density 
    84  
    85          IF( l_trdtra )   THEN                    !* Save ta and sa trends 
     87         ! 
     88         CALL wrk_alloc( jpi, jpj, jpk, zn2 )    ! N2 
     89         CALL wrk_alloc( jpi, jpj, jpk, 2, zab ) ! Alpha and Beta 
     90         CALL wrk_alloc( jpk, 2, zvts, zvab )    ! 1D column vector at point ji,jj 
     91         CALL wrk_alloc( jpk, zvn2 )             ! 1D column vector at point ji,jj 
     92 
     93         IF( l_trdtra )   THEN                    !* Save initial after fields 
    8694            CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
    8795            ztrdt(:,:,:) = tsa(:,:,:,jp_tem)  
     
    8997         ENDIF 
    9098 
    91          !                                                ! =============== 
    92          DO jj = 1, jpj                                   !  Vertical slab 
    93             !                                             ! =============== 
    94             !  Static instability pointer  
    95             ! ---------------------------- 
    96             DO jk = 1, jpkm1 
    97                DO ji = 1, jpi 
    98                   zwx(ji,jk) = ( zrhop(ji,jj,jk) - zrhop(ji,jj,jk+1) ) * tmask(ji,jj,jk+1) 
    99                END DO 
    100             END DO 
    101  
    102             ! 1.1 do not consider the boundary points 
    103  
    104             ! even if east-west cyclic b. c. do not considere ji=1 or jpi 
    105             DO jk = 1, jpkm1 
    106                zwx( 1 ,jk) = 0.e0 
    107                zwx(jpi,jk) = 0.e0 
    108             END DO 
    109             ! even if south-symmetric b. c. used, do not considere jj=1 
    110             IF( jj == 1 )   zwx(:,:) = 0.e0 
    111  
    112             DO jk = 1, jpkm1 
    113                DO ji = 1, jpi 
    114                   zwx(ji,jk) = 1. 
    115                   IF( zwx(ji,jk) < 1.e-5 ) zwx(ji,jk) = 0.e0 
    116                END DO 
    117             END DO 
    118  
    119             zwy(:,1) = 0.e0 
    120             DO ji = 1, jpi 
    121                DO jk = 1, jpkm1 
    122                   zwy(ji,1) = zwy(ji,1) + zwx(ji,jk) 
    123                END DO 
    124             END DO 
    125  
    126             zwz(1,1) = 0.e0 
    127             DO ji = 1, jpi 
    128                zwz(1,1) = zwz(1,1) + zwy(ji,1) 
    129             END DO 
    130  
    131             inpcc = inpcc + NINT( zwz(1,1) ) 
    132  
    133  
    134             ! 2. Vertical mixing for each instable portion of the density profil 
    135             ! ------------------------------------------------------------------ 
    136  
    137             IF( zwz(1,1) /= 0.e0 ) THEN         ! -->> the density profil is statically instable : 
    138                DO ji = 1, jpi 
    139                   IF( zwy(ji,1) /= 0.e0 ) THEN 
     99         !LB debug: 
     100         IF( lwp .AND. l_LB_debug ) THEN 
     101            WRITE(numout,*) 
     102            WRITE(numout,*) 'LOLO: entering tra_npc, kt, narea =', kt, narea 
     103         ENDIF 
     104         !LBdebug: Monitoring of 1 column subject to convection... 
     105         IF( l_LB_debug ) THEN 
     106            ! Location of 1 known convection spot to follow what's happening in the water column 
     107            ilc1 = 54 ;  jlc1 = 15 ; ! Labrador ORCA1 4x4 cpus: 
     108            nncpu = 15  ; ! the CPU domain contains the convection spot 
     109            !ilc1 = 14 ;  jlc1 = 13 ; ! Labrador ORCA1 8x8 cpus: 
     110            !nncpu = 54  ; ! the CPU domain contains the convection spot 
     111            klc1 =  mbkt(ilc1,jlc1) ! bottom of the ocean for debug point...           
     112         ENDIF 
     113         !LBdebug. 
     114 
     115         CALL eos_rab( tsa, zab )         ! after alpha and beta 
     116         CALL bn2    ( tsa, zab, zn2 )    ! after Brunt-Vaisala 
     117         
     118         inpcc = 0 
     119 
     120         DO jj = 2, jpjm1                 ! interior column only 
     121            DO ji = fs_2, fs_jpim1 
     122               ! 
     123               IF( tmask(ji,jj,2) == 1 ) THEN      ! At least 2 ocean points 
     124                  !                                     ! consider one ocean column  
     125                  zvts(:,jp_tem) = tsa(ji,jj,:,jp_tem)      ! temperature 
     126                  zvts(:,jp_sal) = tsa(ji,jj,:,jp_sal)      ! salinity 
     127 
     128                  zvab(:,jp_tem)  = zab(ji,jj,:,jp_tem)     ! Alpha  
     129                  zvab(:,jp_sal)  = zab(ji,jj,:,jp_sal)     ! Beta   
     130                  zvn2(:)         = zn2(ji,jj,:)            ! N^2  
     131                  
     132                  IF( l_LB_debug ) THEN                  !LB debug: 
     133                     lp_monitor_point = .FALSE. 
     134                     IF( ( ji == ilc1 ).AND.( jj == jlc1 ) ) lp_monitor_point = .TRUE. 
     135                     ! writing only if on CPU domain where conv region is: 
     136                     lp_monitor_point = (narea == nncpu).AND.lp_monitor_point  
     137                      
     138                     IF(lp_monitor_point) THEN 
     139                        WRITE(numout,*) '' ;WRITE(numout,*) '' ; 
     140                       WRITE(numout,'("Time step = ",i6.6," !!!")')  kt 
     141                        WRITE(numout,'(" *** BEFORE anything, N^2 for point ",i3,",",i3,":" )') ji,jj 
     142                        DO jk = 1, klc1 
     143                           WRITE(numout,*) jk, zvn2(jk) 
     144                        END DO 
     145                        WRITE(numout,*) ' ' 
     146                     ENDIF 
     147                  ENDIF                                  !LB debug  end 
     148 
     149                  ikbot = mbkt(ji,jj)   ! ikbot: ocean bottom T-level 
     150                  ik = 1                ! because N2 is irrelevant at the surface level (will start at ik=2) 
     151                  ilayer = 0 
     152                  jiter  = 0 
     153                  l_column_treated = .FALSE. 
     154                  
     155                  DO WHILE ( .NOT. l_column_treated ) 
    140156                     ! 
    141                      ikbot = mbkt(ji,jj)        ! ikbot: ocean bottom T-level 
     157                     jiter = jiter + 1 
     158                     
     159                     IF( jiter >= 400 ) EXIT 
     160                     
     161                     l_bottom_reached = .FALSE. 
     162 
     163                     DO WHILE ( .NOT. l_bottom_reached ) 
     164 
     165                        ik = ik + 1 
     166                        
     167                        !! Checking level ik for instability 
     168                        !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
     169 
     170                        IF( zvn2(ik) < 0. ) THEN ! Instability found! 
     171 
     172                           ikm  = ik              ! first level whith negative N2 
     173                           ilayer = ilayer + 1    ! yet another layer found.... 
     174                           IF(jiter == 1) inpcc = inpcc + 1 
     175 
     176                           IF(l_LB_debug .AND. lp_monitor_point) & 
     177                              & WRITE(numout,*) 'Negative N2 at ik =', ikm, ' layer nb.', ilayer, & 
     178                              & ' inpcc =', inpcc 
     179 
     180                           !! Case we mix with upper regions where N2==0: 
     181                           !! All the points above ikup where N2 == 0 must also be mixed => we go 
     182                           !! upward to find a new ikup, where the layer doesn't have N2==0 
     183                           ikup = ikm 
     184                           DO jk = ikm, 2, -1 
     185                              ikup = ikup - 1 
     186                              IF( (zvn2(jk-1) > 0.).OR.(ikup == 1) ) EXIT 
     187                           END DO 
     188                           
     189                           ! adjusting ikup if the upper part of the unstable column was neutral (N2=0) 
     190                           IF((zvn2(ikup+1) == 0.).AND.(ikup /= 1)) ikup = ikup+1 ; 
     191 
     192                           
     193                           IF( lp_monitor_point )   WRITE(numout,*) ' => ikup is =', ikup, ' layer nb.', ilayer 
     194                           
     195                           zsum_temp = 0._wp 
     196                           zsum_sali = 0._wp 
     197                           zsum_alfa = 0._wp 
     198                           zsum_beta = 0._wp 
     199                           zsum_z    = 0._wp 
     200                                                     
     201                           DO jk = ikup, ikbot+1      ! Inside the instable (and overlying neutral) portion of the column 
     202                              ! 
     203                              IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '     -> summing for jk =', jk 
     204                              ! 
     205                              zdz       = fse3t(ji,jj,jk) 
     206                              zsum_temp = zsum_temp + zvts(jk,jp_tem)*zdz 
     207                              zsum_sali = zsum_sali + zvts(jk,jp_sal)*zdz 
     208                              zsum_alfa = zsum_alfa + zvab(jk,jp_tem)*zdz 
     209                              zsum_beta = zsum_beta + zvab(jk,jp_sal)*zdz 
     210                              zsum_z    = zsum_z    + zdz 
     211                              ! 
     212                              !! EXIT if we found the bottom of the unstable portion of the water column     
     213                              IF( (zvn2(jk+1) > 0.).OR.(jk == ikbot ).OR.((jk==ikm).AND.(zvn2(jk+1) == 0.)) )   EXIT 
     214                           END DO 
     215                           
     216                           !ik     = jk !LB remove? 
     217                           ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative N2 
     218                           
     219                           IF(l_LB_debug .AND. lp_monitor_point) & 
     220                              &    WRITE(numout,*) '  => ikdown =', ikdown, '  layer nb.', ilayer 
     221                           
     222                           ! Mixing Temperature and salinity between ikup and ikdown: 
     223                           zta   = zsum_temp/zsum_z 
     224                           zsa   = zsum_sali/zsum_z 
     225                           zalfa = zsum_alfa/zsum_z 
     226                           zbeta = zsum_beta/zsum_z 
     227 
     228                           IF(l_LB_debug .AND. lp_monitor_point) THEN 
     229                              WRITE(numout,*) '  => Mean temp. in that portion =', zta 
     230                              WRITE(numout,*) '  => Mean sali. in that portion =', zsa 
     231                              WRITE(numout,*) '  => Mean Alpha in that portion =', zalfa 
     232                              WRITE(numout,*) '  => Mean Beta  in that portion =', zbeta 
     233                           ENDIF 
     234 
     235                           !! Homogenaizing the temperature, salinity, alpha and beta in this portion of the column 
     236                           DO jk = ikup, ikdown 
     237                              zvts(jk,jp_tem) = zta 
     238                              zvts(jk,jp_sal) = zsa 
     239                              zvab(jk,jp_tem) = zalfa 
     240                              zvab(jk,jp_sal) = zbeta 
     241                           END DO 
     242                           ! 
     243                           !! Before updating N2, it is possible that another unstable 
     244                           !! layer exists underneath the one we just homogeneized! 
     245                           ik = ikdown 
     246                           !  
     247                        ENDIF  ! IF( zvn2(ik+1) < 0. ) THEN 
     248                        ! 
     249                        IF( ik == ikbot ) l_bottom_reached = .TRUE. 
     250                        ! 
     251                     END DO ! DO WHILE ( .NOT. l_bottom_reached ) 
     252 
     253                     IF( ik /= ikbot )   STOP 'ERROR: tranpc.F90 => PROBLEM #1' 
     254                     
     255                     ! ******* At this stage ik == ikbot ! ******* 
     256                     
     257                     IF( ilayer > 0 ) THEN 
     258                        !! least an unstable layer has been found 
     259                        !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion 
     260                        !! => Need to re-compute N2! will use Alpha and Beta! 
     261                        ! 
     262                        DO jk = ikup+1, ikdown+1   ! we must go 1 point deeper than ikdown!      
     263                           !! Doing exactly as in eosbn2.F90: 
     264                           !! * Except that we only are interested in the sign of N2 !!! 
     265                           !!   => just considering the vertical gradient of density 
     266                           zrw =   (fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk)) & 
     267                               & / (fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk)) 
     268                           zaw = zvab(jk,jp_tem) * (1._wp - zrw) + zvab(jk-1,jp_tem) * zrw 
     269                           zbw = zvab(jk,jp_sal) * (1._wp - zrw) + zvab(jk-1,jp_sal) * zrw 
     270                           
     271                           !zvn2(jk) = grav*( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     & 
     272                           !     &           - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )  & 
     273                           !     &       / fse3w(ji,jj,jk) * tmask(ji,jj,jk) 
     274                           zvn2(jk) = (  zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     & 
     275                                &      - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )   
     276                        END DO 
     277 
     278                        IF(l_LB_debug .AND. lp_monitor_point) THEN 
     279                           WRITE(numout, '(" *** After iteration #",i3.3,", N^2 for point ",i3,",",i3,":" )') & 
     280                              & jiter, ji,jj 
     281                           DO jk = 1, klc1 
     282                              WRITE(numout,*) jk, zvn2(jk) 
     283                           END DO 
     284                           WRITE(numout,*) ' ' 
     285                        ENDIF 
     286 
     287                        ik     = 1  ! starting again at the surface for the next iteration 
     288                        ilayer = 0 
     289                     ENDIF 
    142290                     ! 
    143                      DO jiter = 1, jpk          ! vertical iteration 
    144                         ! 
    145                         ! search of ikup : the first static instability from the sea surface 
    146                         ! 
    147                         ik = 0 
    148 220                     CONTINUE 
    149                         ik = ik + 1 
    150                         IF( ik >= ikbot ) GO TO 200 
    151                         zwx(ji,ik) = zrhop(ji,jj,ik) - zrhop(ji,jj,ik+1) 
    152                         IF( zwx(ji,ik) <= 0.e0 ) GO TO 220 
    153                         ikup = ik 
    154                         ! the density profil is instable below ikup 
    155                         ! ikdown : bottom of the instable portion of the density profil 
    156                         ! search of ikdown and vertical mixing from ikup to ikdown 
    157                         ! 
    158                         ze3tot= fse3t(ji,jj,ikup) 
    159                         zta   = tsa  (ji,jj,ikup,jp_tem) 
    160                         zsa   = tsa  (ji,jj,ikup,jp_sal) 
    161                         zraua = zrhop(ji,jj,ikup) 
    162                         ! 
    163                         DO jkdown = ikup+1, ikbot-1 
    164                            IF( zraua <= zrhop(ji,jj,jkdown) ) THEN 
    165                               ikdown = jkdown 
    166                               GO TO 240 
    167                            ENDIF 
    168                            ze3dwn =  fse3t(ji,jj,jkdown) 
    169                            ze3tot =  ze3tot + ze3dwn 
    170                            zta   = ( zta*(ze3tot-ze3dwn) + tsa(ji,jj,jkdown,jp_tem)*ze3dwn )/ze3tot 
    171                            zsa   = ( zsa*(ze3tot-ze3dwn) + tsa(ji,jj,jkdown,jp_sal)*ze3dwn )/ze3tot 
    172                            zraua = ( zraua*(ze3tot-ze3dwn) + zrhop(ji,jj,jkdown)*ze3dwn )/ze3tot 
    173                            inpci = inpci+1 
    174                         END DO 
    175                         ikdown = ikbot-1 
    176 240                     CONTINUE 
    177                         ! 
    178                         DO jkp = ikup, ikdown-1 
    179                            tsa  (ji,jj,jkp,jp_tem) = zta 
    180                            tsa  (ji,jj,jkp,jp_sal) = zsa 
    181                            zrhop(ji,jj,jkp       ) = zraua 
    182                         END DO 
    183                         IF (ikdown == ikbot-1 .AND. zraua >= zrhop(ji,jj,ikdown) ) THEN 
    184                            tsa  (ji,jj,jkp,jp_tem) = zta 
    185                            tsa  (ji,jj,jkp,jp_sal) = zsa 
    186                            zrhop(ji,jj,ikdown    ) = zraua 
    187                         ENDIF 
    188                      END DO 
    189                   ENDIF 
    190 200               CONTINUE 
    191                END DO 
    192                ! <<-- no more static instability on slab jj 
    193             ENDIF 
    194             !                                             ! =============== 
    195          END DO                                           !   End of slab 
    196          !                                                ! =============== 
    197          !  
    198          IF( l_trdtra )   THEN         ! save the Non penetrative mixing trends for diagnostic 
    199             ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    200             ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
    201             CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_npc, ztrdt ) 
    202             CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_npc, ztrds ) 
     291                     IF( ik >= ikbot ) THEN 
     292                        IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '    --- exiting jiter loop ---' 
     293                        l_column_treated = .TRUE. 
     294                     ENDIF 
     295                     ! 
     296                  END DO ! DO WHILE ( .NOT. l_column_treated ) 
     297 
     298                  !! Updating tsa: 
     299                  tsa(ji,jj,:,jp_tem) = zvts(:,jp_tem) 
     300                  tsa(ji,jj,:,jp_sal) = zvts(:,jp_sal) 
     301 
     302                  !! lolo:  Should we update something else???? 
     303                  !! => like alpha and beta? 
     304 
     305                  IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '' 
     306 
     307               ENDIF ! IF( tmask(ji,jj,3) == 1 ) THEN 
     308 
     309            END DO ! ji 
     310         END DO ! jj 
     311         ! 
     312         IF( l_trdtra ) THEN         ! send the Non penetrative mixing trends for diagnostic 
     313            z1_r2dt = 1._wp / (2._wp * rdt) 
     314            ztrdt(:,:,:) = ( tsa(:,:,:,jp_tem) - ztrdt(:,:,:) ) * z1_r2dt 
     315            ztrds(:,:,:) = ( tsa(:,:,:,jp_sal) - ztrds(:,:,:) ) * z1_r2dt 
     316            CALL trd_tra( kt, 'TRA', jp_tem, jptra_npc, ztrdt ) 
     317            CALL trd_tra( kt, 'TRA', jp_sal, jptra_npc, ztrds ) 
    203318            CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
    204319         ENDIF 
    205        
    206          ! Lateral boundary conditions on ( ta, sa )   ( Unchanged sign) 
    207          ! ------------------------------============ 
     320         ! 
    208321         CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ;   CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) 
    209        
    210  
    211          !  2. non penetrative convective scheme statistics 
    212          !  ----------------------------------------------- 
    213          IF( nn_npcp /= 0 .AND. MOD( kt, nn_npcp ) == 0 ) THEN 
    214             IF(lwp) WRITE(numout,*)' kt=',kt, ' number of statically instable',   & 
    215                &                   ' water column : ',inpcc, ' number of iteration : ',inpci 
    216          ENDIF 
    217          ! 
    218       ENDIF 
    219       ! 
    220       CALL wrk_dealloc(jpi, jpj, jpk, zrhop ) 
    221       CALL wrk_dealloc(jpi, jpk, zwx, zwy, zwz ) 
     322         ! 
     323         IF(lwp) THEN 
     324            WRITE(numout,*) 'LOLO: exiting tra_npc, kt =', kt 
     325            WRITE(numout,*)' => number of statically instable water column : ',inpcc 
     326            WRITE(numout,*) '' ; WRITE(numout,*) '' 
     327         ENDIF 
     328         ! 
     329         CALL wrk_dealloc(jpi, jpj, jpk, zn2 ) 
     330         CALL wrk_dealloc(jpi, jpj, jpk, 2, zab ) 
     331         CALL wrk_dealloc(jpk, zvn2 ) 
     332         CALL wrk_dealloc(jpk, 2, zvts, zvab ) 
     333         ! 
     334      ENDIF   ! IF( MOD( kt, nn_npc ) == 0 ) THEN 
    222335      ! 
    223336      IF( nn_timing == 1 )  CALL timing_stop('tra_npc') 
Note: See TracChangeset for help on using the changeset viewer.