Context Navigation

tranpc.F90 @ 8763

Last change on this file since 8763 was 6140, checked in by timgraham, 8 years ago

Merge of branches/2015/dev_merge_2015 back into trunk. Merge excludes NEMOGCM/TOOLS/OBSTOOLS/ for now due to issues with the change of file type. Will sort these manually with further commits.

Branch merged as follows:
In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
Small conflicts due to bug fixes applied to trunk since the dev_merge_2015 was copied. Bug fixes were applied to the branch as well so these were easy to resolve.
Branch committed at this stage

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2015/dev_merge_2015
to merge the branch into the trunk and then commit - no conflicts at this stage.

Property svn:keywords set to Id

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90 @ 8763

Download in other formats: