source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 20 months ago

GMED 450 add flush after prints

File size: 17.7 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 wrk_nemo        ! Memory Allocation
29   USE timing          ! Timing
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   tra_npc    ! routine called by step.F90
35
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/OPA 3.6 , NEMO Consortium (2014)
41   !! $Id$
42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE tra_npc( kt )
47      !!----------------------------------------------------------------------
48      !!                  ***  ROUTINE tranpc  ***
49      !!
50      !! ** Purpose : Non-penetrative convective adjustment scheme. solve
51      !!      the static instability of the water column on after fields
52      !!      while conserving heat and salt contents.
53      !!
54      !! ** Method  : updated algorithm able to deal with non-linear equation of state
55      !!              (i.e. static stability computed locally)
56      !!
57      !! ** Action  : - (ta,sa) after the application od the npc scheme
58      !!              - send the associated trends for on-line diagnostics (l_trdtra=T)
59      !!
60      !! References :     Madec, et al., 1991, JPO, 21, 9, 1349-1371.
61      !!----------------------------------------------------------------------
62      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
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_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"
81      !!----------------------------------------------------------------------
82      !
83      IF( nn_timing == 1 )  CALL timing_start('tra_npc')
84      !
85      IF( MOD( kt, nn_npc ) == 0 ) THEN
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
93            CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
94            ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
95            ztrds(:,:,:) = tsa(:,:,:,jp_sal)
96         ENDIF
97
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 )
136                     !
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                              IF(lflush) CALL flush(numout)
167                           ENDIF
168                           
169
170                           IF( jiter == 1 )   inpcc = inpcc + 1 
171
172                           IF( lp_monitor_point )   THEN
173                              WRITE(numout, *) 'Negative N2 at ikp =',ikp,' for layer #', ilayer
174                              IF(lflush) CALL flush(numout)
175                           ENDIF
176
177                           !! ikup is the uppermost point where mixing will start:
178                           ikup = ikp - 1 ! ikup is always "at most at ikp-1", less if neutral levels overlying
179                           
180                           !! If the points above ikp-1 have N2 == 0 they must also be mixed:
181                           IF( ikp > 2 ) THEN
182                              DO jk = ikp-1, 2, -1
183                                 IF( ABS(zvn2(jk)) < zn2_zero ) THEN
184                                    ikup = ikup - 1  ! 1 more upper level has N2=0 and must be added for the mixing
185                                 ELSE
186                                    EXIT
187                                 ENDIF
188                              END DO
189                           ENDIF
190                           
191                           IF( ikup < 1 )   CALL ctl_stop( 'tra_npc :  PROBLEM #1')
192
193                           zsum_temp = 0._wp
194                           zsum_sali = 0._wp
195                           zsum_alfa = 0._wp
196                           zsum_beta = 0._wp
197                           zsum_z    = 0._wp
198                                                   
199                           DO jk = ikup, ikbot      ! Inside the instable (and overlying neutral) portion of the column
200                              !
201                              zdz       = fse3t(ji,jj,jk)
202                              zsum_temp = zsum_temp + zvts(jk,jp_tem)*zdz
203                              zsum_sali = zsum_sali + zvts(jk,jp_sal)*zdz
204                              zsum_alfa = zsum_alfa + zvab(jk,jp_tem)*zdz
205                              zsum_beta = zsum_beta + zvab(jk,jp_sal)*zdz
206                              zsum_z    = zsum_z    + zdz
207                              !                             
208                              IF( jk == ikbot ) EXIT ! avoid array-index overshoot in case ikbot = jpk, cause we're calling jk+1 next line
209                              !! EXIT when we have reached the last layer that is instable (N2<0) or neutral (N2=0):
210                              IF( zvn2(jk+1) > zn2_zero ) EXIT
211                           END DO
212                         
213                           ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative or neutral N2
214                           IF( ikup == ikdown )   CALL ctl_stop( 'tra_npc :  PROBLEM #2')
215
216                           ! Mixing Temperature, salinity, alpha and beta from ikup to ikdown included:
217                           zta   = zsum_temp/zsum_z
218                           zsa   = zsum_sali/zsum_z
219                           zalfa = zsum_alfa/zsum_z
220                           zbeta = zsum_beta/zsum_z
221
222                           IF( lp_monitor_point ) THEN
223                              WRITE(numout,*) 'MIXED T, S, alfa and beta between ikup =',ikup,   &
224                                 &            ' and ikdown =',ikdown,', in layer #',ilayer
225                              WRITE(numout,*) '  => Mean temp. in that portion =', zta
226                              WRITE(numout,*) '  => Mean sali. in that portion =', zsa
227                              WRITE(numout,*) '  => Mean Alfa  in that portion =', zalfa
228                              WRITE(numout,*) '  => Mean Beta  in that portion =', zbeta
229                              IF(lflush) CALL flush(numout)
230                           ENDIF
231
232                           !! Homogenaizing the temperature, salinity, alpha and beta in this portion of the column
233                           DO jk = ikup, ikdown
234                              zvts(jk,jp_tem) = zta
235                              zvts(jk,jp_sal) = zsa
236                              zvab(jk,jp_tem) = zalfa
237                              zvab(jk,jp_sal) = zbeta
238                           END DO
239                           
240                           
241                           !! Updating N2 in the relvant portion of the water column
242                           !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion
243                           !! => Need to re-compute N2! will use Alpha and Beta!
244                           
245                           ikup   = MAX(2,ikup)         ! ikup can never be 1 !
246                           ik_low = MIN(ikdown+1,ikbot) ! we must go 1 point deeper than ikdown!
247                           
248                           DO jk = ikup, ik_low              ! we must go 1 point deeper than ikdown!
249
250                              !! Interpolating alfa and beta at W point:
251                              zrw =  (fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk)) &
252                                 & / (fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk))
253                              zaw = zvab(jk,jp_tem) * (1._wp - zrw) + zvab(jk-1,jp_tem) * zrw
254                              zbw = zvab(jk,jp_sal) * (1._wp - zrw) + zvab(jk-1,jp_sal) * zrw
255
256                              !! N2 at W point, doing exactly as in eosbn2.F90:
257                              zvn2(jk) = grav*( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     &
258                                 &            - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )  &
259                                 &       / fse3w(ji,jj,jk) * tmask(ji,jj,jk)
260
261                              !! OR, faster  => just considering the vertical gradient of density
262                              !! as only the signa maters...
263                              !zvn2(jk) = (  zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     &
264                              !     &      - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )
265
266                           END DO
267                       
268                           ikp = MIN(ikdown+1,ikbot)
269                           
270
271                        ENDIF  !IF( zvn2(ikp) < 0. )
272
273
274                        IF( ikp == ikbot ) l_bottom_reached = .TRUE.
275                        !
276                     END DO ! DO WHILE ( .NOT. l_bottom_reached )
277
278                     IF( ikp /= ikbot )   CALL ctl_stop( 'tra_npc :  PROBLEM #3')
279                   
280                     ! ******* At this stage ikp == ikbot ! *******
281                   
282                     IF( ilayer > 0 ) THEN      !! least an unstable layer has been found
283                        !
284                        IF( lp_monitor_point ) THEN
285                           WRITE(numout,*)
286                           WRITE(numout,*) 'After ',jiter,' iteration(s), we neutralized ',ilayer,' instable layer(s)'
287                           WRITE(numout,*) '   ==> N2 at i,j=',ji,',',jj,' now looks like this:'
288                           DO jk = 1, klc1
289                              WRITE(numout,*) jk, zvn2(jk)
290                           END DO
291                           WRITE(numout,*)
292                           IF(lflush) CALL flush(numout)
293                        ENDIF
294                        !
295                        ikp    = 1     ! starting again at the surface for the next iteration
296                        ilayer = 0
297                     ENDIF
298                     !
299                     IF( ikp >= ikbot )   l_column_treated = .TRUE.
300                     !
301                  END DO ! DO WHILE ( .NOT. l_column_treated )
302
303                  !! Updating tsa:
304                  tsa(ji,jj,:,jp_tem) = zvts(:,jp_tem)
305                  tsa(ji,jj,:,jp_sal) = zvts(:,jp_sal)
306
307                  !! LB:  Potentially some other global variable beside theta and S can be treated here
308                  !!      like BGC tracers.
309
310                  IF( lp_monitor_point )   WRITE(numout,*)
311
312               ENDIF ! IF( tmask(ji,jj,3) == 1 ) THEN
313
314            END DO ! ji
315         END DO ! jj
316         !
317         IF( l_trdtra ) THEN         ! send the Non penetrative mixing trends for diagnostic
318            z1_r2dt = 1._wp / (2._wp * rdt)
319            ztrdt(:,:,:) = ( tsa(:,:,:,jp_tem) - ztrdt(:,:,:) ) * z1_r2dt
320            ztrds(:,:,:) = ( tsa(:,:,:,jp_sal) - ztrds(:,:,:) ) * z1_r2dt
321            CALL trd_tra( kt, 'TRA', jp_tem, jptra_npc, ztrdt )
322            CALL trd_tra( kt, 'TRA', jp_sal, jptra_npc, ztrds )
323            CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
324         ENDIF
325         !
326         CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ;   CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )
327         !
328         IF( lwp .AND. l_LB_debug ) THEN
329            WRITE(numout,*) 'Exiting tra_npc , kt = ',kt,', => numb. of statically instable water-columns: ', inpcc
330            WRITE(numout,*)
331            IF(lflush) CALL flush(numout)
332         ENDIF
333         !
334         CALL wrk_dealloc(jpi, jpj, jpk, zn2 )
335         CALL wrk_dealloc(jpi, jpj, jpk, 2, zab )
336         CALL wrk_dealloc(jpk, zvn2 )
337         CALL wrk_dealloc(jpk, 2, zvts, zvab )
338         !
339      ENDIF   ! IF( MOD( kt, nn_npc ) == 0 ) THEN
340      !
341      IF( nn_timing == 1 )  CALL timing_stop('tra_npc')
342      !
343   END SUBROUTINE tra_npc
344
345   !!======================================================================
346END MODULE tranpc
Note: See TracBrowser for help on using the repository browser.