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_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tranpc.F90 @ 10954

Last change on this file since 10954 was 10954, checked in by acc, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert TRA modules and all knock on effects of these conversions. SETTE tested

  • Property svn:keywords set to Id
File size: 17.2 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
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)
[4990]24   !
[6140]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
[3]29
30   IMPLICIT NONE
31   PRIVATE
32
[4990]33   PUBLIC   tra_npc    ! routine called by step.F90
[3]34
35   !! * Substitutions
[4990]36#  include "vectopt_loop_substitute.h90"
[3]37   !!----------------------------------------------------------------------
[9598]38   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[4990]39   !! $Id$
[10068]40   !! Software governed by the CeCILL license (see ./LICENSE)
[3]41   !!----------------------------------------------------------------------
42CONTAINS
43
[10946]44   SUBROUTINE tra_npc( kt, Kmm, Krhs )
[3]45      !!----------------------------------------------------------------------
46      !!                  ***  ROUTINE tranpc  ***
47      !!
[4990]48      !! ** Purpose : Non-penetrative convective adjustment scheme. solve
[1111]49      !!      the static instability of the water column on after fields
[3]50      !!      while conserving heat and salt contents.
51      !!
[4990]52      !! ** Method  : updated algorithm able to deal with non-linear equation of state
53      !!              (i.e. static stability computed locally)
[3]54      !!
[6140]55      !! ** Action  : - tsa: after tracers with the application of the npc scheme
[4990]56      !!              - send the associated trends for on-line diagnostics (l_trdtra=T)
[3]57      !!
[4990]58      !! References :     Madec, et al., 1991, JPO, 21, 9, 1349-1371.
[503]59      !!----------------------------------------------------------------------
60      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
[10946]61      INTEGER, INTENT(in) ::   Kmm, Krhs  ! time level indices
[2715]62      !
[503]63      INTEGER  ::   ji, jj, jk   ! dummy loop indices
64      INTEGER  ::   inpcc        ! number of statically instable water column
[5386]65      INTEGER  ::   jiter, ikbot, ikp, ikup, ikdown, ilayer, ik_low   ! local integers
[4990]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
[9019]69      REAL(wp), PARAMETER ::   zn2_zero = 1.e-14_wp      ! acceptance criteria for neutrality (N2==0)
70      REAL(wp), DIMENSION(        jpk     ) ::   zvn2         ! vertical profile of N2 at 1 given point...
71      REAL(wp), DIMENSION(        jpk,jpts) ::   zvts, zvab   ! vertical profile of T & S , and  alpha & betaat 1 given point
72      REAL(wp), DIMENSION(jpi,jpj,jpk     ) ::   zn2          ! N^2
73      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::   zab          ! alpha and beta
74      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds   ! 3D workspace
[4990]75      !
[5386]76      LOGICAL, PARAMETER :: l_LB_debug = .FALSE. ! set to true if you want to follow what is
77      INTEGER :: ilc1, jlc1, klc1, nncpu         ! actually happening in a water column at point "ilc1, jlc1"
78      LOGICAL :: lp_monitor_point = .FALSE.      ! in CPU domain "nncpu"
[3]79      !!----------------------------------------------------------------------
[3294]80      !
[9019]81      IF( ln_timing )   CALL timing_start('tra_npc')
[3294]82      !
[1537]83      IF( MOD( kt, nn_npc ) == 0 ) THEN
[4990]84         !
85         IF( l_trdtra )   THEN                    !* Save initial after fields
[9019]86            ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
[10954]87            ztrdt(:,:,:) = ts(:,:,:,jp_tem,Krhs) 
88            ztrds(:,:,:) = ts(:,:,:,jp_sal,Krhs)
[216]89         ENDIF
[9019]90         !
[4990]91         IF( l_LB_debug ) THEN
[5386]92            ! Location of 1 known convection site to follow what's happening in the water column
93            ilc1 = 45 ;  jlc1 = 3 ; !  ORCA2 4x4, Antarctic coast, more than 2 unstable portions in the  water column...           
94            nncpu = 1  ;            ! the CPU domain contains the convection spot
[4990]95            klc1 =  mbkt(ilc1,jlc1) ! bottom of the ocean for debug point...         
96         ENDIF
[9019]97         !
[10954]98         CALL eos_rab( ts(:,:,:,:,Krhs), zab, Kmm )         ! after alpha and beta (given on T-points)
99         CALL bn2    ( ts(:,:,:,:,Krhs), zab, zn2, Kmm )    ! after Brunt-Vaisala  (given on W-points)
[9019]100         !
[4990]101         inpcc = 0
[9019]102         !
[4990]103         DO jj = 2, jpjm1                 ! interior column only
104            DO ji = fs_2, fs_jpim1
105               !
106               IF( tmask(ji,jj,2) == 1 ) THEN      ! At least 2 ocean points
107                  !                                     ! consider one ocean column
[10954]108                  zvts(:,jp_tem) = ts(ji,jj,:,jp_tem,Krhs)      ! temperature
109                  zvts(:,jp_sal) = ts(ji,jj,:,jp_sal,Krhs)      ! salinity
[6140]110                  !
[4990]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
[6140]114                  !
[4990]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:
[5386]119                     lp_monitor_point = (narea == nncpu).AND.lp_monitor_point                     
[4990]120                  ENDIF                                  !LB debug  end
[6140]121                  !
[4990]122                  ikbot = mbkt(ji,jj)   ! ikbot: ocean bottom T-level
[5386]123                  ikp = 1                  ! because N2 is irrelevant at the surface level (will start at ikp=2)
[4990]124                  ilayer = 0
125                  jiter  = 0
126                  l_column_treated = .FALSE.
[6140]127                  !
[4990]128                  DO WHILE ( .NOT. l_column_treated )
129                     !
130                     jiter = jiter + 1
[6140]131                     !
[4990]132                     IF( jiter >= 400 ) EXIT
[6140]133                     !
[4990]134                     l_bottom_reached = .FALSE.
[6140]135                     !
[4990]136                     DO WHILE ( .NOT. l_bottom_reached )
[6140]137                        !
[5386]138                        ikp = ikp + 1
[6140]139                        !
[5386]140                        !! Testing level ikp for instability
[4990]141                        !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[5386]142                        IF( zvn2(ikp) <  -zn2_zero ) THEN ! Instability found!
[6140]143                           !
[5386]144                           ilayer = ilayer + 1    ! yet another instable portion of the water column found....
[6140]145                           !
[5386]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
[6140]160                           !
[5386]161                           IF( jiter == 1 )   inpcc = inpcc + 1 
[6140]162                           !
[5386]163                           IF( lp_monitor_point )   WRITE(numout, *) 'Negative N2 at ikp =',ikp,' for layer #', ilayer
[6140]164                           !
[5386]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
[6140]167                           !
[5386]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
[6140]178                           !
[5386]179                           IF( ikup < 1 )   CALL ctl_stop( 'tra_npc :  PROBLEM #1')
[6140]180                           !
[4990]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                                                   
[5386]187                           DO jk = ikup, ikbot      ! Inside the instable (and overlying neutral) portion of the column
[4990]188                              !
[10954]189                              zdz       = e3t(ji,jj,jk,Kmm)
[4990]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
[5386]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
[4990]199                           END DO
200                         
[5386]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:
[4990]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
[5386]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
[4990]213                              WRITE(numout,*) '  => Mean temp. in that portion =', zta
214                              WRITE(numout,*) '  => Mean sali. in that portion =', zsa
[5386]215                              WRITE(numout,*) '  => Mean Alfa  in that portion =', zalfa
[4990]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
[5386]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:
[10954]238                              zrw =  (gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm)) &
239                                 & / (gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm))
[5386]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) )  )  &
[10954]246                                 &       / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
[5386]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.
[503]262                        !
[4990]263                     END DO ! DO WHILE ( .NOT. l_bottom_reached )
264
[5386]265                     IF( ikp /= ikbot )   CALL ctl_stop( 'tra_npc :  PROBLEM #3')
[4990]266                   
[5386]267                     ! ******* At this stage ikp == ikbot ! *******
[4990]268                   
[5386]269                     IF( ilayer > 0 ) THEN      !! least an unstable layer has been found
[503]270                        !
[5386]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:'
[4990]275                           DO jk = 1, klc1
276                              WRITE(numout,*) jk, zvn2(jk)
277                           END DO
[5386]278                           WRITE(numout,*)
[3]279                        ENDIF
[5386]280                        !
281                        ikp    = 1     ! starting again at the surface for the next iteration
[4990]282                        ilayer = 0
283                     ENDIF
284                     !
[5386]285                     IF( ikp >= ikbot )   l_column_treated = .TRUE.
[4990]286                     !
287                  END DO ! DO WHILE ( .NOT. l_column_treated )
288
289                  !! Updating tsa:
[10954]290                  ts(ji,jj,:,jp_tem,Krhs) = zvts(:,jp_tem)
291                  ts(ji,jj,:,jp_sal,Krhs) = zvts(:,jp_sal)
[4990]292
[5386]293                  !! LB:  Potentially some other global variable beside theta and S can be treated here
294                  !!      like BGC tracers.
[4990]295
[5386]296                  IF( lp_monitor_point )   WRITE(numout,*)
[4990]297
298               ENDIF ! IF( tmask(ji,jj,3) == 1 ) THEN
299
300            END DO ! ji
301         END DO ! jj
302         !
303         IF( l_trdtra ) THEN         ! send the Non penetrative mixing trends for diagnostic
304            z1_r2dt = 1._wp / (2._wp * rdt)
[10954]305            ztrdt(:,:,:) = ( ts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) ) * z1_r2dt
306            ztrds(:,:,:) = ( ts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:) ) * z1_r2dt
[10946]307            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_npc, ztrdt )
308            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_npc, ztrds )
[9019]309            DEALLOCATE( ztrdt, ztrds )
[216]310         ENDIF
[4990]311         !
[10954]312         CALL lbc_lnk_multi( 'tranpc', ts(:,:,:,jp_tem,Krhs), 'T', 1., ts(:,:,:,jp_sal,Krhs), 'T', 1. )
[4990]313         !
[5386]314         IF( lwp .AND. l_LB_debug ) THEN
315            WRITE(numout,*) 'Exiting tra_npc , kt = ',kt,', => numb. of statically instable water-columns: ', inpcc
316            WRITE(numout,*)
[3]317         ENDIF
[503]318         !
[4990]319      ENDIF   ! IF( MOD( kt, nn_npc ) == 0 ) THEN
[503]320      !
[9019]321      IF( ln_timing )   CALL timing_stop('tra_npc')
[3294]322      !
[3]323   END SUBROUTINE tra_npc
324
325   !!======================================================================
326END MODULE tranpc
Note: See TracBrowser for help on using the repository browser.