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.
tranpc.F90 in NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRA – NEMO

source: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRA/tranpc.F90 @ 12928

Last change on this file since 12928 was 12928, checked in by smueller, 4 years ago

Synchronizing with /NEMO/trunk@12925 (ticket #2170)

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