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_r11943_MERGE_2019/src/OCE/TRA – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/tranpc.F90 @ 11949

Last change on this file since 11949 was 11949, checked in by acc, 5 years ago
Merge in changes from 2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps. This just creates a fresh copy of this branch to use as the merge base. See ticket #2341
Property svn:keywords set to `Id`
File size: 17.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: