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 branches/UKMO/dev_r8183_GC_couple_pkg/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/UKMO/dev_r8183_GC_couple_pkg/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90 @ 8730

Last change on this file since 8730 was 8730, checked in by dancopsey, 6 years ago

Cleared out SVN keywords.

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