Context Navigation

limmsh_2.F90 @ 2375

Last change on this file since 2375 was 2370, checked in by gm, 14 years ago
v3.3beta: ice-ocean stress at kt with VP & EVP (LIM-2 and -3)
Property svn:keywords set to `Id`
File size: 12.0 KB

Rev	Line
[821]	1	MODULE limmsh_2
[3]	2	!!======================================================================
[821]	3	!! * MODULE limmsh_2 *
	4	!! LIM 2.0 ice model : definition of the ice mesh parameters
[3]	5	!!======================================================================
[2370]	6	!! History : - ! 2001-04 (LIM) original code
	7	!! 1.0 ! 2002-08 (C. Ethe, G. Madec) F90, module
	8	!! 3.3 ! 2009-05 (G. Garric, C. Bricaud) addition of the lim2_evp case
	9	!!----------------------------------------------------------------------
[821]	10	#if defined key_lim2
[3]	11	!!----------------------------------------------------------------------
[821]	12	!! 'key_lim2' LIM 2.0sea-ice model
[58]	13	!!----------------------------------------------------------------------
[821]	14	!! lim_msh_2 : definition of the ice mesh
[3]	15	!!----------------------------------------------------------------------
	16	USE phycst
	17	USE dom_oce
[821]	18	USE dom_ice_2
[3]	19	USE lbclnk
[58]	20	USE in_out_manager
[3]	21
	22	IMPLICIT NONE
	23	PRIVATE
	24
[821]	25	PUBLIC lim_msh_2 ! routine called by ice_ini_2.F90
[3]	26
	27	!!----------------------------------------------------------------------
[2287]	28	!! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
[1156]	29	!! $Id$
[2370]	30	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	31	!!----------------------------------------------------------------------
	32	CONTAINS
	33
[821]	34	SUBROUTINE lim_msh_2
[3]	35	!!-------------------------------------------------------------------
[821]	36	!! * ROUTINE lim_msh_2 *
[3]	37	!!
	38	!! ** Purpose : Definition of the charact. of the numerical grid
	39	!!
	40	!! ** Action : - Initialisation of some variables
	41	!! - Definition of some constants linked with the grid
	42	!! - Definition of the metric coef. for the sea/ice
	43	!! - Initialization of the ice masks (tmsk, umsk)
	44	!!
	45	!! ** Refer. : Deleersnijder et al. Ocean Modelling 100, 7-10
	46	!!---------------------------------------------------------------------
[12]	47	INTEGER :: ji, jj ! dummy loop indices
[2319]	48	REAL(wp) :: zusden ! local scalars
	49	#if defined key_lim2_vp
	50	REAL(wp) :: zusden2 ! local scalars
	51	REAL(wp) :: zh1p , zh2p ! - -
	52	REAL(wp) :: zd2d1p, zd1d2p ! - -
	53	REAL(wp), DIMENSION(jpi,jpj) :: zd2d1 , zd1d2 ! 2D workspace
	54	#endif
[3]	55	!!---------------------------------------------------------------------
[58]	56
	57	IF(lwp) THEN
	58	WRITE(numout,*)
[821]	59	WRITE(numout,*) 'lim_msh_2 : LIM 2.0 sea-ice model, mesh initialization'
	60	WRITE(numout,*) '~~~~~~~~~'
[58]	61	ENDIF
[3]	62
[2004]	63	IF( jphgr_msh == 2 .OR. jphgr_msh == 3 .OR. jphgr_msh == 5 ) &
	64	& CALL ctl_stop(' Coriolis parameter in LIM not set for f- or beta-plane' )
	65
[3]	66	!----------------------------------------------------------
	67	! Initialization of local and some global (common) variables
	68	!------------------------------------------------------------------
	69
[58]	70	njeq = INT( jpj / 2 ) !i bug mpp potentiel
	71	njeqm1 = njeq - 1
[3]	72
[2319]	73	fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad ) ! coriolis factor at T-point
[3]	74
[58]	75	!i DO jj = 1, jpj
	76	!i zmsk(jj) = SUM( tmask(:,jj,:) ) ! = 0 if land everywhere on a j-line
	77	!!ii write(numout,*) jj, zind(jj)
	78	!i END DO
[3]	79
[58]	80	IF( fcor(1,1) * fcor(1,nlcj) < 0.e0 ) THEN ! local domain include both hemisphere
	81	l_jeq = .TRUE.
	82	njeq = 1
	83	DO WHILE ( njeq <= jpj .AND. fcor(1,njeq) < 0.e0 )
	84	njeq = njeq + 1
	85	END DO
	86	IF(lwp ) WRITE(numout,*) ' the equator is inside the domain at about njeq = ', njeq
	87	ELSEIF( fcor(1,1) < 0.e0 ) THEN
	88	l_jeq = .FALSE.
