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.

limdyn_2.F90 in branches/dev_001_SBC/NEMO/LIM_SRC – NEMO

Context Navigation

source: branches/dev_001_SBC/NEMO/LIM_SRC/limdyn_2.F90 @ 882

Last change on this file since 882 was 882, checked in by ctlod, 16 years ago
dev_001_SBC: Step II: change effectively the modules names in addining _2 extension, see ticket: #110
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 11.4 KB

Rev	Line
[881]	1	MODULE limdyn_2
[3]	2	!!======================================================================
[881]	3	!! * MODULE limdyn_2 *
[3]	4	!! Sea-Ice dynamics :
	5	!!======================================================================
[717]	6	!! History : 1.0 ! 01-04 (LIM) Original code
	7	!! 2.0 ! 02-08 (C. Ethe, G. Madec) F90, mpp
	8	!! 2.0 ! 03-08 (C. Ethe) add lim_dyn_init
	9	!! 2.0 ! 06-07 (G. Madec) Surface module
	10	!!---------------------------------------------------------------------
[881]	11	#if defined key_lim2
[3]	12	!!----------------------------------------------------------------------
[881]	13	!! 'key_lim2' : LIM 2.0 sea-ice model
[3]	14	!!----------------------------------------------------------------------
[881]	15	!! lim_dyn_2 : computes ice velocities
	16	!! lim_dyn_init_2 : initialization and namelist read
[3]	17	!!----------------------------------------------------------------------
[717]	18	USE dom_oce ! ocean space and time domain
	19	USE sbc_oce !
	20	USE phycst !
[881]	21	USE ice_2 !
[717]	22	USE ice_oce !
[881]	23	USE dom_ice_2 !
	24	USE iceini_2 !
	25	USE limistate_2 !
	26	USE limrhg_2 ! ice rheology
[3]	27
[717]	28	USE lbclnk !
	29	USE lib_mpp !
	30	USE in_out_manager ! I/O manager
	31	USE prtctl ! Print control
	32
[3]	33	IMPLICIT NONE
	34	PRIVATE
	35
[881]	36	PUBLIC lim_dyn_2 ! routine called by sbc_ice_lim
[3]	37
[76]	38	REAL(wp) :: rone = 1.e0 ! constant value
[3]	39
[716]	40	# include "vectopt_loop_substitute.h90"
[3]	41	!!----------------------------------------------------------------------
[717]	42	!! LIM 2.0, UCL-LOCEAN-IPSL (2006)
[699]	43	!! $Id$
[717]	44	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[3]	45	!!----------------------------------------------------------------------
	46
	47	CONTAINS
	48
[881]	49	SUBROUTINE lim_dyn_2( kt )
[3]	50	!!-------------------------------------------------------------------
[881]	51	!! * ROUTINE lim_dyn_2 *
[3]	52	!!
[879]	53	!! ** Purpose : compute ice velocity and ocean-ice friction velocity
[3]	54	!!
	55	!! ** Method :
	56	!!
	57	!! ** Action : - Initialisation
	58	!! - Call of the dynamic routine for each hemisphere
[879]	59	!! - computation of the friction velocity at the sea-ice base
[3]	60	!! - treatment of the case if no ice dynamic
	61	!!---------------------------------------------------------------------
[508]	62	INTEGER, INTENT(in) :: kt ! number of iteration
[717]	63	!!
	64	INTEGER :: ji, jj ! dummy loop indices
	65	INTEGER :: i_j1, i_jpj ! Starting/ending j-indices for rheology
	66	REAL(wp) :: zcoef ! temporary scalar
	67	REAL(wp), DIMENSION(jpj) :: zind ! i-averaged indicator of sea-ice
	68	REAL(wp), DIMENSION(jpj) :: zmsk ! i-averaged of tmask
	69	REAL(wp), DIMENSION(jpi,jpj) :: zu_io, zv_io ! ice-ocean velocity
[3]	70	!!---------------------------------------------------------------------
	71
[881]	72	IF( kt == nit000 ) CALL lim_dyn_init_2 ! Initialization (first time-step only)
[3]	73
[717]	74	IF( ln_limdyn ) THEN
[879]	75	!
[3]	76	! Mean ice and snow thicknesses.
	77	hsnm(:,:) = ( 1.0 - frld(:,:) ) * hsnif(:,:)
	78	hicm(:,:) = ( 1.0 - frld(:,:) ) * hicif(:,:)
