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.

limcons.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90 @ 5176

Last change on this file since 5176 was 5176, checked in by clem, 9 years ago
LIM3 online diagnostics (for debugging). Improvements of the online check of heat, salt and mass. The thresholds used to print errors on conservation (in ocean.output) now depend on the ice area (and therefore are independent on the simulation.
Property svn:keywords set to `Id`
File size: 16.2 KB

Rev	Line
[825]	1	MODULE limcons
[2715]	2	!!======================================================================
	3	!! * MODULE limcons *
	4	!! LIM-3 Sea Ice : conservation check
	5	!!======================================================================
	6	!! History : - ! Original code from William H. Lipscomb, LANL
	7	!! 3.0 ! 2004-06 (M. Vancoppenolle) Energy Conservation
[5123]	8	!! 3.5 ! 2011-02 (G. Madec) add mpp considerations
[4688]	9	!! - ! 2014-05 (C. Rousset) add lim_cons_hsm
[5176]	10	!! - ! 2015-03 (C. Rousset) add lim_cons_final
[2715]	11	!!----------------------------------------------------------------------
[834]	12	#if defined key_lim3
	13	!!----------------------------------------------------------------------
[3625]	14	!! 'key_lim3' LIM-3 sea-ice model
[834]	15	!!----------------------------------------------------------------------
[3625]	16	!! lim_cons : checks whether energy, mass and salt are conserved
[825]	17	!!----------------------------------------------------------------------
[4688]	18	USE phycst ! physical constants
[3625]	19	USE ice ! LIM-3 variables
	20	USE dom_ice ! LIM-3 domain
	21	USE dom_oce ! ocean domain
	22	USE in_out_manager ! I/O manager
	23	USE lib_mpp ! MPP library
	24	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[5167]	25	USE sbc_oce , ONLY : sfx ! Surface boundary condition: ocean fields
[825]	26
	27	IMPLICIT NONE
	28	PRIVATE
	29
[2715]	30	PUBLIC lim_column_sum
	31	PUBLIC lim_column_sum_energy
	32	PUBLIC lim_cons_check
[4688]	33	PUBLIC lim_cons_hsm
[5167]	34	PUBLIC lim_cons_final
[825]	35
	36	!!----------------------------------------------------------------------
[4161]	37	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
[1156]	38	!! $Id$
[2715]	39	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	40	!!----------------------------------------------------------------------
	41	CONTAINS
	42
[2715]	43	SUBROUTINE lim_column_sum( ksum, pin, pout )
	44	!!-------------------------------------------------------------------
	45	!! * ROUTINE lim_column_sum *
	46	!!
	47	!! ** Purpose : Compute the sum of xin over nsum categories
	48	!!
	49	!! ** Method : Arithmetics
	50	!!
	51	!! ** Action : Gets xin(ji,jj,jl) and computes xout(ji,jj)
	52	!!---------------------------------------------------------------------
	53	INTEGER , INTENT(in ) :: ksum ! number of categories/layers
	54	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pin ! input field
	55	REAL(wp), DIMENSION(:,:) , INTENT( out) :: pout ! output field
	56	!
	57	INTEGER :: jl ! dummy loop indices
	58	!!---------------------------------------------------------------------
	59	!
	60	pout(:,:) = pin(:,:,1)
	61	DO jl = 2, ksum
	62	pout(:,:) = pout(:,:) + pin(:,:,jl)
	63	END DO
	64	!
[825]	65	END SUBROUTINE lim_column_sum
	66
	67
[2715]	68	SUBROUTINE lim_column_sum_energy( ksum, klay, pin, pout)
[825]	69	!!-------------------------------------------------------------------
	70	!! * ROUTINE lim_column_sum_energy *
	71	!!
	72	!! ** Purpose : Compute the sum of xin over nsum categories
	73	!! and nlay layers
	74	!!
	75	!! ** Method : Arithmetics
	76	!!---------------------------------------------------------------------
[4873]	77	INTEGER , INTENT(in ) :: ksum !: number of categories
	78	INTEGER , INTENT(in ) :: klay !: number of vertical layers
	79	REAL(wp), DIMENSION(jpi,jpj,nlay_i+1,jpl), INTENT(in ) :: pin !: input field
	80	REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pout !: output field
[2715]	81	!
	82	INTEGER :: jk, jl ! dummy loop indices
[825]	83	!!---------------------------------------------------------------------
[2715]	84	!
[2777]	85	pout(:,:) = 0._wp
[2715]	86	DO jl = 1, ksum
	87	DO jk = 2, klay
	88	pout(:,:) = pout(:,:) + pin(:,:,jk,jl)
	89	END DO
	90	END DO
	91	!
[825]	92	END SUBROUTINE lim_column_sum_energy
	93
[921]	94
[2715]	95	SUBROUTINE lim_cons_check( px1, px2, pmax_err, cd_fieldid )
[825]	96	!!-------------------------------------------------------------------
	97	!! * ROUTINE lim_cons_check *
	98	!!
	99	!! ** Purpose : Test the conservation of a certain variable
	100	!! For each physical grid cell, check that initial
	101	!! and final values
	102	!! of a conserved field are equal to within a small value.
	103	!!
	104	!! ** Method :
	105	!!---------------------------------------------------------------------
[2715]	106	REAL(wp), DIMENSION(:,:), INTENT(in ) :: px1 !: initial field
	107	REAL(wp), DIMENSION(:,:), INTENT(in ) :: px2 !: final field
	108	REAL(wp) , INTENT(in ) :: pmax_err !: max allowed error
	109	CHARACTER(len=15) , INTENT(in ) :: cd_fieldid !: field identifyer
	110	!
	111	INTEGER :: ji, jj ! dummy loop indices
	112	INTEGER :: inb_error ! number of g.c where there is a cons. error
	113	LOGICAL :: llconserv_err ! = .true. if conservation check failed
	114	REAL(wp) :: zmean_error ! mean error on error points
[825]	115	!!---------------------------------------------------------------------
[2715]	116	!
	117	IF(lwp) WRITE(numout,*) ' lim_cons_check '
	118	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
[825]	119
[2715]	120	llconserv_err = .FALSE.
	121	inb_error = 0
	122	zmean_error = 0._wp
	123	IF( MAXVAL( px2(:,:) - px1(:,:) ) > pmax_err ) llconserv_err = .TRUE.
[825]	124
[2715]	125	IF( llconserv_err ) THEN
[825]	126	DO jj = 1, jpj
	127	DO ji = 1, jpi
[2715]	128	IF( ABS( px2(ji,jj) - px1(ji,jj) ) > pmax_err ) THEN
	129	inb_error = inb_error + 1
	130	zmean_error = zmean_error + ABS( px2(ji,jj) - px1(ji,jj) )
	131	!
	132	IF(lwp) THEN
	133	WRITE (numout,*) ' ALERTE 99 '
	134	WRITE (numout,*) ' Conservation error: ', cd_fieldid
	135	WRITE (numout,*) ' Point : ', ji, jj
	136	WRITE (numout,*) ' lat, lon : ', gphit(ji,jj), glamt(ji,jj)
	137	WRITE (numout,*) ' Initial value : ', px1(ji,jj)
	138	WRITE (numout,*) ' Final value : ', px2(ji,jj)
	139	WRITE (numout,*) ' Difference : ', px2(ji,jj) - px1(ji,jj)
	140	ENDIF
[825]	141	ENDIF
	142	END DO
	143	END DO
[2715]	144	!
	145	ENDIF
	146	IF(lk_mpp) CALL mpp_sum( inb_error )
	147	IF(lk_mpp) CALL mpp_sum( zmean_error )
	148	!
	149	IF( inb_error > 0 .AND. lwp ) THEN
	150	zmean_error = zmean_error / REAL( inb_error, wp )
	151	WRITE(numout,*) ' Conservation check for : ', cd_fieldid
	152	WRITE(numout,*) ' Number of error points : ', inb_error
	153	WRITE(numout,*) ' Mean error on these pts: ', zmean_error
	154	ENDIF
	155	!
[825]	156	END SUBROUTINE lim_cons_check
	157
[4688]	158
	159	SUBROUTINE lim_cons_hsm( icount, cd_routine, zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b )
[5176]	160	!!--------------------------------------------------------------------------------------------------------
	161	!! * ROUTINE lim_cons_hsm *
[4688]	162	!!
[5176]	163	!! ** Purpose : Test the conservation of heat, salt and mass for each ice routine
	164	!! + test if ice concentration and volume are > 0
[4688]	165	!!
[5176]	166	!! ** Method : This is an online diagnostics which can be activated with ln_limdiahsb=true
	167	!! It prints in ocean.output if there is a violation of conservation at each time-step
	168	!! The thresholds (zv_sill, zs_sill, zh_sill) which determine violations are set to
	169	!! a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years.
	170	!! For salt and heat thresholds, ice is considered to have a salinity of 10
	171	!! and a heat content of 3e5 J/kg (=latent heat of fusion)
	172	!!--------------------------------------------------------------------------------------------------------
	173	INTEGER , INTENT(in) :: icount ! determine wether this is the beggining of the routine (0) or the end (1)
	174	CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine
[4688]	175	REAL(wp) , INTENT(inout) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
	176	REAL(wp) :: zvi, zsmv, zei, zfs, zfw, zft
	177	REAL(wp) :: zvmin, zamin, zamax
[5167]	178	REAL(wp) :: zvtrp, zetrp
[5176]	179	REAL(wp) :: zarea, zv_sill, zs_sill, zh_sill
	180	REAL(wp), PARAMETER :: zconv = 1.e-9 ! convert W to GW and kg to Mt
[4688]	181
	182	IF( icount == 0 ) THEN
	183
[5176]	184	! salt flux
[5123]	185	zfs_b = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + &
	186	& sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) &