[192]	89	njeq = jpj
[58]	90	IF(lwp ) WRITE(numout,*) ' the model domain is entirely in the southern hemisphere: njeq = ', njeq
	91	ELSE
	92	l_jeq = .FALSE.
[192]	93	njeq = 2
[58]	94	IF(lwp ) WRITE(numout,*) ' the model domain is entirely in the northern hemisphere: njeq = ', njeq
	95	ENDIF
	96
	97	njeqm1 = njeq - 1
	98
	99
[3]	100	! For each grid, definition of geometric tables
	101	!------------------------------------------------------------------
	102
	103	!-------------------
[58]	104	! Conventions : !
[3]	105	!-------------------
	106	! indices 1 \ 2 <-> localisation in the 2 direction x \ y
	107	! 3rd indice <-> localisation on the mesh :
	108	! 0 = Centre ; 1 = corner W x(i-1/2) ; 2 = corner S y(j-1/2) ;
	109	! 3 = corner SW x(i-1/2),y(j-1/2)
	110	!-------------------
[58]	111	!!ibug ???
	112	wght(:,:,:,:) = 0.e0
[2319]	113	tmu(:,:) = 0.e0
	114	#if defined key_lim2_vp
	115	akappa(:,:,:,:) = 0.e0
[58]	116	alambd(:,:,:,:,:,:) = 0.e0
[2319]	117	#else
	118	tmv(:,:) = 0.e0
	119	tmf(:,:) = 0.e0
	120	#endif
[58]	121	!!i
[3]	122
[2319]	123
	124	#if defined key_lim2_vp
[3]	125	! metric coefficients for sea ice dynamic
	126	!----------------------------------------
	127	! ! akappa
[12]	128	DO jj = 2, jpj
[58]	129	zd1d2(:,jj) = e1v(:,jj) - e1v(:,jj-1)
[3]	130	END DO
[12]	131	CALL lbc_lnk( zd1d2, 'T', -1. )
[3]	132
[58]	133	DO ji = 2, jpi
	134	zd2d1(ji,:) = e2u(ji,:) - e2u(ji-1,:)
[3]	135	END DO
[12]	136	CALL lbc_lnk( zd2d1, 'T', -1. )
[3]	137
[58]	138	akappa(:,:,1,1) = 1.0 / ( 2.0 * e1t(:,:) )
	139	akappa(:,:,1,2) = zd1d2(:,:) / ( 4.0 * e1t(:,:) * e2t(:,:) )
	140	akappa(:,:,2,1) = zd2d1(:,:) / ( 4.0 * e1t(:,:) * e2t(:,:) )
	141	akappa(:,:,2,2) = 1.0 / ( 2.0 * e2t(:,:) )
[3]	142
	143	! ! weights (wght)
[12]	144	DO jj = 2, jpj
	145	DO ji = 2, jpi
[3]	146	zusden = 1. / ( ( e1t(ji,jj) + e1t(ji-1,jj ) ) &
	147	& * ( e2t(ji,jj) + e2t(ji ,jj-1) ) )
	148	wght(ji,jj,1,1) = zusden * e1t(ji ,jj) * e2t(ji,jj )
	149	wght(ji,jj,1,2) = zusden * e1t(ji ,jj) * e2t(ji,jj-1)
	150	wght(ji,jj,2,1) = zusden * e1t(ji-1,jj) * e2t(ji,jj )
	151	wght(ji,jj,2,2) = zusden * e1t(ji-1,jj) * e2t(ji,jj-1)
	152	END DO
	153	END DO
	154	CALL lbc_lnk( wght(:,:,1,1), 'I', 1. ) ! CAUTION: even with the lbc_lnk at ice U-V-point
	155	CALL lbc_lnk( wght(:,:,1,2), 'I', 1. ) ! the value of wght at jpj is wrong
	156	CALL lbc_lnk( wght(:,:,2,1), 'I', 1. ) ! but it is never used
	157	CALL lbc_lnk( wght(:,:,2,2), 'I', 1. )