[879]	79	!
	80	! ! Rheology (ice dynamics)
	81	! ! ========
[76]	82
	83	! Define the j-limits where ice rheology is computed
	84	! ---------------------------------------------------
	85
[531]	86	IF( lk_mpp .OR. nbit_cmp == 1 ) THEN ! mpp: compute over the whole domain
[76]	87	i_j1 = 1
	88	i_jpj = jpj
[717]	89	IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn : i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
[881]	90	CALL lim_rhg_2( i_j1, i_jpj )
[879]	91	!
[76]	92	ELSE ! optimization of the computational area
[879]	93	!
[76]	94	DO jj = 1, jpj
	95	zind(jj) = SUM( frld (:,jj ) ) ! = FLOAT(jpj) if ocean everywhere on a j-line
	96	zmsk(jj) = SUM( tmask(:,jj,1) ) ! = 0 if land everywhere on a j-line
	97	END DO
[879]	98	!
[76]	99	IF( l_jeq ) THEN ! local domain include both hemisphere
	100	! ! Rheology is computed in each hemisphere
	101	! ! only over the ice cover latitude strip
	102	! Northern hemisphere
	103	i_j1 = njeq
	104	i_jpj = jpj
	105	DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
	106	i_j1 = i_j1 + 1
	107	END DO
	108	i_j1 = MAX( 1, i_j1-1 )
[258]	109	IF(ln_ctl) WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
[879]	110	!
[881]	111	CALL lim_rhg_2( i_j1, i_jpj )
[879]	112	!
[76]	113	! Southern hemisphere
	114	i_j1 = 1
	115	i_jpj = njeq
	116	DO WHILE ( i_jpj >= 1 .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
	117	i_jpj = i_jpj - 1
	118	END DO
	119	i_jpj = MIN( jpj, i_jpj+2 )
[258]	120	IF(ln_ctl) WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
[879]	121	!
[881]	122	CALL lim_rhg_2( i_j1, i_jpj )
[879]	123	!
[76]	124	ELSE ! local domain extends over one hemisphere only
	125	! ! Rheology is computed only over the ice cover
	126	! ! latitude strip
	127	i_j1 = 1
	128	DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
	129	i_j1 = i_j1 + 1
	130	END DO
	131	i_j1 = MAX( 1, i_j1-1 )
	132
	133	i_jpj = jpj
	134	DO WHILE ( i_jpj >= 1 .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
	135	i_jpj = i_jpj - 1
	136	END DO
	137	i_jpj = MIN( jpj, i_jpj+2)
	138
[258]	139	IF(ln_ctl) WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
[879]	140	!
[881]	141	CALL lim_rhg_2( i_j1, i_jpj )
[879]	142	!
[76]	143	ENDIF
[879]	144	!
[76]	145	ENDIF
	146
[717]	147	IF(ln_ctl) CALL prt_ctl(tab2d_1=ui_ice , clinfo1=' lim_dyn : ui_ice :', tab2d_2=vi_ice , clinfo2=' vi_ice :')
[3]	148
[717]	149	! computation of friction velocity
	150	! --------------------------------
	151	! ice-ocean velocity at U & V-points (ui_ice vi_ice at I-point ; ssu_m, ssv_m at U- & V-points)
	152
	153	DO jj = 1, jpjm1
	154	DO ji = 1, fs_jpim1 ! vector opt.
	155	zu_io(ji,jj) = 0.5 * ( ui_ice(ji+1,jj+1) + ui_ice(ji+1,jj ) ) - ssu_m(ji,jj)
	156	zv_io(ji,jj) = 0.5 * ( vi_ice(ji+1,jj+1) + vi_ice(ji ,jj+1) ) - ssv_m(ji,jj)
[3]	157	END DO
	158	END DO
[717]	159	! frictional velocity at T-point
[3]	160	DO jj = 2, jpjm1
[717]	161	DO ji = fs_2, fs_jpim1 ! vector opt.
	162	ust2s(ji,jj) = 0.5 * cw &
	163	& * ( zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj) &
	164	& + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) * tms(ji,jj)
[3]	165	END DO
	166	END DO
[717]	167	!
	168	ELSE ! no ice dynamics : transmit directly the atmospheric stress to the ocean
	169	!
	170	zcoef = SQRT( 0.5 ) / rau0
	171	DO jj = 2, jpjm1
	172	DO ji = fs_2, fs_jpim1 ! vector opt.
	173	ust2s(ji,jj) = zcoef * tms(ji,jj) * SQRT( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) &
	174	& + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) )
[3]	175	END DO
	176	END DO
[717]	177	!
[3]	178	ENDIF
[879]	179	!
[3]	180	CALL lbc_lnk( ust2s, 'T', 1. ) ! T-point
[879]	181	!
[717]	182	IF(ln_ctl) CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn : ust2s :')
[3]	183
[881]	184	END SUBROUTINE lim_dyn_2
[3]	185
[76]	186
[881]	187	SUBROUTINE lim_dyn_init_2
[3]	188	!!-------------------------------------------------------------------
[881]	189	!! * ROUTINE lim_dyn_init_2 *
[3]	190	!!
[717]	191	!! ** Purpose : Physical constants and parameters linked to the ice
	192	!! dynamics
[3]	193	!!
[717]	194	!! ** Method : Read the namicedyn namelist and check the ice-dynamic
	195	!! parameter values
[3]	196	!!
	197	!! ** input : Namelist namicedyn
	198	!!-------------------------------------------------------------------
[12]	199	NAMELIST/namicedyn/ epsd, alpha, &
[3]	200	& dm, nbiter, nbitdr, om, resl, cw, angvg, pstar, &
	201	& c_rhg, etamn, creepl, ecc, ahi0
	202	!!-------------------------------------------------------------------
	203
[717]	204	REWIND ( numnam_ice ) ! Read Namelist namicedyn
	205	READ ( numnam_ice , namicedyn )
[3]	206
[717]	207	IF(lwp) THEN ! Control print
[3]	208	WRITE(numout,*)
[881]	209	WRITE(numout,*) 'lim_dyn_init_2: ice parameters for ice dynamics '
	210	WRITE(numout,*) '~~~~~~~~~~~~~~'
[76]	211	WRITE(numout,*) ' tolerance parameter epsd = ', epsd
	212	WRITE(numout,*) ' coefficient for semi-implicit coriolis alpha = ', alpha
	213	WRITE(numout,*) ' diffusion constant for dynamics dm = ', dm
	214	WRITE(numout,*) ' number of sub-time steps for relaxation nbiter = ', nbiter
	215	WRITE(numout,*) ' maximum number of iterations for relaxation nbitdr = ', nbitdr
	216	WRITE(numout,*) ' relaxation constant om = ', om
	217	WRITE(numout,*) ' maximum value for the residual of relaxation resl = ', resl
	218	WRITE(numout,*) ' drag coefficient for oceanic stress cw = ', cw
[717]	219	WRITE(numout,*) ' turning angle for oceanic stress angvg = ', angvg, ' degrees'
[76]	220	WRITE(numout,*) ' first bulk-rheology parameter pstar = ', pstar
	221	WRITE(numout,*) ' second bulk-rhelogy parameter c_rhg = ', c_rhg
	222	WRITE(numout,*) ' minimun value for viscosity etamn = ', etamn
	223	WRITE(numout,*) ' creep limit creepl = ', creepl
	224	WRITE(numout,*) ' eccentricity of the elliptical yield curve ecc = ', ecc
	225	WRITE(numout,*) ' horizontal diffusivity coeff. for sea-ice ahi0 = ', ahi0
[3]	226	ENDIF
	227
	228	usecc2 = 1.0 / ( ecc * ecc )
	229	rhoco = rau0 * cw
[717]	230	angvg = angvg * rad ! convert angvg from degree to radian
[3]	231	sangvg = SIN( angvg )
	232	cangvg = COS( angvg )
	233	pstarh = pstar / 2.0
[717]	234	!
	235	ahiu(:,:) = ahi0 * umask(:,:,1) ! Ice eddy Diffusivity coefficients.
[3]	236	ahiv(:,:) = ahi0 * vmask(:,:,1)
[717]	237	!
[881]	238	END SUBROUTINE lim_dyn_init_2
[3]	239
	240	#else
	241	!!----------------------------------------------------------------------
[881]	242	!! Default option Empty module NO LIM 2.0 sea-ice model
[3]	243	!!----------------------------------------------------------------------
	244	CONTAINS
[881]	245	SUBROUTINE lim_dyn_2 ! Empty routine
	246	END SUBROUTINE lim_dyn_2
[3]	247	#endif
	248
	249	!!======================================================================
[881]	250	END MODULE limdyn_2

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

Download in other formats: