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/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/tranpc.F90 @ 13819

Last change on this file since 13819 was 13819, checked in by hadcv, 3 years ago

#2365: Final fixes and changes

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