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

Ignore:
Timestamp:
2015-07-10T13:28:53+02:00 (9 years ago)
Author:
timgraham
Message:

Merged head of trunk into branch

File:
1 edited

Legend:

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

    r4313 r5581  
    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.6  ! 2015-05  (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, ikp, ikup, ikdown, ilayer, ik_low   ! 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), PARAMETER :: zn2_zero = 1.e-14_wp       ! acceptance criteria for neutrality (N2==0) 
     71      REAL(wp), POINTER, DIMENSION(:)       ::   zvn2   ! vertical profile of N2 at 1 given point... 
     72      REAL(wp), POINTER, DIMENSION(:,:)     ::   zvts   ! vertical profile of T and S at 1 given point... 
     73      REAL(wp), POINTER, DIMENSION(:,:)     ::   zvab   ! vertical profile of alpha and beta 
     74      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zn2    ! N^2  
     75      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zab    ! alpha and beta 
     76      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   ztrdt, ztrds   ! 3D workspace 
     77      ! 
     78      LOGICAL, PARAMETER :: l_LB_debug = .FALSE. ! set to true if you want to follow what is 
     79      INTEGER :: ilc1, jlc1, klc1, nncpu         ! actually happening in a water column at point "ilc1, jlc1" 
     80      LOGICAL :: lp_monitor_point = .FALSE.      ! in CPU domain "nncpu" 
    7181      !!---------------------------------------------------------------------- 
    7282      ! 
    7383      IF( nn_timing == 1 )  CALL timing_start('tra_npc') 
    7484      ! 
    75       CALL wrk_alloc(jpi, jpj, jpk, zrhop ) 
    76       CALL wrk_alloc(jpi, jpk, zwx, zwy, zwz ) 
    77       ! 
    7885      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 
     86         ! 
     87         CALL wrk_alloc( jpi, jpj, jpk, zn2 )    ! N2 
     88         CALL wrk_alloc( jpi, jpj, jpk, 2, zab ) ! Alpha and Beta 
     89         CALL wrk_alloc( jpk, 2, zvts, zvab )    ! 1D column vector at point ji,jj 
     90         CALL wrk_alloc( jpk, zvn2 )             ! 1D column vector at point ji,jj 
     91 
     92         IF( l_trdtra )   THEN                    !* Save initial after fields 
    8693            CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
    8794            ztrdt(:,:,:) = tsa(:,:,:,jp_tem)  
     
    8996         ENDIF 
    9097 
    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 
     98         IF( l_LB_debug ) THEN 
     99            ! Location of 1 known convection site to follow what's happening in the water column 
     100            ilc1 = 45 ;  jlc1 = 3 ; !  ORCA2 4x4, Antarctic coast, more than 2 unstable portions in the  water column...            
     101            nncpu = 1  ;            ! the CPU domain contains the convection spot 
     102            klc1 =  mbkt(ilc1,jlc1) ! bottom of the ocean for debug point...           
     103         ENDIF 
     104          
     105         CALL eos_rab( tsa, zab )         ! after alpha and beta (given on T-points) 
     106         CALL bn2    ( tsa, zab, zn2 )    ! after Brunt-Vaisala  (given on W-points) 
     107         
     108         inpcc = 0 
     109 
     110         DO jj = 2, jpjm1                 ! interior column only 
     111            DO ji = fs_2, fs_jpim1 
     112               ! 
     113               IF( tmask(ji,jj,2) == 1 ) THEN      ! At least 2 ocean points 
     114                  !                                     ! consider one ocean column  
     115                  zvts(:,jp_tem) = tsa(ji,jj,:,jp_tem)      ! temperature 
     116                  zvts(:,jp_sal) = tsa(ji,jj,:,jp_sal)      ! salinity 
     117 
     118                  zvab(:,jp_tem)  = zab(ji,jj,:,jp_tem)     ! Alpha  
     119                  zvab(:,jp_sal)  = zab(ji,jj,:,jp_sal)     ! Beta   
     120                  zvn2(:)         = zn2(ji,jj,:)            ! N^2  
     121                  
     122                  IF( l_LB_debug ) THEN                  !LB debug: 
     123                     lp_monitor_point = .FALSE. 
     124                     IF( ( ji == ilc1 ).AND.( jj == jlc1 ) ) lp_monitor_point = .TRUE. 
     125                     ! writing only if on CPU domain where conv region is: 
     126                     lp_monitor_point = (narea == nncpu).AND.lp_monitor_point                       
     127                  ENDIF                                  !LB debug  end 
     128 
     129                  ikbot = mbkt(ji,jj)   ! ikbot: ocean bottom T-level 
     130                  ikp = 1                  ! because N2 is irrelevant at the surface level (will start at ikp=2) 
     131                  ilayer = 0 
     132                  jiter  = 0 
     133                  l_column_treated = .FALSE. 
     134                  
     135                  DO WHILE ( .NOT. l_column_treated ) 
    140136                     ! 
    141                      ikbot = mbkt(ji,jj)        ! ikbot: ocean bottom T-level 
     137                     jiter = jiter + 1 
     138                     
     139                     IF( jiter >= 400 ) EXIT 
     140                     
     141                     l_bottom_reached = .FALSE. 
     142 
     143                     DO WHILE ( .NOT. l_bottom_reached ) 
     144 
     145                        ikp = ikp + 1 
     146                        
     147                        !! Testing level ikp for instability 
     148                        !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
     149                        IF( zvn2(ikp) <  -zn2_zero ) THEN ! Instability found! 
     150 
     151                           ilayer = ilayer + 1    ! yet another instable portion of the water column found.... 
     152 
     153                           IF( lp_monitor_point ) THEN  
     154                              WRITE(numout,*) 
     155                              IF( ilayer == 1 .AND. jiter == 1 ) THEN   ! first time a column is spoted with an instability 
     156                                 WRITE(numout,*) 
     157                                 WRITE(numout,*) 'Time step = ',kt,' !!!' 
     158                              ENDIF 
     159                              WRITE(numout,*)  ' * Iteration #',jiter,': found instable portion #',ilayer,   & 
     160                                 &                                    ' in column! Starting at ikp =', ikp 
     161                              WRITE(numout,*)  ' *** N2 for point (i,j) = ',ji,' , ',jj 
     162                              DO jk = 1, klc1 
     163                                 WRITE(numout,*) jk, zvn2(jk) 
     164                              END DO 
     165                              WRITE(numout,*) 
     166                           ENDIF 
     167                            
     168 
     169                           IF( jiter == 1 )   inpcc = inpcc + 1  
     170 
     171                           IF( lp_monitor_point )   WRITE(numout, *) 'Negative N2 at ikp =',ikp,' for layer #', ilayer 
     172 
     173                           !! ikup is the uppermost point where mixing will start: 
     174                           ikup = ikp - 1 ! ikup is always "at most at ikp-1", less if neutral levels overlying 
     175                            
     176                           !! If the points above ikp-1 have N2 == 0 they must also be mixed: 
     177                           IF( ikp > 2 ) THEN 
     178                              DO jk = ikp-1, 2, -1 
     179                                 IF( ABS(zvn2(jk)) < zn2_zero ) THEN 
     180                                    ikup = ikup - 1  ! 1 more upper level has N2=0 and must be added for the mixing 
     181                                 ELSE 
     182                                    EXIT 
     183                                 ENDIF 
     184                              END DO 
     185                           ENDIF 
     186                            
     187                           IF( ikup < 1 )   CALL ctl_stop( 'tra_npc :  PROBLEM #1') 
     188 
     189                           zsum_temp = 0._wp 
     190                           zsum_sali = 0._wp 
     191                           zsum_alfa = 0._wp 
     192                           zsum_beta = 0._wp 
     193                           zsum_z    = 0._wp 
     194                                                     
     195                           DO jk = ikup, ikbot      ! Inside the instable (and overlying neutral) portion of the column 
     196                              ! 
     197                              zdz       = fse3t(ji,jj,jk) 
     198                              zsum_temp = zsum_temp + zvts(jk,jp_tem)*zdz 
     199                              zsum_sali = zsum_sali + zvts(jk,jp_sal)*zdz 
     200                              zsum_alfa = zsum_alfa + zvab(jk,jp_tem)*zdz 
     201                              zsum_beta = zsum_beta + zvab(jk,jp_sal)*zdz 
     202                              zsum_z    = zsum_z    + zdz 
     203                              !                               
     204                              IF( jk == ikbot ) EXIT ! avoid array-index overshoot in case ikbot = jpk, cause we're calling jk+1 next line 
     205                              !! EXIT when we have reached the last layer that is instable (N2<0) or neutral (N2=0): 
     206                              IF( zvn2(jk+1) > zn2_zero ) EXIT 
     207                           END DO 
     208                           
     209                           ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative or neutral N2 
     210                           IF( ikup == ikdown )   CALL ctl_stop( 'tra_npc :  PROBLEM #2') 
     211 
     212                           ! Mixing Temperature, salinity, alpha and beta from ikup to ikdown included: 
     213                           zta   = zsum_temp/zsum_z 
     214                           zsa   = zsum_sali/zsum_z 
     215                           zalfa = zsum_alfa/zsum_z 
     216                           zbeta = zsum_beta/zsum_z 
     217 
     218                           IF( lp_monitor_point ) THEN 
     219                              WRITE(numout,*) 'MIXED T, S, alfa and beta between ikup =',ikup,   & 
     220                                 &            ' and ikdown =',ikdown,', in layer #',ilayer 
     221                              WRITE(numout,*) '  => Mean temp. in that portion =', zta 
     222                              WRITE(numout,*) '  => Mean sali. in that portion =', zsa 
     223                              WRITE(numout,*) '  => Mean Alfa  in that portion =', zalfa 
     224                              WRITE(numout,*) '  => Mean Beta  in that portion =', zbeta 
     225                           ENDIF 
     226 
     227                           !! Homogenaizing the temperature, salinity, alpha and beta in this portion of the column 
     228                           DO jk = ikup, ikdown 
     229                              zvts(jk,jp_tem) = zta 
     230                              zvts(jk,jp_sal) = zsa 
     231                              zvab(jk,jp_tem) = zalfa 
     232                              zvab(jk,jp_sal) = zbeta 
     233                           END DO 
     234                            
     235                            
     236                           !! Updating N2 in the relvant portion of the water column 
     237                           !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion 
     238                           !! => Need to re-compute N2! will use Alpha and Beta! 
     239                            
     240                           ikup   = MAX(2,ikup)         ! ikup can never be 1 ! 
     241                           ik_low = MIN(ikdown+1,ikbot) ! we must go 1 point deeper than ikdown! 
     242                            
     243                           DO jk = ikup, ik_low              ! we must go 1 point deeper than ikdown! 
     244 
     245                              !! Interpolating alfa and beta at W point: 
     246                              zrw =  (fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk)) & 
     247                                 & / (fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk)) 
     248                              zaw = zvab(jk,jp_tem) * (1._wp - zrw) + zvab(jk-1,jp_tem) * zrw 
     249                              zbw = zvab(jk,jp_sal) * (1._wp - zrw) + zvab(jk-1,jp_sal) * zrw 
     250 
     251                              !! N2 at W point, doing exactly as in eosbn2.F90: 
     252                              zvn2(jk) = grav*( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     & 
     253                                 &            - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )  & 
     254                                 &       / fse3w(ji,jj,jk) * tmask(ji,jj,jk) 
     255 
     256                              !! OR, faster  => just considering the vertical gradient of density 
     257                              !! as only the signa maters... 
     258                              !zvn2(jk) = (  zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     & 
     259                              !     &      - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  ) 
     260 
     261                           END DO 
     262                         
     263                           ikp = MIN(ikdown+1,ikbot) 
     264                            
     265 
     266                        ENDIF  !IF( zvn2(ikp) < 0. ) 
     267 
     268 
     269                        IF( ikp == ikbot ) l_bottom_reached = .TRUE. 
     270                        ! 
     271                     END DO ! DO WHILE ( .NOT. l_bottom_reached ) 
     272 
     273                     IF( ikp /= ikbot )   CALL ctl_stop( 'tra_npc :  PROBLEM #3') 
     274                     
     275                     ! ******* At this stage ikp == ikbot ! ******* 
     276                     
     277                     IF( ilayer > 0 ) THEN      !! least an unstable layer has been found 
     278                        ! 
     279                        IF( lp_monitor_point ) THEN 
     280                           WRITE(numout,*) 
     281                           WRITE(numout,*) 'After ',jiter,' iteration(s), we neutralized ',ilayer,' instable layer(s)' 
     282                           WRITE(numout,*) '   ==> N2 at i,j=',ji,',',jj,' now looks like this:' 
     283                           DO jk = 1, klc1 
     284                              WRITE(numout,*) jk, zvn2(jk) 
     285                           END DO 
     286                           WRITE(numout,*) 
     287                        ENDIF 
     288                        ! 
     289                        ikp    = 1     ! starting again at the surface for the next iteration 
     290                        ilayer = 0 
     291                     ENDIF 
    142292                     ! 
    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 ) 
     293                     IF( ikp >= ikbot )   l_column_treated = .TRUE. 
     294                     ! 
     295                  END DO ! DO WHILE ( .NOT. l_column_treated ) 
     296 
     297                  !! Updating tsa: 
     298                  tsa(ji,jj,:,jp_tem) = zvts(:,jp_tem) 
     299                  tsa(ji,jj,:,jp_sal) = zvts(:,jp_sal) 
     300 
     301                  !! LB:  Potentially some other global variable beside theta and S can be treated here 
     302                  !!      like BGC tracers. 
     303 
     304                  IF( lp_monitor_point )   WRITE(numout,*) 
     305 
     306               ENDIF ! IF( tmask(ji,jj,3) == 1 ) THEN 
     307 
     308            END DO ! ji 
     309         END DO ! jj 
     310         ! 
     311         IF( l_trdtra ) THEN         ! send the Non penetrative mixing trends for diagnostic 
     312            z1_r2dt = 1._wp / (2._wp * rdt) 
     313            ztrdt(:,:,:) = ( tsa(:,:,:,jp_tem) - ztrdt(:,:,:) ) * z1_r2dt 
     314            ztrds(:,:,:) = ( tsa(:,:,:,jp_sal) - ztrds(:,:,:) ) * z1_r2dt 
     315            CALL trd_tra( kt, 'TRA', jp_tem, jptra_npc, ztrdt ) 
     316            CALL trd_tra( kt, 'TRA', jp_sal, jptra_npc, ztrds ) 
    203317            CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
    204318         ENDIF 
    205        
    206          ! Lateral boundary conditions on ( ta, sa )   ( Unchanged sign) 
    207          ! ------------------------------============ 
     319         ! 
    208320         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 
     321         ! 
     322         IF( lwp .AND. l_LB_debug ) THEN 
     323            WRITE(numout,*) 'Exiting tra_npc , kt = ',kt,', => numb. of statically instable water-columns: ', inpcc 
     324            WRITE(numout,*) 
    216325         ENDIF 
    217326         ! 
    218       ENDIF 
    219       ! 
    220       CALL wrk_dealloc(jpi, jpj, jpk, zrhop ) 
    221       CALL wrk_dealloc(jpi, jpk, zwx, zwy, zwz ) 
     327         CALL wrk_dealloc(jpi, jpj, jpk, zn2 ) 
     328         CALL wrk_dealloc(jpi, jpj, jpk, 2, zab ) 
     329         CALL wrk_dealloc(jpk, zvn2 ) 
     330         CALL wrk_dealloc(jpk, 2, zvts, zvab ) 
     331         ! 
     332      ENDIF   ! IF( MOD( kt, nn_npc ) == 0 ) THEN 
    222333      ! 
    223334      IF( nn_timing == 1 )  CALL timing_stop('tra_npc') 
Note: See TracChangeset for help on using the changeset viewer.