source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90 @ 5443

Last change on this file since 5443 was 5443, checked in by davestorkey, 5 years ago

Update 2015/dev_r5021_UKMO1_CICE_coupling branch to revision 5442 of the trunk.

File size: 17.4 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                           ENDIF
167                           
168
169                           IF( jiter == 1 )   inpcc = inpcc + 1 
170
171                           IF( lp_monitor_point )   WRITE(numout, *) 'Negative N2 at ikp =',ikp,' for layer #', ilayer
172
173                           !! ikup is the uppermost point where mixing will start:
174                           ikup = ikp - 1 ! ikup is always "at most at ikp-1", less if neutral levels overlying
175                           
176                           !! If the points above ikp-1 have N2 == 0 they must also be mixed:
177                           IF( ikp > 2 ) THEN
178                              DO jk = ikp-1, 2, -1
179                                 IF( ABS(zvn2(jk)) < zn2_zero ) THEN
180                                    ikup = ikup - 1  ! 1 more upper level has N2=0 and must be added for the mixing
181                                 ELSE
182                                    EXIT
183                                 ENDIF
184                              END DO
185                           ENDIF
186                           
187                           IF( ikup < 1 )   CALL ctl_stop( 'tra_npc :  PROBLEM #1')
188
189                           zsum_temp = 0._wp
190                           zsum_sali = 0._wp
191                           zsum_alfa = 0._wp
192                           zsum_beta = 0._wp
193                           zsum_z    = 0._wp
194                                                   
195                           DO jk = ikup, ikbot      ! Inside the instable (and overlying neutral) portion of the column
196                              !
197                              zdz       = fse3t(ji,jj,jk)
198                              zsum_temp = zsum_temp + zvts(jk,jp_tem)*zdz
199                              zsum_sali = zsum_sali + zvts(jk,jp_sal)*zdz
200                              zsum_alfa = zsum_alfa + zvab(jk,jp_tem)*zdz
201                              zsum_beta = zsum_beta + zvab(jk,jp_sal)*zdz
202                              zsum_z    = zsum_z    + zdz
203                              !                             
204                              IF( jk == ikbot ) EXIT ! avoid array-index overshoot in case ikbot = jpk, cause we're calling jk+1 next line
205                              !! EXIT when we have reached the last layer that is instable (N2<0) or neutral (N2=0):
206                              IF( zvn2(jk+1) > zn2_zero ) EXIT
207                           END DO
208                         
209                           ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative or neutral N2
210                           IF( ikup == ikdown )   CALL ctl_stop( 'tra_npc :  PROBLEM #2')
211
212                           ! Mixing Temperature, salinity, alpha and beta from ikup to ikdown included:
213                           zta   = zsum_temp/zsum_z
214                           zsa   = zsum_sali/zsum_z
215                           zalfa = zsum_alfa/zsum_z
216                           zbeta = zsum_beta/zsum_z
217
218                           IF( lp_monitor_point ) THEN
219                              WRITE(numout,*) 'MIXED T, S, alfa and beta between ikup =',ikup,   &
220                                 &            ' and ikdown =',ikdown,', in layer #',ilayer
221                              WRITE(numout,*) '  => Mean temp. in that portion =', zta
222                              WRITE(numout,*) '  => Mean sali. in that portion =', zsa
223                              WRITE(numout,*) '  => Mean Alfa  in that portion =', zalfa
224                              WRITE(numout,*) '  => Mean Beta  in that portion =', zbeta
225                           ENDIF
226
227                           !! Homogenaizing the temperature, salinity, alpha and beta in this portion of the column
228                           DO jk = ikup, ikdown
229                              zvts(jk,jp_tem) = zta
230                              zvts(jk,jp_sal) = zsa
231                              zvab(jk,jp_tem) = zalfa
232                              zvab(jk,jp_sal) = zbeta
233                           END DO
234                           
235                           
236                           !! Updating N2 in the relvant portion of the water column
237                           !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion
238                           !! => Need to re-compute N2! will use Alpha and Beta!
239                           
240                           ikup   = MAX(2,ikup)         ! ikup can never be 1 !
241                           ik_low = MIN(ikdown+1,ikbot) ! we must go 1 point deeper than ikdown!
242                           
243                           DO jk = ikup, ik_low              ! we must go 1 point deeper than ikdown!
244
245                              !! Interpolating alfa and beta at W point:
246                              zrw =  (fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk)) &
247                                 & / (fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk))
248                              zaw = zvab(jk,jp_tem) * (1._wp - zrw) + zvab(jk-1,jp_tem) * zrw
249                              zbw = zvab(jk,jp_sal) * (1._wp - zrw) + zvab(jk-1,jp_sal) * zrw
250
251                              !! N2 at W point, doing exactly as in eosbn2.F90:
252                              zvn2(jk) = grav*( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     &
253                                 &            - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )  &
254                                 &       / fse3w(ji,jj,jk) * tmask(ji,jj,jk)
255
256                              !! OR, faster  => just considering the vertical gradient of density
257                              !! as only the signa maters...
258                              !zvn2(jk) = (  zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     &
259                              !     &      - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )
260
261                           END DO
262                       
263                           ikp = MIN(ikdown+1,ikbot)
264                           
265
266                        ENDIF  !IF( zvn2(ikp) < 0. )
267
268
269                        IF( ikp == ikbot ) l_bottom_reached = .TRUE.
270                        !
271                     END DO ! DO WHILE ( .NOT. l_bottom_reached )
272
273                     IF( ikp /= ikbot )   CALL ctl_stop( 'tra_npc :  PROBLEM #3')
274                   
275                     ! ******* At this stage ikp == ikbot ! *******
276                   
277                     IF( ilayer > 0 ) THEN      !! least an unstable layer has been found
278                        !
279                        IF( lp_monitor_point ) THEN
280                           WRITE(numout,*)
281                           WRITE(numout,*) 'After ',jiter,' iteration(s), we neutralized ',ilayer,' instable layer(s)'
282                           WRITE(numout,*) '   ==> N2 at i,j=',ji,',',jj,' now looks like this:'
283                           DO jk = 1, klc1
284                              WRITE(numout,*) jk, zvn2(jk)
285                           END DO
286                           WRITE(numout,*)
287                        ENDIF
288                        !
289                        ikp    = 1     ! starting again at the surface for the next iteration
290                        ilayer = 0
291                     ENDIF
292                     !
293                     IF( ikp >= ikbot )   l_column_treated = .TRUE.
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                  !! LB:  Potentially some other global variable beside theta and S can be treated here
302                  !!      like BGC tracers.
303
304                  IF( 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
315            CALL trd_tra( kt, 'TRA', jp_tem, jptra_npc, ztrdt )
316            CALL trd_tra( kt, 'TRA', jp_sal, jptra_npc, ztrds )
317            CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
318         ENDIF
319         !
320         CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ;   CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )
321         !
322         IF( lwp .AND. l_LB_debug ) THEN
323            WRITE(numout,*) 'Exiting tra_npc , kt = ',kt,', => numb. of statically instable water-columns: ', inpcc
324            WRITE(numout,*)
325         ENDIF
326         !
327         CALL wrk_dealloc(jpi, jpj, jpk, zn2 )
328         CALL wrk_dealloc(jpi, jpj, jpk, 2, zab )
329         CALL wrk_dealloc(jpk, zvn2 )
330         CALL wrk_dealloc(jpk, 2, zvts, zvab )
331         !
332      ENDIF   ! IF( MOD( kt, nn_npc ) == 0 ) THEN
333      !
334      IF( nn_timing == 1 )  CALL timing_stop('tra_npc')
335      !
336   END SUBROUTINE tra_npc
337
338   !!======================================================================
339END MODULE tranpc
Note: See TracBrowser for help on using the repository browser.