[2319]	158	#else
	159	! metric coefficients for sea ice dynamic (EVP rheology)
	160	!----------------------------------------
	161	DO jj = 1, jpjm1 ! weights (wght) at F-points
	162	DO ji = 1, jpim1
	163	zusden = 1. / ( ( e1t(ji+1,jj ) + e1t(ji,jj) ) &
	164	& * ( e2t(ji ,jj+1) + e2t(ji,jj) ) )
	165	wght(ji,jj,1,1) = zusden * e1t(ji+1,jj) * e2t(ji,jj+1)
	166	wght(ji,jj,1,2) = zusden * e1t(ji+1,jj) * e2t(ji,jj )
	167	wght(ji,jj,2,1) = zusden * e1t(ji ,jj) * e2t(ji,jj+1)
	168	wght(ji,jj,2,2) = zusden * e1t(ji ,jj) * e2t(ji,jj )
	169	END DO
	170	END DO
	171	CALL lbc_lnk( wght(:,:,1,1), 'F', 1. ) ; CALL lbc_lnk( wght(:,:,1,2),'F', 1. ) ! lateral boundary cond.
	172	CALL lbc_lnk( wght(:,:,2,1), 'F', 1. ) ; CALL lbc_lnk( wght(:,:,2,2),'F', 1. )
	173	#endif
[3]	174
	175	! Coefficients for divergence of the stress tensor
	176	!-------------------------------------------------
	177
[2319]	178	#if defined key_lim2_vp
[12]	179	DO jj = 2, jpj
[1694]	180	DO ji = 2, jpi ! NO vector opt.
[3]	181	zh1p = e1t(ji ,jj ) * wght(ji,jj,2,2) &
	182	& + e1t(ji-1,jj ) * wght(ji,jj,1,2) &
	183	& + e1t(ji ,jj-1) * wght(ji,jj,2,1) &
	184	& + e1t(ji-1,jj-1) * wght(ji,jj,1,1)
	185
	186	zh2p = e2t(ji ,jj ) * wght(ji,jj,2,2) &
	187	& + e2t(ji-1,jj ) * wght(ji,jj,1,2) &
	188	& + e2t(ji ,jj-1) * wght(ji,jj,2,1) &
	189	& + e2t(ji-1,jj-1) * wght(ji,jj,1,1)
	190
[58]	191	! better written but change the last digit and thus solver in less than 100 timestep
[3]	192	! zh1p = e1t(ji-1,jj ) * wght(ji,jj,1,2) + e1t(ji,jj ) * wght(ji,jj,2,2) &
	193	! & + e1t(ji-1,jj-1) * wght(ji,jj,1,1) + e1t(ji,jj-1) * wght(ji,jj,2,1)
	194
	195	! zh2p = e2t(ji-1,jj ) * wght(ji,jj,1,2) + e2t(ji,jj ) * wght(ji,jj,2,2) &
	196	! & + e2t(ji-1,jj-1) * wght(ji,jj,1,1) + e2t(ji,jj-1) * wght(ji,jj,2,1)
	197
[58]	198	!!ibug =0 zusden = 1.0 / ( zh1p * zh2p * 4.e0 )
	199	zusden = 1.0 / MAX( zh1p * zh2p * 4.e0 , 1.e-20 )
[3]	200	zusden2 = zusden * 2.0
	201
	202	zd1d2p = zusden * 0.5 * ( -e1t(ji-1,jj-1) + e1t(ji-1,jj ) - e1t(ji,jj-1) + e1t(ji ,jj) )
	203	zd2d1p = zusden * 0.5 * ( e2t(ji ,jj-1) - e2t(ji-1,jj-1) + e2t(ji,jj ) - e2t(ji-1,jj) )
	204
	205	alambd(ji,jj,2,2,2,1) = zusden2 * e2t(ji ,jj-1)
	206	alambd(ji,jj,2,2,2,2) = zusden2 * e2t(ji ,jj )
	207	alambd(ji,jj,2,2,1,1) = zusden2 * e2t(ji-1,jj-1)
	208	alambd(ji,jj,2,2,1,2) = zusden2 * e2t(ji-1,jj )
	209
	210	alambd(ji,jj,1,1,2,1) = zusden2 * e1t(ji ,jj-1)
	211	alambd(ji,jj,1,1,2,2) = zusden2 * e1t(ji ,jj )
	212	alambd(ji,jj,1,1,1,1) = zusden2 * e1t(ji-1,jj-1)
	213	alambd(ji,jj,1,1,1,2) = zusden2 * e1t(ji-1,jj )
	214
	215	alambd(ji,jj,1,2,2,1) = zd1d2p
	216	alambd(ji,jj,1,2,2,2) = zd1d2p
	217	alambd(ji,jj,1,2,1,1) = zd1d2p
	218	alambd(ji,jj,1,2,1,2) = zd1d2p
	219
	220	alambd(ji,jj,2,1,2,1) = zd2d1p
	221	alambd(ji,jj,2,1,2,2) = zd2d1p
	222	alambd(ji,jj,2,1,1,1) = zd2d1p
	223	alambd(ji,jj,2,1,1,2) = zd2d1p
	224	END DO
	225	END DO
	226
	227	CALL lbc_lnk( alambd(:,:,2,2,2,1), 'I', 1. ) ! CAUTION: even with the lbc_lnk at ice U-V point
	228	CALL lbc_lnk( alambd(:,:,2,2,2,2), 'I', 1. ) ! the value of wght at jpj is wrong
	229	CALL lbc_lnk( alambd(:,:,2,2,1,1), 'I', 1. ) ! but it is never used
	230	CALL lbc_lnk( alambd(:,:,2,2,1,2), 'I', 1. ) !
	231
	232	CALL lbc_lnk( alambd(:,:,1,1,2,1), 'I', 1. ) ! CAUTION: idem
	233	CALL lbc_lnk( alambd(:,:,1,1,2,2), 'I', 1. ) !
	234	CALL lbc_lnk( alambd(:,:,1,1,1,1), 'I', 1. ) !
	235	CALL lbc_lnk( alambd(:,:,1,1,1,2), 'I', 1. ) !
	236
	237	CALL lbc_lnk( alambd(:,:,1,2,2,1), 'I', 1. ) ! CAUTION: idem
	238	CALL lbc_lnk( alambd(:,:,1,2,2,2), 'I', 1. ) !
	239	CALL lbc_lnk( alambd(:,:,1,2,1,1), 'I', 1. ) !
	240	CALL lbc_lnk( alambd(:,:,1,2,1,2), 'I', 1. ) !
	241
	242	CALL lbc_lnk( alambd(:,:,2,1,2,1), 'I', 1. ) ! CAUTION: idem
	243	CALL lbc_lnk( alambd(:,:,2,1,2,2), 'I', 1. ) !
	244	CALL lbc_lnk( alambd(:,:,2,1,1,1), 'I', 1. ) !
	245	CALL lbc_lnk( alambd(:,:,2,1,1,2), 'I', 1. ) !
[2319]	246	#endif
[3]	247
	248
	249	! Initialization of ice masks
	250	!----------------------------
	251
	252	tms(:,:) = tmask(:,:,1) ! ice T-point : use surface tmask
	253
[2319]	254	#if defined key_lim2_vp
	255	! VP rheology : ice velocity point is I-point
[3]	256	!i here we can use umask with a i and j shift of -1,-1
	257	tmu(:,1) = 0.e0
	258	tmu(1,:) = 0.e0
	259	DO jj = 2, jpj ! ice U.V-point: computed from ice T-point mask
[1694]	260	DO ji = 2, jpim1 ! NO vector opt.
[3]	261	tmu(ji,jj) = tms(ji,jj) * tms(ji-1,jj) * tms(ji,jj-1) * tms(ji-1,jj-1)
	262	END DO
	263	END DO
[2370]	264	CALL lbc_lnk( tmu(:,:), 'I', 1. ) !--lateral boundary conditions
[2319]	265	#else
	266	! EVP rheology : ice velocity point are U- & V-points ; ice vorticity
	267	! point is F-point
	268	tmu(:,:) = umask(:,:,1)
	269	tmv(:,:) = vmask(:,:,1)
	270	tmf(:,:) = 0.e0 ! used of fmask except its special value along the coast (rn_shlat)
	271	WHERE( fmask(:,:,1) == 1.e0 ) tmf(:,:) = 1.e0
	272	#endif
[2370]	273	!
[3]	274	! unmasked and masked area of T-grid cell
	275	area(:,:) = e1t(:,:) * e2t(:,:)
[2370]	276	!
[821]	277	END SUBROUTINE lim_msh_2
[58]	278
[3]	279	#else
[58]	280	!!----------------------------------------------------------------------
	281	!! Default option Dummy Module NO LIM sea-ice model
	282	!!----------------------------------------------------------------------
[3]	283	CONTAINS
[821]	284	SUBROUTINE lim_msh_2 ! Dummy routine
	285	END SUBROUTINE lim_msh_2
[58]	286	#endif
[3]	287
	288	!!======================================================================
[821]	289	END MODULE limmsh_2

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

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

Context Navigation

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90 @ 2375

Download in other formats: