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.

limtrp_2.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_2/limtrp_2.F90 @ 3837

Last change on this file since 3837 was 3558, checked in by rblod, 12 years ago
Fix issues when using key_nosignedzeo, see ticket #996
Property svn:keywords set to `Id`
File size: 16.9 KB

Rev	Line
[821]	1	MODULE limtrp_2
[3]	2	!!======================================================================
[821]	3	!! * MODULE limtrp_2 *
	4	!! LIM 2.0 transport ice model : sea-ice advection/diffusion
[3]	5	!!======================================================================
[1922]	6	!! History : LIM ! 2000-01 (UCL) Original code
	7	!! 2.0 ! 2001-05 (G. Madec, R. Hordoir) opa norm
	8	!! - ! 2004-01 (G. Madec, C. Ethe) F90, mpp
[2528]	9	!! 3.3 ! 2009-05 (G. Garric, C. Bricaud) addition of the lim2_evp case
[1922]	10	!!----------------------------------------------------------------------
[821]	11	#if defined key_lim2
[3]	12	!!----------------------------------------------------------------------
[821]	13	!! 'key_lim2' : LIM 2.0 sea-ice model
[3]	14	!!----------------------------------------------------------------------
[821]	15	!! lim_trp_2 : advection/diffusion process of sea ice
	16	!! lim_trp_init_2 : initialization and namelist read
[3]	17	!!----------------------------------------------------------------------
[1922]	18	USE phycst ! physical constant
	19	USE sbc_oce ! ocean surface boundary condition
	20	USE dom_oce ! ocean domain
[3]	21	USE in_out_manager ! I/O manager
[1922]	22	USE dom_ice_2 ! LIM-2 domain
	23	USE ice_2 ! LIM-2 variables
	24	USE limistate_2 ! LIM-2 initial state
	25	USE limadv_2 ! LIM-2 advection
	26	USE limhdf_2 ! LIM-2 horizontal diffusion
	27	USE lbclnk ! lateral boundary conditions -- MPP exchanges
	28	USE lib_mpp ! MPP library
[3294]	29	USE wrk_nemo ! work arrays
[3558]	30	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[3]	31
	32	IMPLICIT NONE
	33	PRIVATE
	34
[1922]	35	PUBLIC lim_trp_2 ! called by sbc_ice_lim_2
[3]	36
[2528]	37	REAL(wp), PUBLIC :: bound = 0.e0 !: boundary condit. (0.0 no-slip, 1.0 free-slip)
[12]	38
[2528]	39	REAL(wp) :: epsi06 = 1.e-06 ! constant values
	40	REAL(wp) :: epsi03 = 1.e-03
	41	REAL(wp) :: epsi16 = 1.e-16
	42	REAL(wp) :: rzero = 0.e0
	43	REAL(wp) :: rone = 1.e0
[3]	44
	45	!! * Substitution
	46	# include "vectopt_loop_substitute.h90"
	47	!!----------------------------------------------------------------------
[2528]	48	!! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
[1156]	49	!! $Id$
[2528]	50	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	51	!!----------------------------------------------------------------------
	52
	53	CONTAINS
	54
[821]	55	SUBROUTINE lim_trp_2( kt )
[3]	56	!!-------------------------------------------------------------------
[821]	57	!! * ROUTINE lim_trp_2 *
[3]	58	!!
	59	!! ** purpose : advection/diffusion process of sea ice
	60	!!
	61	!! ** method : variables included in the process are scalar,
	62	!! other values are considered as second order.
	63	!! For advection, a second order Prather scheme is used.
	64	!!
	65	!! ** action :
	66	!!---------------------------------------------------------------------
[508]	67	INTEGER, INTENT(in) :: kt ! number of iteration
[1922]	68	!!
	69	INTEGER :: ji, jj, jk ! dummy loop indices
	70	INTEGER :: initad ! number of sub-timestep for the advection
	71	REAL(wp) :: zindb , zindsn , zindic, zacrith ! local scalars
	72	REAL(wp) :: zusvosn, zusvoic, zignm , zindhe ! - -
	73	REAL(wp) :: zvbord , zcfl , zusnit ! - -
	74	REAL(wp) :: zrtt , ztsn , ztic1 , ztic2 ! - -
[3294]	75	REAL(wp), POINTER, DIMENSION(:,:) :: zui_u , zvi_v , zsm ! 2D workspace
	76	REAL(wp), POINTER, DIMENSION(:,:) :: zs0ice, zs0sn , zs0a ! - -
	77	REAL(wp), POINTER, DIMENSION(:,:) :: zs0c0 , zs0c1 , zs0c2 , zs0st ! - -
[3]	78	!---------------------------------------------------------------------
	79
[3294]	80	CALL wrk_alloc( jpi, jpj, zui_u , zvi_v , zsm, zs0ice, zs0sn , zs0a, zs0c0 , zs0c1 , zs0c2 , zs0st )
[2715]	81
[821]	82	IF( kt == nit000 ) CALL lim_trp_init_2 ! Initialization (first time-step only)
[3]	83
	84	zsm(:,:) = area(:,:)
	85
[76]	86	IF( ln_limdyn ) THEN
[3]	87	!-------------------------------------!
	88	! Advection of sea ice properties !
	89	!-------------------------------------!
	90
	91	! ice velocities at ocean U- and V-points (zui_u,zvi_v)
	92	! ---------------------------------------
[2528]	93	IF( lk_lim2_vp ) THEN ! VP rheology : B-grid sea-ice dynamics (I-point ice velocity)
	94	zvbord = 1._wp + ( 1._wp - bound ) ! zvbord=2 no-slip, =0 free slip boundary conditions
	95	DO jj = 1, jpjm1
	96	DO ji = 1, jpim1 ! NO vector opt.
	97	zui_u(ji,jj) = ( u_ice(ji+1,jj) + u_ice(ji+1,jj+1) ) / ( MAX( tmu(ji+1,jj)+tmu(ji+1,jj+1), zvbord ) )
	98	zvi_v(ji,jj) = ( v_ice(ji,jj+1) + v_ice(ji+1,jj+1) ) / ( MAX( tmu(ji,jj+1)+tmu(ji+1,jj+1), zvbord ) )
	99	END DO
[3]	100	END DO
[2528]	101	CALL lbc_lnk( zui_u, 'U', -1. ) ; CALL lbc_lnk( zvi_v, 'V', -1. ) ! Lateral boundary conditions
	102	!
	103	ELSE ! EVP rheology : C-grid sea-ice dynamics (u- & v-points ice velocity)
	104	zui_u(:,:) = u_ice(:,:) ! EVP rheology: ice (u,v) at u- and v-points
	105	zvi_v(:,:) = v_ice(:,:)
	106	ENDIF
[3]	107
	108	! CFL test for stability
	109	! ----------------------
[2528]	110	zcfl = 0._wp
[3]	111	zcfl = MAX( zcfl, MAXVAL( ABS( zui_u(1:jpim1, : ) ) * rdt_ice / e1u(1:jpim1, : ) ) )
	112	zcfl = MAX( zcfl, MAXVAL( ABS( zvi_v( : ,1:jpjm1) ) * rdt_ice / e2v( : ,1:jpjm1) ) )
[1922]	113	!
	114	IF(lk_mpp) CALL mpp_max( zcfl )
	115	!
	116	IF( zcfl > 0.5 .AND. lwp ) WRITE(numout,*) 'lim_trp_2 : violation of cfl criterion the ',nday,'th day, cfl = ', zcfl
[3]	117
	118	! content of properties
	119	! ---------------------
[2528]	120	zs0sn (:,:) = hsnm(:,:) * area (:,:) ! Snow volume.
	121	zs0ice(:,:) = hicm(:,:) * area (:,:) ! Ice volume.
[1922]	122	zs0a (:,:) = ( 1.0 - frld(:,:) ) * area (:,:) ! Surface covered by ice.
	123	zs0c0 (:,:) = tbif(:,:,1) / rt0_snow * zs0sn (:,:) ! Heat content of the snow layer.
[3]	124	zs0c1 (:,:) = tbif(:,:,2) / rt0_ice * zs0ice(:,:) ! Heat content of the first ice layer.
	125	zs0c2 (:,:) = tbif(:,:,3) / rt0_ice * zs0ice(:,:) ! Heat content of the second ice layer.
[1922]	126	zs0st (:,:) = qstoif(:,:) / xlic * zs0a (:,:) ! Heat reservoir for brine pockets.
[3]	127
	128
[1922]	129	! Advection (Prather scheme)
[3]	130	! ---------
[1922]	131	initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) ) ! If ice drift field is too fast,
	132	zusnit = 1.0 / REAL( initad ) ! split the ice time step in two
	133	!
	134	IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0) THEN !== odd ice time step: adv_x then adv_y ==!
	135	DO jk = 1, initad
[821]	136	CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
	137	CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
	138	CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0sn , sxsn , sxxsn , sysn , syysn , sxysn )
	139	CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0sn , sxsn , sxxsn , sysn , syysn , sxysn )
	140	CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0a , sxa , sxxa , sya , syya , sxya )
	141	CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0a , sxa , sxxa , sya , syya , sxya )
	142	CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0 )
	143	CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0 )
	144	CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1 )
	145	CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1 )
	146	CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2 )
	147	CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2 )
	148	CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0st , sxst , sxxst , syst , syyst , sxyst )
	149	CALL lim_adv_y_2( zusnit, zvi_v, rzero, zsm, zs0st , sxst , sxxst , syst , syyst , sxyst )
[3]	150	END DO
[1922]	151	ELSE !== even ice time step: adv_x then adv_y ==!
[3]	152	DO jk = 1, initad
[821]	153	CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
	154	CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
	155	CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0sn , sxsn , sxxsn , sysn , syysn , sxysn )
	156	CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0sn , sxsn , sxxsn , sysn , syysn , sxysn )
	157	CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0a , sxa , sxxa , sya , syya , sxya )
	158	CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0a , sxa , sxxa , sya , syya , sxya )
	159	CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0 )
	160	CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0c0 , sxc0 , sxxc0 , syc0 , syyc0 , sxyc0 )
	161	CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1 )
	162	CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0c1 , sxc1 , sxxc1 , syc1 , syyc1 , sxyc1 )
	163	CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2 )
	164	CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0c2 , sxc2 , sxxc2 , syc2 , syyc2 , sxyc2 )
	165	CALL lim_adv_y_2( zusnit, zvi_v, rone , zsm, zs0st , sxst , sxxst , syst , syyst , sxyst )
	166	CALL lim_adv_x_2( zusnit, zui_u, rzero, zsm, zs0st , sxst , sxxst , syst , syyst , sxyst )
[3]	167	END DO
	168	ENDIF
	169
	170	! recover the properties from their contents
	171	! ------------------------------------------
[1922]	172	!!gm Define in limmsh one for all area = 1 /area (CPU time saved !)
[3]	173	zs0ice(:,:) = zs0ice(:,:) / area(:,:)
	174	zs0sn (:,:) = zs0sn (:,:) / area(:,:)
	175	zs0a (:,:) = zs0a (:,:) / area(:,:)
	176	zs0c0 (:,:) = zs0c0 (:,:) / area(:,:)
	177	zs0c1 (:,:) = zs0c1 (:,:) / area(:,:)
	178	zs0c2 (:,:) = zs0c2 (:,:) / area(:,:)
	179	zs0st (:,:) = zs0st (:,:) / area(:,:)
	180
	181
	182	!-------------------------------------!
	183	! Diffusion of sea ice properties !
	184	!-------------------------------------!
	185
	186	! Masked eddy diffusivity coefficient at ocean U- and V-points
	187	! ------------------------------------------------------------
	188	DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row
	189	DO ji = 1 , fs_jpim1 ! vector opt.
	190	pahu(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji ,jj) ) ) ) &
	191	& * ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji+1,jj) ) ) ) * ahiu(ji,jj)
	192	pahv(ji,jj) = ( 1.0 - MAX( rzero, SIGN( rone, -zs0a(ji,jj ) ) ) ) &
	193	& * ( 1.0 - MAX( rzero, SIGN( rone,- zs0a(ji,jj+1) ) ) ) * ahiv(ji,jj)
	194	END DO
	195	END DO
[1922]	196	!!gm more readable coding: (and avoid an error in F90 with sign of zero)
	197	! DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row
	198	! DO ji = 1 , fs_jpim1 ! vector opt.
[2528]	199	! IF( MIN( zs0a(ji,jj) , zs0a(ji+1,jj) ) == 0.e0 ) pahu(ji,jj) = 0._wp
	200	! IF( MIN( zs0a(ji,jj) , zs0a(ji,jj+1) ) == 0.e0 ) pahv(ji,jj) = 0._wp
[1922]	201	! END DO
	202	! END DO
	203	!!gm end
[3]	204
	205	! diffusion
	206	! ---------
[821]	207	CALL lim_hdf_2( zs0ice )
	208	CALL lim_hdf_2( zs0sn )
	209	CALL lim_hdf_2( zs0a )
	210	CALL lim_hdf_2( zs0c0 )
	211	CALL lim_hdf_2( zs0c1 )
	212	CALL lim_hdf_2( zs0c2 )
	213	CALL lim_hdf_2( zs0st )
[3]	214
[1922]	215	!!gm see comment this can be skipped
	216	zs0ice(:,:) = MAX( rzero, zs0ice(:,:) * area(:,:) ) !!bug: useless
	217	zs0sn (:,:) = MAX( rzero, zs0sn (:,:) * area(:,:) ) !!bug: cf /area just below
	218	zs0a (:,:) = MAX( rzero, zs0a (:,:) * area(:,:) ) !! caution: the suppression of the 2 changes
	219	zs0c0 (:,:) = MAX( rzero, zs0c0 (:,:) * area(:,:) ) !! the last digit of the results
[3]	220	zs0c1 (:,:) = MAX( rzero, zs0c1 (:,:) * area(:,:) )
	221	zs0c2 (:,:) = MAX( rzero, zs0c2 (:,:) * area(:,:) )
	222	zs0st (:,:) = MAX( rzero, zs0st (:,:) * area(:,:) )
	223
	224
[1922]	225	!-------------------------------------------------------------------!
	226	! Updating and limitation of sea ice properties after transport !
	227	!-------------------------------------------------------------------!
[3]	228	DO jj = 1, jpj
[76]	229	zindhe = MAX( 0.e0, SIGN( 1.e0, fcor(1,jj) ) ) ! = 0 for SH, =1 for NH
[3]	230	DO ji = 1, jpi
[1922]	231	!
[3]	232	! Recover mean values over the grid squares.
	233	zs0sn (ji,jj) = MAX( rzero, zs0sn (ji,jj)/area(ji,jj) )
	234	zs0ice(ji,jj) = MAX( rzero, zs0ice(ji,jj)/area(ji,jj) )
	235	zs0a (ji,jj) = MAX( rzero, zs0a (ji,jj)/area(ji,jj) )
	236	zs0c0 (ji,jj) = MAX( rzero, zs0c0 (ji,jj)/area(ji,jj) )
	237	zs0c1 (ji,jj) = MAX( rzero, zs0c1 (ji,jj)/area(ji,jj) )
	238	zs0c2 (ji,jj) = MAX( rzero, zs0c2 (ji,jj)/area(ji,jj) )
	239	zs0st (ji,jj) = MAX( rzero, zs0st (ji,jj)/area(ji,jj) )
	240
	241	! Recover in situ values.
	242	zindb = MAX( rzero, SIGN( rone, zs0a(ji,jj) - epsi06 ) )
	243	zacrith = 1.0 - ( zindhe * acrit(1) + ( 1.0 - zindhe ) * acrit(2) )
	244	zs0a (ji,jj) = zindb * MIN( zs0a(ji,jj), zacrith )
	245	hsnif(ji,jj) = zindb * ( zs0sn(ji,jj) /MAX( zs0a(ji,jj), epsi16 ) )
	246	hicif(ji,jj) = zindb * ( zs0ice(ji,jj)/MAX( zs0a(ji,jj), epsi16 ) )
	247	zindsn = MAX( rzero, SIGN( rone, hsnif(ji,jj) - epsi06 ) )
	248	zindic = MAX( rzero, SIGN( rone, hicif(ji,jj) - epsi03 ) )
	249	zindb = MAX( zindsn, zindic )
	250	zs0a (ji,jj) = zindb * zs0a(ji,jj)
	251	frld (ji,jj) = 1.0 - zs0a(ji,jj)
	252	hsnif(ji,jj) = zindsn * hsnif(ji,jj)
	253	hicif(ji,jj) = zindic * hicif(ji,jj)
	254	zusvosn = 1.0/MAX( hsnif(ji,jj) * zs0a(ji,jj), epsi16 )
	255	zusvoic = 1.0/MAX( hicif(ji,jj) * zs0a(ji,jj), epsi16 )
	256	zignm = MAX( rzero, SIGN( rone, hsndif - hsnif(ji,jj) ) )
	257	zrtt = 173.15 * rone
	258	ztsn = zignm * tbif(ji,jj,1) &
	259	+ ( 1.0 - zignm ) * MIN( MAX( zrtt, rt0_snow * zusvosn * zs0c0(ji,jj)) , tfu(ji,jj) )
	260	ztic1 = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c1(ji,jj) ) , tfu(ji,jj) )
	261	ztic2 = MIN( MAX( zrtt, rt0_ice * zusvoic * zs0c2(ji,jj) ) , tfu(ji,jj) )
	262
	263	tbif(ji,jj,1) = zindsn * ztsn + ( 1.0 - zindsn ) * tfu(ji,jj)
	264	tbif(ji,jj,2) = zindic * ztic1 + ( 1.0 - zindic ) * tfu(ji,jj)
	265	tbif(ji,jj,3) = zindic * ztic2 + ( 1.0 - zindic ) * tfu(ji,jj)
	266	qstoif(ji,jj) = zindb * xlic * zs0st(ji,jj) / MAX( zs0a(ji,jj), epsi16 )
	267	END DO
	268	END DO
[1922]	269	!
[3]	270	ENDIF
[1922]	271	!
[3294]	272	CALL wrk_dealloc( jpi, jpj, zui_u , zvi_v , zsm, zs0ice, zs0sn , zs0a, zs0c0 , zs0c1 , zs0c2 , zs0st )
[2715]	273	!
[821]	274	END SUBROUTINE lim_trp_2
[3]	275
	276
[821]	277	SUBROUTINE lim_trp_init_2
[3]	278	!!-------------------------------------------------------------------
[821]	279	!! * ROUTINE lim_trp_init_2 *
[3]	280	!!
	281	!! ** Purpose : initialization of ice advection parameters
	282	!!
[1922]	283	!! ** Method : Read the namicetrp namelist and check the parameter
	284	!! values called at the first timestep (nit000)
[3]	285	!!
	286	!! ** input : Namelist namicetrp
	287	!!-------------------------------------------------------------------
	288	NAMELIST/namicetrp/ bound
	289	!!-------------------------------------------------------------------
[1922]	290	!
	291	REWIND ( numnam_ice ) ! Read Namelist namicetrp
[3]	292	READ ( numnam_ice , namicetrp )
	293	IF(lwp) THEN
	294	WRITE(numout,*)
[821]	295	WRITE(numout,*) 'lim_trp_init_2 : Ice parameters for advection '
	296	WRITE(numout,*) '~~~~~~~~~~~~~~'
[76]	297	WRITE(numout,*) ' boundary conditions (0. no-slip, 1. free-slip) bound = ', bound
[3]	298	ENDIF
[1922]	299	!
[821]	300	END SUBROUTINE lim_trp_init_2
[3]	301
	302	#else
	303	!!----------------------------------------------------------------------
	304	!! Default option Empty Module No sea-ice model
	305	!!----------------------------------------------------------------------
	306	CONTAINS
[821]	307	SUBROUTINE lim_trp_2 ! Empty routine
	308	END SUBROUTINE lim_trp_2
[3]	309	#endif
	310
	311	!!======================================================================
[821]	312	END MODULE limtrp_2

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

Download in other formats: