Changeset 4680


Ignore:
Timestamp:
2014-06-20T12:46:45+02:00 (6 years ago)
Author:
gm
Message:

#1350 : new version of tranpc.F90 based local static stability

Location:
branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/LIM_SRC_2/limhdf_2.F90

    r3625 r4680  
    8686      zdiv0(:, 1 ) = 0._wp 
    8787      zdiv0(:,jpj) = 0._wp 
    88       IF( .NOT.lk_vopt_loop ) THEN 
    89          zflu (jpi,:) = 0._wp    
    90          zflv (jpi,:) = 0._wp 
    91          zdiv0(1,  :) = 0._wp 
    92          zdiv0(jpi,:) = 0._wp 
    93       ENDIF 
     88      zflu (jpi,:) = 0._wp    
     89      zflv (jpi,:) = 0._wp 
     90      zdiv0(1,  :) = 0._wp 
     91      zdiv0(jpi,:) = 0._wp 
    9492 
    9593      zconv = 1._wp           !==  horizontal diffusion using a Crant-Nicholson scheme  ==! 
  • branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/LIM_SRC_2/limistate_2.F90

    r4147 r4680  
    1414   !!   'key_lim2' :                                  LIM 2.0 sea-ice model 
    1515   !!---------------------------------------------------------------------- 
    16    !!---------------------------------------------------------------------- 
    1716   !!   lim_istate_2      :  Initialisation of diagnostics ice variables 
    1817   !!   lim_istate_init_2 :  initialization of ice state and namelist read 
     
    3433   PUBLIC lim_istate_2      ! routine called by lim_init_2.F90 
    3534 
    36    !!! **  namelist (namiceini) ** 
    37    LOGICAL  ::   ln_limini   !: Ice initialization state 
     35   !                        !! **  namelist (namiceini) ** 
     36   LOGICAL  ::   ln_limini   ! Ice initialization state 
    3837   REAL(wp) ::   ttest       ! threshold water temperature for initial sea ice 
    3938   REAL(wp) ::   hninn       ! initial snow thickness in the north 
     
    5150   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    5251   !!---------------------------------------------------------------------- 
    53  
    5452CONTAINS 
    5553 
     
    7169      IF( .NOT. ln_limini ) THEN   
    7270          
    73          tfu(:,:) = tfreez( tsn(:,:,1,jp_sal) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius] 
     71         tfu(:,:) = eos_fzp( tsn(:,:,1,jp_sal) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius] 
    7472 
    7573         DO jj = 1, jpj 
  • branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90

    r4333 r4680  
    8383      zdiv0(:, 1 ) = 0._wp 
    8484      zdiv0(:,jpj) = 0._wp 
    85       IF( .NOT.lk_vopt_loop ) THEN 
    86          zflu (jpi,:) = 0._wp    
    87          zflv (jpi,:) = 0._wp 
    88          zdiv0(1,  :) = 0._wp 
    89          zdiv0(jpi,:) = 0._wp 
    90       ENDIF 
     85      zflu (jpi,:) = 0._wp    
     86      zflv (jpi,:) = 0._wp 
     87      zdiv0(1,  :) = 0._wp 
     88      zdiv0(jpi,:) = 0._wp 
    9189 
    9290      zconv = 1._wp           !==  horizontal diffusion using a Crant-Nicholson scheme  ==! 
  • branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

    r4335 r4680  
    3636   PUBLIC   lim_istate      ! routine called by lim_init.F90 
    3737 
    38    !! * Module variables 
    3938   !                          !!** init namelist (namiceini) ** 
    4039   REAL(wp) ::   ttest   ! threshold water temperature for initial sea ice 
     
    5352   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    5453   !!---------------------------------------------------------------------- 
    55  
    5654CONTAINS 
    5755 
     
    121119 
    122120      ! Basal temperature is set to the freezing point of seawater in Celsius 
    123       t_bo(:,:) = tfreez( tsn(:,:,1,jp_sal) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius] 
     121      t_bo(:,:) = eos_fzp( tsn(:,:,1,jp_sal) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius] 
    124122 
    125123      DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest 
  • branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90

    r4619 r4680  
    551551      DO jj = 1, jpjm1 
    552552         DO ji = 1, jpim1 
    553             mgrhu(ji,jj) = INT(  SIGN( 1.e0, fsdept_0(ji+1,jj,mbkt(ji+1,jj)) - fsdept_0(ji,jj,mbkt(ji,jj)) )  ) 
    554             mgrhv(ji,jj) = INT(  SIGN( 1.e0, fsdept_0(ji,jj+1,mbkt(ji,jj+1)) - fsdept_0(ji,jj,mbkt(ji,jj)) )  ) 
     553            mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  ) 
     554            mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  ) 
    555555         END DO 
    556556      END DO 
     
    558558      DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point 
    559559         DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0) 
    560             e3u_bbl_0(ji,jj) = MIN( fse3u_0(ji,jj,mbkt(ji+1,jj  )), fse3u_0(ji,jj,mbkt(ji,jj)) ) 
    561             e3v_bbl_0(ji,jj) = MIN( fse3v_0(ji,jj,mbkt(ji  ,jj+1)), fse3v_0(ji,jj,mbkt(ji,jj)) ) 
     560            e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) ) 
     561            e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) ) 
    562562         END DO 
    563563      END DO 
  • branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90

    r4619 r4680  
    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 
    1921   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) 
    2224   USE lbclnk          ! lateral boundary conditions (or mpp link) 
    2325   USE in_out_manager  ! I/O manager 
     
    2931   PRIVATE 
    3032 
    31    PUBLIC   tra_npc       ! routine called by step.F90 
     33   PUBLIC   tra_npc    ! routine called by step.F90 
    3234 
    3335   !! * Substitutions 
    3436#  include "domzgr_substitute.h90" 
    35    !!---------------------------------------------------------------------- 
    36    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    37    !! $Id$  
     37#  include "vectopt_loop_substitute.h90" 
     38   !!---------------------------------------------------------------------- 
     39   !! NEMO/OPA 3.7 , NEMO Consortium (2014) 
     40   !! $Id$ 
    3841   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    3942   !!---------------------------------------------------------------------- 
     
    4447      !!                  ***  ROUTINE tranpc  *** 
    4548      !! 
    46       !! ** Purpose :   Non penetrative convective adjustment scheme. solve  
     49      !! ** Purpose : Non-penetrative convective adjustment scheme. solve 
    4750      !!      the static instability of the water column on after fields 
    4851      !!      while conserving heat and salt contents. 
    4952      !! 
    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. 
     53      !! ** Method  : updated algorithm able to deal with non-linear equation of state 
     54      !!              (i.e. static stability computed locally) 
    5455      !! 
    5556      !! ** Action  : - (ta,sa) after the application od the npc scheme 
    56       !!              - save the associated trends (ttrd,strd) ('key_trdtra') 
     57      !!              - send the associated trends for on-line diagnostics (l_trdtra=T) 
    5758      !! 
    58       !! References : Madec, et al., 1991, JPO, 21, 9, 1349-1371. 
     59      !! References :     Madec, et al., 1991, JPO, 21, 9, 1349-1371. 
    5960      !!---------------------------------------------------------------------- 
    60       ! 
    6161      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    6262      ! 
    6363      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    6464      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 
     65      INTEGER  ::   jiter, ikbot, ik, ikup, ikdown, ilayer, ikm   ! local integers 
     66      LOGICAL  ::   l_bottom_reached, l_column_treated 
     67      REAL(wp) ::   zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z 
     68      REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_r2dt 
     69      REAL(wp), POINTER, DIMENSION(:)       ::   zvn2   ! vertical profile of N2 at 1 given point... 
     70      REAL(wp), POINTER, DIMENSION(:,:)     ::   zvts   ! vertical profile of T and S at 1 given point... 
     71      REAL(wp), POINTER, DIMENSION(:,:)     ::   zvab   ! vertical profile of alpha and beta 
     72      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zn2    ! N^2  
     73      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zab    ! alpha and beta 
     74      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   ztrdt, ztrds   ! 3D workspace 
     75      ! 
     76      !!LB debug: 
     77      LOGICAL, PARAMETER :: l_LB_debug = .FALSE. 
     78      INTEGER :: ilc1, jlc1, klc1, nncpu 
     79      LOGICAL :: lp_monitor_point 
     80      !!LB debug. 
    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         !LB debug: 
     99         IF( lwp .AND. l_LB_debug ) THEN 
     100            WRITE(numout,*) '' ; WRITE(numout,*) '' 
     101            WRITE(numout,*) 'LOLO: entering tra_npc, kt, narea =', kt, narea 
     102         ENDIF 
     103         !LBdebug: Monitoring of 1 column subject to convection... 
     104         IF( l_LB_debug ) THEN 
     105            ! Location of 1 known convection spot to follow what's happening in the water column 
     106            ilc1 = 54 ;  jlc1 = 15 ; ! Labrador ORCA1 4x4 cpus: 
     107            nncpu = 15  ; ! the CPU domain contains the convection spot 
     108            !ilc1 = 14 ;  jlc1 = 13 ; ! Labrador ORCA1 8x8 cpus: 
     109            !nncpu = 54  ; ! the CPU domain contains the convection spot 
     110            klc1 =  mbkt(ilc1,jlc1) ! bottom of the ocean for debug point...           
     111         ENDIF 
     112         !LBdebug. 
     113 
     114         CALL eos_rab( tsa, zab )         ! after alpha and beta 
     115         CALL bn2    ( tsa, zab, zn2 )    ! after Brunt-Vaisala 
     116         
     117         inpcc = 0 
     118 
     119         DO jj = 2, jpjm1                 ! interior column only 
     120            DO ji = fs_2, fs_jpim1 
     121               ! 
     122               IF( tmask(ji,jj,2) == 1 ) THEN      ! At least 2 ocean points 
     123                  !                                     ! consider one ocean column  
     124                  zvts(:,jp_tem) = tsa(ji,jj,:,jp_tem)      ! temperature 
     125                  zvts(:,jp_sal) = tsa(ji,jj,:,jp_sal)      ! salinity 
     126 
     127                  zvab(:,jp_tem)  = zab(ji,jj,:,jp_tem)     ! Alpha  
     128                  zvab(:,jp_sal)  = zab(ji,jj,:,jp_sal)     ! Beta   
     129                  zvn2(:)         = zn2(ji,jj,:)            ! N^2  
     130                  
     131                  IF( l_LB_debug ) THEN                  !LB debug: 
     132                     lp_monitor_point = .FALSE. 
     133                     IF( ( ji == ilc1 ).AND.( jj == jlc1 ) ) lp_monitor_point = .TRUE. 
     134                     ! writing only if on CPU domain where conv region is: 
     135                     lp_monitor_point = (narea == nncpu).AND.lp_monitor_point  
     136                      
     137                     IF(lp_monitor_point) THEN 
     138                        WRITE(numout,*) '' ;WRITE(numout,*) '' ; 
     139                       WRITE(numout,'("Time step = ",i6.6," !!!")')  kt 
     140                        WRITE(numout,'(" *** BEFORE anything, N^2 for point ",i3,",",i3,":" )') , ji,jj 
     141                        DO jk = 1, klc1 
     142                           WRITE(numout,*) jk, zvn2(jk) 
     143                        END DO 
     144                        WRITE(numout,*) ' ' 
     145                     ENDIF 
     146                  ENDIF                                  !LB debug  end 
     147 
     148                  ikbot = mbkt(ji,jj)   ! ikbot: ocean bottom T-level 
     149                  ik = 1                ! because N2 is irrelevant at the surface level (will start at ik=2) 
     150                  ilayer = 0 
     151                  jiter  = 0 
     152                  l_column_treated = .FALSE. 
     153                  
     154                  DO WHILE ( .NOT. l_column_treated ) 
    140155                     ! 
    141                      ikbot = mbkt(ji,jj)        ! ikbot: ocean bottom T-level 
     156                     jiter = jiter + 1 
     157                     
     158                     IF( jiter >= 400 ) EXIT 
     159                     
     160                     l_bottom_reached = .FALSE. 
     161 
     162                     DO WHILE ( .NOT. l_bottom_reached ) 
     163 
     164                        ik = ik + 1 
     165                        
     166                        !! Checking level ik for instability 
     167                        !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
     168 
     169                        IF( zvn2(ik) < 0. ) THEN ! Instability found! 
     170 
     171                           ikm  = ik              ! first level whith negative N2 
     172                           ilayer = ilayer + 1    ! yet another layer found.... 
     173                           IF(jiter == 1) inpcc = inpcc + 1 
     174 
     175                           IF(l_LB_debug .AND. lp_monitor_point) & 
     176                              & WRITE(numout,*) 'Negative N2 at ik =', ikm, ' layer nb.', ilayer, & 
     177                              & ' inpcc =', inpcc 
     178 
     179                           !! Case we mix with upper regions where N2==0: 
     180                           !! All the points above ikup where N2 == 0 must also be mixed => we go 
     181                           !! upward to find a new ikup, where the layer doesn't have N2==0 
     182                           ikup = ikm 
     183                           DO jk = ikm, 2, -1 
     184                              ikup = ikup - 1 
     185                              IF( (zvn2(jk-1) > 0.).OR.(ikup == 1) ) EXIT 
     186                           END DO 
     187                           
     188                           ! adjusting ikup if the upper part of the unstable column was neutral (N2=0) 
     189                           IF((zvn2(ikup+1) == 0.).AND.(ikup /= 1)) ikup = ikup+1 ; 
     190 
     191                           
     192                           IF(lp_monitor_point) WRITE(numout,*) ' => ikup is =', ikup, ' layer nb.', ilayer 
     193                           
     194                           zsum_temp = 0._wp 
     195                           zsum_sali = 0._wp 
     196                           zsum_alfa = 0._wp 
     197                           zsum_beta = 0._wp 
     198                           zsum_z    = 0._wp 
     199                                                     
     200                           DO jk = ikup, ikbot+1      ! Inside the instable (and overlying neutral) portion of the column 
     201                              ! 
     202                              IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '     -> summing for jk =', jk 
     203                              ! 
     204                              zdz       = fse3t(ji,jj,jk) 
     205                              zsum_temp = zsum_temp + zvts(jk,jp_tem)*zdz 
     206                              zsum_sali = zsum_sali + zvts(jk,jp_sal)*zdz 
     207                              zsum_alfa = zsum_alfa + zvab(jk,jp_tem)*zdz 
     208                              zsum_beta = zsum_beta + zvab(jk,jp_sal)*zdz 
     209                              zsum_z    = zsum_z    + zdz 
     210                              ! 
     211                              !! EXIT if we found the bottom of the unstable portion of the water column     
     212                              IF( (zvn2(jk+1) > 0.).OR.(jk == ikbot ).OR.((jk==ikm).AND.(zvn2(jk+1) == 0.)) )   EXIT 
     213                           END DO 
     214                           
     215                           !ik     = jk !LB remove? 
     216                           ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative N2 
     217                           
     218                           IF(l_LB_debug .AND. lp_monitor_point) & 
     219                              &    WRITE(numout,*) '  => ikdown =', ikdown, '  layer nb.', ilayer 
     220                           
     221                           ! Mixing Temperature and salinity between ikup and ikdown: 
     222                           zta   = zsum_temp/zsum_z 
     223                           zsa   = zsum_sali/zsum_z 
     224                           zalfa = zsum_alfa/zsum_z 
     225                           zbeta = zsum_beta/zsum_z 
     226 
     227                           IF(l_LB_debug .AND. lp_monitor_point) THEN 
     228                              WRITE(numout,*) '  => Mean temp. in that portion =', zta 
     229                              WRITE(numout,*) '  => Mean sali. in that portion =', zsa 
     230                              WRITE(numout,*) '  => Mean Alpha in that portion =', zalfa 
     231                              WRITE(numout,*) '  => Mean Beta  in that portion =', zbeta 
     232                           ENDIF 
     233 
     234                           !! Homogenaizing the temperature, salinity, alpha and beta in this portion of the column 
     235                           DO jk = ikup, ikdown 
     236                              zvts(jk,jp_tem) = zta 
     237                              zvts(jk,jp_sal) = zsa 
     238                              zvab(jk,jp_tem) = zalfa 
     239                              zvab(jk,jp_sal) = zbeta 
     240                           END DO 
     241                           ! 
     242                           !! Before updating N2, it is possible that another unstable 
     243                           !! layer exists underneath the one we just homogeneized! 
     244                           ik = ikdown 
     245                           !  
     246                        ENDIF  ! IF( zvn2(ik+1) < 0. ) THEN 
     247                        ! 
     248                        IF( ik == ikbot ) l_bottom_reached = .TRUE. 
     249                        ! 
     250                     END DO ! DO WHILE ( .NOT. l_bottom_reached ) 
     251 
     252                     IF( ik /= ikbot )   STOP 'ERROR: tranpc.F90 => PROBLEM #1' 
     253                     
     254                     ! ******* At this stage ik == ikbot ! ******* 
     255                     
     256                     IF( ilayer > 0 ) THEN 
     257                        !! least an unstable layer has been found 
     258                        !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion 
     259                        !! => Need to re-compute N2! will use Alpha and Beta! 
     260                        ! 
     261                        DO jk = ikup+1, ikdown+1   ! we must go 1 point deeper than ikdown!      
     262                           !! Doing exactly as in eosbn2.F90: 
     263                           !! * Except that we only are interested in the sign of N2 !!! 
     264                           !!   => just considering the vertical gradient of density 
     265                           zrw =   (fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk)) & 
     266                               & / (fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk)) 
     267                           zaw = zvab(jk,jp_tem) * (1._wp - zrw) + zvab(jk-1,jp_tem) * zrw 
     268                           zbw = zvab(jk,jp_sal) * (1._wp - zrw) + zvab(jk-1,jp_sal) * zrw 
     269                           
     270                           !zvn2(jk) = grav*( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     & 
     271                           !     &           - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )  & 
     272                           !     &       / fse3w(ji,jj,jk) * tmask(ji,jj,jk) 
     273                           zvn2(jk) = (  zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     & 
     274                                &      - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )   
     275                        END DO 
     276 
     277                        IF(l_LB_debug .AND. lp_monitor_point) THEN 
     278                           WRITE(numout, '(" *** After iteration #",i3.3,", N^2 for point ",i3,",",i3,":" )') & 
     279                              & jiter, ji,jj 
     280                           DO jk = 1, klc1 
     281                              WRITE(numout,*) jk, zvn2(jk) 
     282                           END DO 
     283                           WRITE(numout,*) ' ' 
     284                        ENDIF 
     285 
     286                        ik     = 1  ! starting again at the surface for the next iteration 
     287                        ilayer = 0 
     288                     ENDIF 
    142289                     ! 
    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(:,:,:) 
     290                     IF( ik >= ikbot ) THEN 
     291                        IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '    --- exiting jiter loop ---' 
     292                        l_column_treated = .TRUE. 
     293                     ENDIF 
     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                  !! lolo:  Should we update something else???? 
     302                  !! => like alpha and beta? 
     303 
     304                  IF(l_LB_debug .AND. 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 
    201315            CALL trd_tra( kt, 'TRA', jp_tem, jptra_npc, ztrdt ) 
    202316            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 
    216          ENDIF 
    217          ! 
    218       ENDIF 
    219       ! 
    220       CALL wrk_dealloc(jpi, jpj, jpk, zrhop ) 
    221       CALL wrk_dealloc(jpi, jpk, zwx, zwy, zwz ) 
     321         ! 
     322         IF(lwp) THEN 
     323            WRITE(numout,*) 'LOLO: exiting tra_npc, kt =', kt 
     324            WRITE(numout,*)' => number of statically instable water column : ',inpcc 
     325            WRITE(numout,*) '' ; WRITE(numout,*) '' 
     326         ENDIF 
     327         ! 
     328         CALL wrk_dealloc(jpi, jpj, jpk, zn2 ) 
     329         CALL wrk_dealloc(jpi, jpj, jpk, 2, zab ) 
     330         CALL wrk_dealloc(jpk, zvn2 ) 
     331         CALL wrk_dealloc(jpk, 2, zvts, zvab ) 
     332         ! 
     333      ENDIF   ! IF( MOD( kt, nn_npc ) == 0 ) THEN 
    222334      ! 
    223335      IF( nn_timing == 1 )  CALL timing_stop('tra_npc') 
Note: See TracChangeset for help on using the changeset viewer.