[5176]	187	& ) * e12t(:,:) * tmask(:,:,1) * zconv )
[4688]	188
[5176]	189	! water flux
[5123]	190	zfw_b = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + &
	191	& wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) &
[5176]	192	& ) * e12t(:,:) * tmask(:,:,1) * zconv )
[5123]	193
[5176]	194	! heat flux
[5123]	195	zft_b = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) &
	196	& - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) &
	197	& ) * e12t(:,:) * tmask(:,:,1) * zconv )
	198
[5176]	199	zvi_b = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e12t * tmask(:,:,1) * zconv )
[5123]	200
[5176]	201	zsmv_b = glob_sum( SUM( smv_i * rhoic , dim=3 ) * e12t * tmask(:,:,1) * zconv )
[5123]	202
	203	zei_b = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) + &
	204	& SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) &
[5176]	205	) * e12t * tmask(:,:,1) * zconv )
[5123]	206
[4688]	207	ELSEIF( icount == 1 ) THEN
	208
[5176]	209	! salt flux
[5123]	210	zfs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + &
	211	& sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) &
[5176]	212	& ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zfs_b
[5123]	213
[5176]	214	! water flux
[5123]	215	zfw = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + &
	216	& wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) &
[5176]	217	& ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zfw_b
[5123]	218
[5176]	219	! heat flux
[5123]	220	zft = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) &
	221	& - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) &
	222	& ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zft_b
[4688]	223
[5176]	224	! outputs
	225	zvi = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) &
	226	& * e12t * tmask(:,:,1) * zconv ) - zvi_b ) * r1_rdtice - zfw ) * rday
[4688]	227
[5176]	228	zsmv = ( ( glob_sum( SUM( smv_i * rhoic , dim=3 ) &
	229	& * e12t * tmask(:,:,1) * zconv ) - zsmv_b ) * r1_rdtice + zfs ) * rday
[5123]	230
	231	zei = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) + &
	232	& SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) &
[5176]	233	& ) * e12t * tmask(:,:,1) * zconv ) * r1_rdtice - zei_b * r1_rdtice + zft
[5123]	234
[5176]	235	! zvtrp and zetrp must be close to 0 if the advection scheme is conservative
	236	zvtrp = glob_sum( ( diag_trp_vi * rhoic + diag_trp_vs * rhosn ) * e12t * tmask(:,:,1) * zconv ) * rday
	237	zetrp = glob_sum( ( diag_trp_ei + diag_trp_es ) * e12t * tmask(:,:,1) * zconv )
	238
[5123]	239	zvmin = glob_min( v_i )
	240	zamax = glob_max( SUM( a_i, dim=3 ) )
	241	zamin = glob_min( a_i )
[5167]	242
[5176]	243	! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)
	244	zarea = glob_sum( SUM( a_i + epsi10, dim=3 ) * e12t * zconv ) ! in 1.e9 m2
	245	zv_sill = zarea * 2.5e-5
	246	zs_sill = zarea * 25.e-5
	247	zh_sill = zarea * 10.e-5
	248
[4688]	249	IF(lwp) THEN
[5176]	250	IF ( ABS( zvi ) > zv_sill ) WRITE(numout,*) 'violation volume [Mt/day] (',cd_routine,') = ',zvi
	251	IF ( ABS( zsmv ) > zs_sill ) WRITE(numout,) 'violation saline [psuMt/day] (',cd_routine,') = ',zsmv
	252	IF ( ABS( zei ) > zh_sill ) WRITE(numout,*) 'violation enthalpy [GW] (',cd_routine,') = ',zei
	253	IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'limtrp' ) THEN
	254	WRITE(numout,*) 'violation vtrp [Mt/day] (',cd_routine,') = ',zvtrp
	255	WRITE(numout,*) 'violation etrp [GW] (',cd_routine,') = ',zetrp
[4688]	256	ENDIF
[5176]	257	IF ( zvmin < -epsi10 ) WRITE(numout,*) 'violation v_i<0 [m] (',cd_routine,') = ',zvmin
	258	IF ( zamax > rn_amax+epsi10 .AND. cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' ) THEN
	259	WRITE(numout,*) 'violation a_i>amax (',cd_routine,') = ',zamax
[5167]	260	ENDIF
[5176]	261	IF ( zamin < -epsi10 ) WRITE(numout,*) 'violation a_i<0 (',cd_routine,') = ',zamin
[4688]	262	ENDIF
	263
	264	ENDIF
	265
	266	END SUBROUTINE lim_cons_hsm
	267
[5167]	268	SUBROUTINE lim_cons_final( cd_routine )
[5176]	269	!!---------------------------------------------------------------------------------------------------------
	270	!! * ROUTINE lim_cons_final *
	271	!!
	272	!! ** Purpose : Test the conservation of heat, salt and mass at the end of each ice time-step
	273	!!
	274	!! ** Method : This is an online diagnostics which can be activated with ln_limdiahsb=true
	275	!! It prints in ocean.output if there is a violation of conservation at each time-step
	276	!! The thresholds (zv_sill, zs_sill, zh_sill) which determine the violation are set to
	277	!! a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years.
	278	!! For salt and heat thresholds, ice is considered to have a salinity of 10
	279	!! and a heat content of 3e5 J/kg (=latent heat of fusion)
	280	!!--------------------------------------------------------------------------------------------------------
	281	CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine
[5167]	282	REAL(wp) :: zhfx, zsfx, zvfx
[5176]	283	REAL(wp) :: zarea, zv_sill, zs_sill, zh_sill
	284	REAL(wp), PARAMETER :: zconv = 1.e-9 ! convert W to GW and kg to Mt
[5167]	285
[5176]	286	! heat flux
	287	zhfx = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es - hfx_sub ) * e12t * tmask(:,:,1) * zconv )
	288	! salt flux
	289	zsfx = glob_sum( ( sfx + diag_smvi ) * e12t * tmask(:,:,1) * zconv ) * rday
	290	! water flux
	291	zvfx = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e12t * tmask(:,:,1) * zconv ) * rday
[5167]	292
[5176]	293	! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)
	294	zarea = glob_sum( SUM( a_i + epsi10, dim=3 ) * e12t * zconv ) ! in 1.e9 m2
	295	zv_sill = zarea * 2.5e-5
	296	zs_sill = zarea * 25.e-5
	297	zh_sill = zarea * 10.e-5
[5167]	298
[5176]	299	IF( ABS( zvfx ) > zv_sill ) WRITE(numout,*) 'violation vfx [Mt/day] (',cd_routine,') = ',(zvfx)
	300	IF( ABS( zsfx ) > zs_sill ) WRITE(numout,) 'violation sfx [psuMt/day] (',cd_routine,') = ',(zsfx)
	301	IF( ABS( zhfx ) > zh_sill ) WRITE(numout,*) 'violation hfx [GW] (',cd_routine,') = ',(zhfx)
	302
[5167]	303	END SUBROUTINE lim_cons_final
	304
[834]	305	#else
	306	!!----------------------------------------------------------------------
	307	!! Default option Empty module NO LIM sea-ice model
	308	!!----------------------------------------------------------------------
	309	#endif
	310	!!======================================================================
	311	END MODULE limcons

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

Download in other formats: