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 @ 5167

Last change on this file since 5167 was 5167, checked in by clem, 9 years ago
LIM3 bug fix. see ticket #1492 (forthcoming update) which also solve ticket #1497
Property svn:keywords set to `Id`
File size: 13.5 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
[2715]	10	!!----------------------------------------------------------------------
[834]	11	#if defined key_lim3
	12	!!----------------------------------------------------------------------
[3625]	13	!! 'key_lim3' LIM-3 sea-ice model
[834]	14	!!----------------------------------------------------------------------
[3625]	15	!! lim_cons : checks whether energy, mass and salt are conserved
[825]	16	!!----------------------------------------------------------------------
[4688]	17	USE phycst ! physical constants
[3625]	18	USE ice ! LIM-3 variables
	19	USE dom_ice ! LIM-3 domain
	20	USE dom_oce ! ocean domain
	21	USE in_out_manager ! I/O manager
	22	USE lib_mpp ! MPP library
	23	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[5167]	24	USE sbc_oce , ONLY : sfx ! Surface boundary condition: ocean fields
[825]	25
	26	IMPLICIT NONE
	27	PRIVATE
	28
[2715]	29	PUBLIC lim_column_sum
	30	PUBLIC lim_column_sum_energy
	31	PUBLIC lim_cons_check
[4688]	32	PUBLIC lim_cons_hsm
[5167]	33	PUBLIC lim_cons_final
[825]	34
	35	!!----------------------------------------------------------------------
[4161]	36	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
[1156]	37	!! $Id$
[2715]	38	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	39	!!----------------------------------------------------------------------
	40	CONTAINS
	41
[2715]	42	SUBROUTINE lim_column_sum( ksum, pin, pout )
	43	!!-------------------------------------------------------------------
	44	!! * ROUTINE lim_column_sum *
	45	!!
	46	!! ** Purpose : Compute the sum of xin over nsum categories
	47	!!
	48	!! ** Method : Arithmetics
	49	!!
	50	!! ** Action : Gets xin(ji,jj,jl) and computes xout(ji,jj)
	51	!!---------------------------------------------------------------------
	52	INTEGER , INTENT(in ) :: ksum ! number of categories/layers
	53	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pin ! input field
	54	REAL(wp), DIMENSION(:,:) , INTENT( out) :: pout ! output field
	55	!
	56	INTEGER :: jl ! dummy loop indices
	57	!!---------------------------------------------------------------------
	58	!
	59	pout(:,:) = pin(:,:,1)
	60	DO jl = 2, ksum
	61	pout(:,:) = pout(:,:) + pin(:,:,jl)
	62	END DO
	63	!
[825]	64	END SUBROUTINE lim_column_sum
	65
	66
[2715]	67	SUBROUTINE lim_column_sum_energy( ksum, klay, pin, pout)
[825]	68	!!-------------------------------------------------------------------
	69	!! * ROUTINE lim_column_sum_energy *
	70	!!
	71	!! ** Purpose : Compute the sum of xin over nsum categories
	72	!! and nlay layers
	73	!!
	74	!! ** Method : Arithmetics
	75	!!---------------------------------------------------------------------
[4873]	76	INTEGER , INTENT(in ) :: ksum !: number of categories
	77	INTEGER , INTENT(in ) :: klay !: number of vertical layers
	78	REAL(wp), DIMENSION(jpi,jpj,nlay_i+1,jpl), INTENT(in ) :: pin !: input field
	79	REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pout !: output field
[2715]	80	!
	81	INTEGER :: jk, jl ! dummy loop indices
[825]	82	!!---------------------------------------------------------------------
[2715]	83	!
[2777]	84	pout(:,:) = 0._wp
[2715]	85	DO jl = 1, ksum
	86	DO jk = 2, klay
	87	pout(:,:) = pout(:,:) + pin(:,:,jk,jl)
	88	END DO
	89	END DO
	90	!
[825]	91	END SUBROUTINE lim_column_sum_energy
	92
[921]	93
[2715]	94	SUBROUTINE lim_cons_check( px1, px2, pmax_err, cd_fieldid )
[825]	95	!!-------------------------------------------------------------------
	96	!! * ROUTINE lim_cons_check *
	97	!!
	98	!! ** Purpose : Test the conservation of a certain variable
	99	!! For each physical grid cell, check that initial
	100	!! and final values
	101	!! of a conserved field are equal to within a small value.
	102	!!
	103	!! ** Method :
	104	!!---------------------------------------------------------------------
[2715]	105	REAL(wp), DIMENSION(:,:), INTENT(in ) :: px1 !: initial field
	106	REAL(wp), DIMENSION(:,:), INTENT(in ) :: px2 !: final field
	107	REAL(wp) , INTENT(in ) :: pmax_err !: max allowed error
	108	CHARACTER(len=15) , INTENT(in ) :: cd_fieldid !: field identifyer
	109	!
	110	INTEGER :: ji, jj ! dummy loop indices
	111	INTEGER :: inb_error ! number of g.c where there is a cons. error
	112	LOGICAL :: llconserv_err ! = .true. if conservation check failed
	113	REAL(wp) :: zmean_error ! mean error on error points
[825]	114	!!---------------------------------------------------------------------
[2715]	115	!
	116	IF(lwp) WRITE(numout,*) ' lim_cons_check '
	117	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
[825]	118
[2715]	119	llconserv_err = .FALSE.
	120	inb_error = 0
	121	zmean_error = 0._wp
	122	IF( MAXVAL( px2(:,:) - px1(:,:) ) > pmax_err ) llconserv_err = .TRUE.
[825]	123
[2715]	124	IF( llconserv_err ) THEN
[825]	125	DO jj = 1, jpj
	126	DO ji = 1, jpi
[2715]	127	IF( ABS( px2(ji,jj) - px1(ji,jj) ) > pmax_err ) THEN
	128	inb_error = inb_error + 1
	129	zmean_error = zmean_error + ABS( px2(ji,jj) - px1(ji,jj) )
	130	!
	131	IF(lwp) THEN
	132	WRITE (numout,*) ' ALERTE 99 '
	133	WRITE (numout,*) ' Conservation error: ', cd_fieldid
	134	WRITE (numout,*) ' Point : ', ji, jj
	135	WRITE (numout,*) ' lat, lon : ', gphit(ji,jj), glamt(ji,jj)
	136	WRITE (numout,*) ' Initial value : ', px1(ji,jj)
	137	WRITE (numout,*) ' Final value : ', px2(ji,jj)
	138	WRITE (numout,*) ' Difference : ', px2(ji,jj) - px1(ji,jj)
	139	ENDIF
[825]	140	ENDIF
	141	END DO
	142	END DO
[2715]	143	!
	144	ENDIF
	145	IF(lk_mpp) CALL mpp_sum( inb_error )
	146	IF(lk_mpp) CALL mpp_sum( zmean_error )
	147	!
	148	IF( inb_error > 0 .AND. lwp ) THEN
	149	zmean_error = zmean_error / REAL( inb_error, wp )
	150	WRITE(numout,*) ' Conservation check for : ', cd_fieldid
	151	WRITE(numout,*) ' Number of error points : ', inb_error
	152	WRITE(numout,*) ' Mean error on these pts: ', zmean_error
	153	ENDIF
	154	!
[825]	155	END SUBROUTINE lim_cons_check
	156
[4688]	157
	158	SUBROUTINE lim_cons_hsm( icount, cd_routine, zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b )
	159	!!-------------------------------------------------------------------
	160	!! * ROUTINE lim_cons_hsm *
	161	!!
	162	!! ** Purpose : Test the conservation of heat, salt and mass for each routine
	163	!!
	164	!! ** Method :
	165	!!---------------------------------------------------------------------
	166	INTEGER , INTENT(in) :: icount ! determine wether this is the beggining of the routine (0) or the end (1)
	167	CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine
	168	REAL(wp) , INTENT(inout) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
	169	REAL(wp) :: zvi, zsmv, zei, zfs, zfw, zft
	170	REAL(wp) :: zvmin, zamin, zamax
[5167]	171	REAL(wp) :: zvtrp, zetrp
	172	REAL(wp), PARAMETER :: zconv = 1.e-9
[4688]	173
	174	IF( icount == 0 ) THEN
	175
[5123]	176	zfs_b = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + &
	177	& sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) &
	178	& ) * e12t(:,:) * tmask(:,:,1) )
[4688]	179
[5123]	180	zfw_b = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + &
	181	& wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) &
	182	& ) * e12t(:,:) * tmask(:,:,1) )
	183
	184	zft_b = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) &
	185	& - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) &
	186	& ) * e12t(:,:) * tmask(:,:,1) * zconv )
	187
	188	zvi_b = glob_sum( SUM( v_i(:,:,:)rhoic + v_s(:,:,:)rhosn, dim=3 ) * e12t(:,:) * tmask(:,:,1) )
	189
[5167]	190	zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) * rhoic )
[5123]	191
	192	zei_b = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) + &
	193	& SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) &
	194	) * e12t(:,:) * tmask(:,:,1) * zconv )
	195
[4688]	196	ELSEIF( icount == 1 ) THEN
	197
[5123]	198	zfs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + &
	199	& sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) &
	200	& ) * e12t(:,:) * tmask(:,:,1) ) - zfs_b
	201
	202	zfw = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + &
	203	& wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) &
	204	& ) * e12t(:,:) * tmask(:,:,1) ) - zfw_b
	205
	206	zft = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) &
	207	& - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) &
	208	& ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zft_b
[4688]	209
[5123]	210	zvi = ( glob_sum( SUM( v_i(:,:,:)rhoic + v_s(:,:,:)rhosn, dim=3 ) &
	211	& * e12t(:,:) * tmask(:,:,1) ) - zvi_b ) * r1_rdtice - zfw
[4688]	212
[5167]	213	zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) * rhoic ) - zsmv_b ) * r1_rdtice + zfs
[5123]	214
	215	zei = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) + &
	216	& SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) &
	217	& ) * e12t(:,:) * tmask(:,:,1) * zconv ) * r1_rdtice - zei_b * r1_rdtice + zft
	218
[5167]	219	zvtrp = glob_sum( ( diag_trp_vi * rhoic + diag_trp_vs * rhosn ) * e12t(:,:) * tmask(:,:,1) )
	220	zetrp = glob_sum( ( diag_trp_ei + diag_trp_es ) * e12t(:,:) * tmask(:,:,1) * zconv )
[5123]	221	zvmin = glob_min( v_i )
	222	zamax = glob_max( SUM( a_i, dim=3 ) )
	223	zamin = glob_min( a_i )
[5167]	224
[4688]	225
	226	IF(lwp) THEN
[5167]	227	IF ( ABS( zvi * rday ) > 0.5 * 1.e9 ) WRITE(numout,) 'violation volume [kg/day] (',cd_routine,') = ',(zvi rday)
	228	IF ( ABS( zsmv * rday ) > 5. * 1.e9 ) WRITE(numout,) 'violation saline [psukg/day] (',cd_routine,') = ',(zsmv * rday)
	229	IF ( ABS( zei ) > 2. * 1.e9 ) WRITE(numout,*) 'violation enthalpy [GW] (',cd_routine,') = ',(zei)
	230	IF ( zvmin < -epsi10 ) WRITE(numout,*) 'violation v_i<0 [m] (',cd_routine,') = ',(zvmin)
[5123]	231	IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > rn_amax+epsi10 ) THEN
[5167]	232	WRITE(numout,*) 'violation a_i>amax (',cd_routine,') = ',zamax
[4688]	233	ENDIF
[5167]	234	IF ( zamin < -epsi10 ) WRITE(numout,*) 'violation a_i<0 (',cd_routine,') = ',zamin
	235	IF( cd_routine == 'limtrp' .AND. ABS( zvtrp * rday ) > 0.5*1.e9 ) THEN
	236	WRITE(numout,) 'violation vtrp [kg/day] (',cd_routine,') = ',(zvtrp rday)
	237	WRITE(numout,*) 'violation etrp [GW] (',cd_routine,') = ',(zetrp )
	238	ENDIF
[4688]	239	ENDIF
	240
	241	ENDIF
	242
	243	END SUBROUTINE lim_cons_hsm
	244
[5167]	245	SUBROUTINE lim_cons_final( cd_routine )
	246	CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine
	247	REAL(wp) :: zhfx, zsfx, zvfx
	248	REAL(wp), PARAMETER :: zconv = 1.e-9
	249
	250	zhfx = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es - hfx_sub ) * e12t(:,:) * tmask(:,:,1) * zconv )
	251	zsfx = glob_sum( ( sfx + diag_smvi ) * e12t(:,:) * tmask(:,:,1) ) * rday
	252	zvfx = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e12t(:,:) * tmask(:,:,1) ) * rday
	253
	254	! if error > 1 mm / 100 years over the Arctic Basin
	255	IF( ABS( zvfx ) > 0.5 * 1.e9 ) WRITE(numout,*) 'violation vfx [kg/day] (',cd_routine,') = ',(zvfx)
	256	! if error > 1 mm / 100 years over the Arctic Basin (ice with latent heat = 3e6 J/kg)
	257	IF( ABS( zhfx ) > 2. * 1.e9 ) WRITE(numout,*) 'violation hfx [GW] (',cd_routine,') = ',(zhfx)
	258	! if error > 1 mm / 100 years over the Arctic Basin (ice of salinity = 10 pss)
	259	IF( ABS( zsfx ) > 5. * 1.e9 ) WRITE(numout,) 'violation sfx [psukg/day] (',cd_routine,') = ',(zsfx)
	260
	261	END SUBROUTINE lim_cons_final
	262
[834]	263	#else
	264	!!----------------------------------------------------------------------
	265	!! Default option Empty module NO LIM sea-ice model
	266	!!----------------------------------------------------------------------
	267	#endif
	268	!!======================================================================
	269	END MODULE limcons

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

Download in other formats: