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/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

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, 9 years ago
Update 2015/dev_r5021_UKMO1_CICE_coupling branch to revision 5442 of the trunk.
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
[5443]	11	!! 3.6 ! 2015-05 (L. Brodeau) new algorithm based on local Brunt-Vaisala freq.
[503]	12	!!----------------------------------------------------------------------
[3]	13
	14	!!----------------------------------------------------------------------
[1537]	15	!! tra_npc : apply the non penetrative convection scheme
[3]	16	!!----------------------------------------------------------------------
[4990]	17	USE oce ! ocean dynamics and active tracers
[3]	18	USE dom_oce ! ocean space and time domain
[4990]	19	USE phycst ! physical constants
[1537]	20	USE zdf_oce ! ocean vertical physics
[4990]	21	USE trd_oce ! ocean active tracer trends
[2715]	22	USE trdtra ! ocean active tracer trends
[4990]	23	USE eosbn2 ! equation of state (eos routine)
	24	!
[3]	25	USE lbclnk ! lateral boundary conditions (or mpp link)
[216]	26	USE in_out_manager ! I/O manager
[2715]	27	USE lib_mpp ! MPP library
[3294]	28	USE wrk_nemo ! Memory Allocation
	29	USE timing ! Timing
[3]	30
	31	IMPLICIT NONE
	32	PRIVATE
	33
[4990]	34	PUBLIC tra_npc ! routine called by step.F90
[3]	35
	36	!! * Substitutions
	37	# include "domzgr_substitute.h90"
[4990]	38	# include "vectopt_loop_substitute.h90"
[3]	39	!!----------------------------------------------------------------------
[4990]	40	!! NEMO/OPA 3.6 , NEMO Consortium (2014)
[5234]	41	!! $Id$
[2528]	42	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	43	!!----------------------------------------------------------------------
	44	CONTAINS
	45
	46	SUBROUTINE tra_npc( kt )
	47	!!----------------------------------------------------------------------
	48	!! * ROUTINE tranpc *
	49	!!
[4990]	50	!! ** Purpose : Non-penetrative convective adjustment scheme. solve
[1111]	51	!! the static instability of the water column on after fields
[3]	52	!! while conserving heat and salt contents.
	53	!!
[4990]	54	!! ** Method : updated algorithm able to deal with non-linear equation of state
	55	!! (i.e. static stability computed locally)
[3]	56	!!
[1111]	57	!! ** Action : - (ta,sa) after the application od the npc scheme
[4990]	58	!! - send the associated trends for on-line diagnostics (l_trdtra=T)
[3]	59	!!
[4990]	60	!! References : Madec, et al., 1991, JPO, 21, 9, 1349-1371.
[503]	61	!!----------------------------------------------------------------------
	62	INTEGER, INTENT(in) :: kt ! ocean time-step index
[2715]	63	!
[503]	64	INTEGER :: ji, jj, jk ! dummy loop indices
	65	INTEGER :: inpcc ! number of statically instable water column
[5443]	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
[5443]	70	REAL(wp), PARAMETER :: zn2_zero = 1.e-14_wp ! acceptance criteria for neutrality (N2==0)
[4990]	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	!
[5443]	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"
[3]	81	!!----------------------------------------------------------------------
[3294]	82	!
	83	IF( nn_timing == 1 ) CALL timing_start('tra_npc')
	84	!
[1537]	85	IF( MOD( kt, nn_npc ) == 0 ) THEN
[4990]	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
[3]	91
[4990]	92	IF( l_trdtra ) THEN !* Save initial after fields
[3294]	93	CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
[2715]	94	ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
	95	ztrds(:,:,:) = tsa(:,:,:,jp_sal)
[216]	96	ENDIF
	97
[4990]	98	IF( l_LB_debug ) THEN
[5443]	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
[4990]	102	klc1 = mbkt(ilc1,jlc1) ! bottom of the ocean for debug point...
	103	ENDIF
[5443]	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)
[4990]	107
	108	inpcc = 0
[3]	109
[4990]	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
[3]	117
[4990]	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:
[5443]	126	lp_monitor_point = (narea == nncpu).AND.lp_monitor_point
[4990]	127	ENDIF !LB debug end
[3]	128
[4990]	129	ikbot = mbkt(ji,jj) ! ikbot: ocean bottom T-level
[5443]	130	ikp = 1 ! because N2 is irrelevant at the surface level (will start at ikp=2)
[4990]	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.
[3]	142
[4990]	143	DO WHILE ( .NOT. l_bottom_reached )
[3]	144
[5443]	145	ikp = ikp + 1
[4990]	146
[5443]	147	!! Testing level ikp for instability
[4990]	148	!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[5443]	149	IF( zvn2(ikp) < -zn2_zero ) THEN ! Instability found!
[3]	150
[5443]	151	ilayer = ilayer + 1 ! yet another instable portion of the water column found....
[3]	152
[5443]	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
[3]	168
[5443]	169	IF( jiter == 1 ) inpcc = inpcc + 1
[4990]	170
[5443]	171	IF( lp_monitor_point ) WRITE(numout, *) 'Negative N2 at ikp =',ikp,' for layer #', ilayer
[4990]	172
[5443]	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
[4990]	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
[5443]	195	DO jk = ikup, ikbot ! Inside the instable (and overlying neutral) portion of the column
[4990]	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
[5443]	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
[4990]	207	END DO
	208
[5443]	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:
[4990]	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
[5443]	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
[4990]	221	WRITE(numout,*) ' => Mean temp. in that portion =', zta
	222	WRITE(numout,*) ' => Mean sali. in that portion =', zsa
[5443]	223	WRITE(numout,*) ' => Mean Alfa in that portion =', zalfa
[4990]	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
[5443]	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.
[503]	270	!
[4990]	271	END DO ! DO WHILE ( .NOT. l_bottom_reached )
	272
[5443]	273	IF( ikp /= ikbot ) CALL ctl_stop( 'tra_npc : PROBLEM #3')
[4990]	274
[5443]	275	! ***** At this stage ikp == ikbot ! *****
[4990]	276
[5443]	277	IF( ilayer > 0 ) THEN !! least an unstable layer has been found
[503]	278	!
[5443]	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:'
[4990]	283	DO jk = 1, klc1
	284	WRITE(numout,*) jk, zvn2(jk)
	285	END DO
[5443]	286	WRITE(numout,*)
[3]	287	ENDIF
[5443]	288	!
	289	ikp = 1 ! starting again at the surface for the next iteration
[4990]	290	ilayer = 0
	291	ENDIF
	292	!
[5443]	293	IF( ikp >= ikbot ) l_column_treated = .TRUE.
[4990]	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
[5443]	301	!! LB: Potentially some other global variable beside theta and S can be treated here
	302	!! like BGC tracers.
[4990]	303
[5443]	304	IF( lp_monitor_point ) WRITE(numout,*)
[4990]	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 )
[3294]	317	CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
[216]	318	ENDIF
[4990]	319	!
[2715]	320	CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )
[4990]	321	!
[5443]	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,*)
[3]	325	ENDIF
[503]	326	!
[4990]	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
[503]	333	!
[3294]	334	IF( nn_timing == 1 ) CALL timing_stop('tra_npc')
	335	!
[3]	336	END SUBROUTINE tra_npc
	337
	338	!!======================================================================
	339	END MODULE tranpc

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

Download in other formats: