source: NEMO/trunk/src/OCE/TRA/tranpc.F90 @ 13237

Last change on this file since 13237 was 13237, checked in by smasson, 3 months ago

trunk: Mid-year merge, merge back KERNEL-06_techene_e3

  • Property svn:keywords set to Id
File size: 16.8 KB
Line 
1MODULE tranpc
2   !!==============================================================================
3   !!                       ***  MODULE  tranpc  ***
4   !! Ocean active tracers:  non penetrative convective adjustment scheme
5   !!==============================================================================
6   !! History :  1.0  ! 1990-09  (G. Madec)  Original code
7   !!                 ! 1996-01  (G. Madec)  statement function for e3
8   !!   NEMO     1.0  ! 2002-06  (G. Madec)  free form F90
9   !!            3.0  ! 2008-06  (G. Madec)  applied on ta, sa and called before tranxt in step.F90
10   !!            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.
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   tra_npc       : apply the non penetrative convection scheme
16   !!----------------------------------------------------------------------
17   USE oce            ! ocean dynamics and active tracers
18   USE dom_oce        ! ocean space and time domain
19   USE phycst         ! physical constants
20   USE zdf_oce        ! ocean vertical physics
21   USE trd_oce        ! ocean active tracer trends
22   USE trdtra         ! ocean active tracer trends
23   USE eosbn2         ! equation of state (eos routine)
24   !
25   USE lbclnk         ! lateral boundary conditions (or mpp link)
26   USE in_out_manager ! I/O manager
27   USE lib_mpp        ! MPP library
28   USE timing         ! Timing
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   tra_npc    ! routine called by step.F90
34
35   !! * Substitutions
36#  include "do_loop_substitute.h90"
37#  include "domzgr_substitute.h90"
38   !!----------------------------------------------------------------------
39   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
40   !! $Id$
41   !! Software governed by the CeCILL license (see ./LICENSE)
42   !!----------------------------------------------------------------------
43CONTAINS
44
45   SUBROUTINE tra_npc( kt, Kmm, Krhs, pts, Kaa )
46      !!----------------------------------------------------------------------
47      !!                  ***  ROUTINE tranpc  ***
48      !!
49      !! ** Purpose : Non-penetrative convective adjustment scheme. solve
50      !!      the static instability of the water column on after fields
51      !!      while conserving heat and salt contents.
52      !!
53      !! ** Method  : updated algorithm able to deal with non-linear equation of state
54      !!              (i.e. static stability computed locally)
55      !!
56      !! ** Action  : - tsa: after tracers with the application of the npc scheme
57      !!              - send the associated trends for on-line diagnostics (l_trdtra=T)
58      !!
59      !! References :     Madec, et al., 1991, JPO, 21, 9, 1349-1371.
60      !!----------------------------------------------------------------------
61      INTEGER,                                   INTENT(in   ) :: kt              ! ocean time-step index
62      INTEGER,                                   INTENT(in   ) :: Kmm, Krhs, Kaa  ! time level indices
63      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts             ! active tracers and RHS of tracer equation
64      !
65      INTEGER  ::   ji, jj, jk   ! dummy loop indices
66      INTEGER  ::   inpcc        ! number of statically instable water column
67      INTEGER  ::   jiter, ikbot, ikp, ikup, ikdown, ilayer, ik_low   ! local integers
68      LOGICAL  ::   l_bottom_reached, l_column_treated
69      REAL(wp) ::   zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z
70      REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_rDt
71      REAL(wp), PARAMETER ::   zn2_zero = 1.e-14_wp             ! acceptance criteria for neutrality (N2==0)
72      REAL(wp), DIMENSION(        jpk     )   ::   zvn2         ! vertical profile of N2 at 1 given point...
73      REAL(wp), DIMENSION(        jpk,jpts)   ::   zvts, zvab   ! vertical profile of T & S , and  alpha & betaat 1 given point
74      REAL(wp), DIMENSION(jpi,jpj,jpk     )   ::   zn2          ! N^2
75      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts)   ::   zab          ! alpha and beta
76      REAL(wp), ALLOCATABLE, 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"
81      !!----------------------------------------------------------------------
82      !
83      IF( ln_timing )   CALL timing_start('tra_npc')
84      !
85      IF( MOD( kt, nn_npc ) == 0 ) THEN
86         !
87         IF( l_trdtra )   THEN                    !* Save initial after fields
88            ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
89            ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kaa) 
90            ztrds(:,:,:) = pts(:,:,:,jp_sal,Kaa)
91         ENDIF
92         !
93         IF( l_LB_debug ) THEN
94            ! Location of 1 known convection site to follow what's happening in the water column
95            ilc1 = 45 ;  jlc1 = 3 ; !  ORCA2 4x4, Antarctic coast, more than 2 unstable portions in the  water column...           
96            nncpu = 1  ;            ! the CPU domain contains the convection spot
97            klc1 =  mbkt(ilc1,jlc1) ! bottom of the ocean for debug point...         
98         ENDIF
99         !
100         CALL eos_rab( pts(:,:,:,:,Kaa), zab, Kmm )         ! after alpha and beta (given on T-points)
101         CALL bn2    ( pts(:,:,:,:,Kaa), zab, zn2, Kmm )    ! after Brunt-Vaisala  (given on W-points)
102         !
103         inpcc = 0
104         !
105         DO_2D_00_00
106            !
107            IF( tmask(ji,jj,2) == 1 ) THEN      ! At least 2 ocean points
108               !                                     ! consider one ocean column
109               zvts(:,jp_tem) = pts(ji,jj,:,jp_tem,Kaa)      ! temperature
110               zvts(:,jp_sal) = pts(ji,jj,:,jp_sal,Kaa)      ! salinity
111               !
112               zvab(:,jp_tem)  = zab(ji,jj,:,jp_tem)     ! Alpha
113               zvab(:,jp_sal)  = zab(ji,jj,:,jp_sal)     ! Beta 
114               zvn2(:)         = zn2(ji,jj,:)            ! N^2
115               !
116               IF( l_LB_debug ) THEN                  !LB debug:
117                  lp_monitor_point = .FALSE.
118                  IF( ( ji == ilc1 ).AND.( jj == jlc1 ) ) lp_monitor_point = .TRUE.
119                  ! writing only if on CPU domain where conv region is:
120                  lp_monitor_point = (narea == nncpu).AND.lp_monitor_point                     
121               ENDIF                                  !LB debug  end
122               !
123               ikbot = mbkt(ji,jj)   ! ikbot: ocean bottom T-level
124               ikp = 1                  ! because N2 is irrelevant at the surface level (will start at ikp=2)
125               ilayer = 0
126               jiter  = 0
127               l_column_treated = .FALSE.
128               !
129               DO WHILE ( .NOT. l_column_treated )
130                  !
131                  jiter = jiter + 1
132                  !
133                  IF( jiter >= 400 ) EXIT
134                  !
135                  l_bottom_reached = .FALSE.
136                  !
137                  DO WHILE ( .NOT. l_bottom_reached )
138                     !
139                     ikp = ikp + 1
140                     !
141                     !! Testing level ikp for instability
142                     !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
143                     IF( zvn2(ikp) <  -zn2_zero ) THEN ! Instability found!
144                        !
145                        ilayer = ilayer + 1    ! yet another instable portion of the water column found....
146                        !
147                        IF( lp_monitor_point ) THEN
148                           WRITE(numout,*)
149                           IF( ilayer == 1 .AND. jiter == 1 ) THEN   ! first time a column is spoted with an instability
150                              WRITE(numout,*)
151                              WRITE(numout,*) 'Time step = ',kt,' !!!'
152                           ENDIF
153                           WRITE(numout,*)  ' * Iteration #',jiter,': found instable portion #',ilayer,   &
154                              &                                    ' in column! Starting at ikp =', ikp
155                           WRITE(numout,*)  ' *** N2 for point (i,j) = ',ji,' , ',jj
156                           DO jk = 1, klc1
157                              WRITE(numout,*) jk, zvn2(jk)
158                           END DO
159                           WRITE(numout,*)
160                        ENDIF
161                        !
162                        IF( jiter == 1 )   inpcc = inpcc + 1 
163                        !
164                        IF( lp_monitor_point )   WRITE(numout, *) 'Negative N2 at ikp =',ikp,' for layer #', ilayer
165                        !
166                        !! ikup is the uppermost point where mixing will start:
167                        ikup = ikp - 1 ! ikup is always "at most at ikp-1", less if neutral levels overlying
168                        !
169                        !! If the points above ikp-1 have N2 == 0 they must also be mixed:
170                        IF( ikp > 2 ) THEN
171                           DO jk = ikp-1, 2, -1
172                              IF( ABS(zvn2(jk)) < zn2_zero ) THEN
173                                 ikup = ikup - 1  ! 1 more upper level has N2=0 and must be added for the mixing
174                              ELSE
175                                 EXIT
176                              ENDIF
177                           END DO
178                        ENDIF
179                        !
180                        IF( ikup < 1 )   CALL ctl_stop( 'tra_npc :  PROBLEM #1')
181                        !
182                        zsum_temp = 0._wp
183                        zsum_sali = 0._wp
184                        zsum_alfa = 0._wp
185                        zsum_beta = 0._wp
186                        zsum_z    = 0._wp
187                                                 
188                        DO jk = ikup, ikbot      ! Inside the instable (and overlying neutral) portion of the column
189                           !
190                           zdz       = e3t(ji,jj,jk,Kmm)
191                           zsum_temp = zsum_temp + zvts(jk,jp_tem)*zdz
192                           zsum_sali = zsum_sali + zvts(jk,jp_sal)*zdz
193                           zsum_alfa = zsum_alfa + zvab(jk,jp_tem)*zdz
194                           zsum_beta = zsum_beta + zvab(jk,jp_sal)*zdz
195                           zsum_z    = zsum_z    + zdz
196                           !                             
197                           IF( jk == ikbot ) EXIT ! avoid array-index overshoot in case ikbot = jpk, cause we're calling jk+1 next line
198                           !! EXIT when we have reached the last layer that is instable (N2<0) or neutral (N2=0):
199                           IF( zvn2(jk+1) > zn2_zero ) EXIT
200                        END DO
201                       
202                        ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative or neutral N2
203                        IF( ikup == ikdown )   CALL ctl_stop( 'tra_npc :  PROBLEM #2')
204
205                        ! Mixing Temperature, salinity, alpha and beta from ikup to ikdown included:
206                        zta   = zsum_temp/zsum_z
207                        zsa   = zsum_sali/zsum_z
208                        zalfa = zsum_alfa/zsum_z
209                        zbeta = zsum_beta/zsum_z
210
211                        IF( lp_monitor_point ) THEN
212                           WRITE(numout,*) 'MIXED T, S, alfa and beta between ikup =',ikup,   &
213                              &            ' and ikdown =',ikdown,', in layer #',ilayer
214                           WRITE(numout,*) '  => Mean temp. in that portion =', zta
215                           WRITE(numout,*) '  => Mean sali. in that portion =', zsa
216                           WRITE(numout,*) '  => Mean Alfa  in that portion =', zalfa
217                           WRITE(numout,*) '  => Mean Beta  in that portion =', zbeta
218                        ENDIF
219
220                        !! Homogenaizing the temperature, salinity, alpha and beta in this portion of the column
221                        DO jk = ikup, ikdown
222                           zvts(jk,jp_tem) = zta
223                           zvts(jk,jp_sal) = zsa
224                           zvab(jk,jp_tem) = zalfa
225                           zvab(jk,jp_sal) = zbeta
226                        END DO
227                       
228                       
229                        !! Updating N2 in the relvant portion of the water column
230                        !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion
231                        !! => Need to re-compute N2! will use Alpha and Beta!
232                       
233                        ikup   = MAX(2,ikup)         ! ikup can never be 1 !
234                        ik_low = MIN(ikdown+1,ikbot) ! we must go 1 point deeper than ikdown!
235                       
236                        DO jk = ikup, ik_low              ! we must go 1 point deeper than ikdown!
237
238                           !! Interpolating alfa and beta at W point:
239                           zrw =  (gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm)) &
240                              & / (gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm))
241                           zaw = zvab(jk,jp_tem) * (1._wp - zrw) + zvab(jk-1,jp_tem) * zrw
242                           zbw = zvab(jk,jp_sal) * (1._wp - zrw) + zvab(jk-1,jp_sal) * zrw
243
244                           !! N2 at W point, doing exactly as in eosbn2.F90:
245                           zvn2(jk) = grav*( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     &
246                              &            - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )  &
247                              &       / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
248
249                           !! OR, faster  => just considering the vertical gradient of density
250                           !! as only the signa maters...
251                           !zvn2(jk) = (  zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     &
252                           !     &      - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )
253
254                        END DO
255                     
256                        ikp = MIN(ikdown+1,ikbot)
257                       
258
259                     ENDIF  !IF( zvn2(ikp) < 0. )
260
261
262                     IF( ikp == ikbot ) l_bottom_reached = .TRUE.
263                     !
264                  END DO ! DO WHILE ( .NOT. l_bottom_reached )
265
266                  IF( ikp /= ikbot )   CALL ctl_stop( 'tra_npc :  PROBLEM #3')
267                 
268                  ! ******* At this stage ikp == ikbot ! *******
269                 
270                  IF( ilayer > 0 ) THEN      !! least an unstable layer has been found
271                     !
272                     IF( lp_monitor_point ) THEN
273                        WRITE(numout,*)
274                        WRITE(numout,*) 'After ',jiter,' iteration(s), we neutralized ',ilayer,' instable layer(s)'
275                        WRITE(numout,*) '   ==> N2 at i,j=',ji,',',jj,' now looks like this:'
276                        DO jk = 1, klc1
277                           WRITE(numout,*) jk, zvn2(jk)
278                        END DO
279                        WRITE(numout,*)
280                     ENDIF
281                     !
282                     ikp    = 1     ! starting again at the surface for the next iteration
283                     ilayer = 0
284                  ENDIF
285                  !
286                  IF( ikp >= ikbot )   l_column_treated = .TRUE.
287                  !
288               END DO ! DO WHILE ( .NOT. l_column_treated )
289
290               !! Updating pts:
291               pts(ji,jj,:,jp_tem,Kaa) = zvts(:,jp_tem)
292               pts(ji,jj,:,jp_sal,Kaa) = zvts(:,jp_sal)
293
294               !! LB:  Potentially some other global variable beside theta and S can be treated here
295               !!      like BGC tracers.
296
297               IF( lp_monitor_point )   WRITE(numout,*)
298
299            ENDIF ! IF( tmask(ji,jj,3) == 1 ) THEN
300
301         END_2D
302         !
303         IF( l_trdtra ) THEN         ! send the Non penetrative mixing trends for diagnostic
304            z1_rDt = 1._wp / (2._wp * rn_Dt)
305            ztrdt(:,:,:) = ( pts(:,:,:,jp_tem,Kaa) - ztrdt(:,:,:) ) * z1_rDt
306            ztrds(:,:,:) = ( pts(:,:,:,jp_sal,Kaa) - ztrds(:,:,:) ) * z1_rDt
307            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_npc, ztrdt )
308            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_npc, ztrds )
309            DEALLOCATE( ztrdt, ztrds )
310         ENDIF
311         !
312         CALL lbc_lnk_multi( 'tranpc', pts(:,:,:,jp_tem,Kaa), 'T', 1.0_wp, pts(:,:,:,jp_sal,Kaa), 'T', 1.0_wp )
313         !
314         IF( lwp .AND. l_LB_debug ) THEN
315            WRITE(numout,*) 'Exiting tra_npc , kt = ',kt,', => numb. of statically instable water-columns: ', inpcc
316            WRITE(numout,*)
317         ENDIF
318         !
319      ENDIF   ! IF( MOD( kt, nn_npc ) == 0 ) THEN
320      !
321      IF( ln_timing )   CALL timing_stop('tra_npc')
322      !
323   END SUBROUTINE tra_npc
324
325   !!======================================================================
326END MODULE tranpc
Note: See TracBrowser for help on using the repository browser.