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.

ldfdyn_c3d.h90 in branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

Context Navigation

source: branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90 @ 7226

Last change on this file since 7226 was 5312, checked in by timgraham, 9 years ago
Reset svn:keywords Id property
Property svn:keywords set to `Id`
File size: 18.5 KB

Rev	Line
[3]	1	!!----------------------------------------------------------------------
[690]	2	!! * ldfdyn_c3d.h90 *
[3]	3	!!----------------------------------------------------------------------
	4
	5	!!----------------------------------------------------------------------
[2528]	6	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[5235]	7	!! $Id$
[2528]	8	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[247]	9	!!----------------------------------------------------------------------
	10
[690]	11	!!----------------------------------------------------------------------
	12	!! 'key_dynldf_c3d' 3D lateral eddy viscosity coefficients
	13	!!----------------------------------------------------------------------
	14
	15	SUBROUTINE ldf_dyn_c3d( ld_print )
[3]	16	!!----------------------------------------------------------------------
[690]	17	!! * ROUTINE ldf_dyn_c3d *
	18	!!
[3]	19	!! ** Purpose : initializations of the horizontal ocean physics
	20	!!
[690]	21	!! ** Method : 3D eddy viscosity coef. ( longitude, latitude, depth )
	22	!! laplacian operator : ahm1, ahm2 defined at T- and F-points
	23	!! ahm2, ahm4 never used
	24	!! bilaplacian operator : ahm1, ahm2 never used
	25	!! : ahm3, ahm4 defined at U- and V-points
	26	!! ??? explanation of the default is missing
[3]	27	!!----------------------------------------------------------------------
[2715]	28	USE ldftra_oce, ONLY : aht0
[2528]	29	!!
[2715]	30	LOGICAL, INTENT (in) :: ld_print ! If true, output arrays on numout
[2528]	31	!!
[2715]	32	INTEGER :: ji, jj, jk ! dummy loop indices
	33	REAL(wp) :: zr = 0.2 ! maximum of the reduction factor at the bottom ocean ( 0 < zr < 1 )
	34	REAL(wp) :: zh = 500. ! depth of at which start the reduction ( > dept(1) )
	35	REAL(wp) :: zd_max ! maximum grid spacing over the global domain
	36	REAL(wp) :: za00, zc, zd, zetmax, zefmax, zeumax, zevmax ! local scalars
[3294]	37	REAL(wp), POINTER, DIMENSION(:) :: zcoef
[3]	38	!!----------------------------------------------------------------------
[3294]	39	!
	40	CALL wrk_alloc( jpk, zcoef )
	41	!
[3]	42	IF(lwp) WRITE(numout,*)
[690]	43	IF(lwp) WRITE(numout,*) 'ldf_dyn_c3d : 3D lateral eddy viscosity coefficient'
[3]	44	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	45
[690]	46
	47	! Set ahm1 and ahm2 ( T- and F- points) (used for laplacian operators
	48	! ================= whatever its orientation is)
[3]	49	IF( ln_dynldf_lap ) THEN
	50	! define ahm1 and ahm2 at the right grid point position
	51	! (USER: modify ahm1 and ahm2 following your desiderata)
	52
[901]	53	zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
	54	IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain
[32]	55
[3]	56	IF(lwp) WRITE(numout,*) ' laplacian operator: ahm proportional to e1'
[901]	57	IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
[3]	58
[901]	59	za00 = ahm0 / zd_max
[690]	60
[3]	61	IF( ln_dynldf_iso ) THEN
	62	IF(lwp) WRITE(numout,*) ' Caution, as implemented now, the isopycnal part of momentum'
	63	IF(lwp) WRITE(numout,*) ' mixing use aht0 as eddy viscosity coefficient. Thus, it is'
	64	IF(lwp) WRITE(numout,*) ' uniform and you must be sure that your ahm is greater than'
	65	IF(lwp) WRITE(numout,*) ' aht0 everywhere in the model domain.'
	66	ENDIF
	67
[690]	68	CALL ldf_zpf( .TRUE. , 1000., 500., 0.25, fsdept(:,:,:), ahm1 ) ! vertical profile
	69	CALL ldf_zpf( .TRUE. , 1000., 500., 0.25, fsdept(:,:,:), ahm2 ) ! vertical profile
	70	DO jk = 1,jpk
[901]	71	DO jj = 1, jpj
	72	DO ji = 1, jpi
	73	zetmax = MAX( e1t(ji,jj), e2t(ji,jj) )
	74	zefmax = MAX( e1f(ji,jj), e2f(ji,jj) )
	75	ahm1(ji,jj,jk) = za00 * zetmax * ahm1(ji,jj,jk)
	76	ahm2(ji,jj,jk) = za00 * zefmax * ahm2(ji,jj,jk)
	77	END DO
	78	END DO
[690]	79	END DO
	80
	81
[2528]	82	! Special case for ORCA R1, R2 and R4 configurations (overwrite the value of ahm1 ahm2)
[3]	83	! ==============================================
[2528]	84	IF( cp_cfg == "orca" .AND. ( jp_cfg == 1 .OR. jp_cfg == 2 .OR. jp_cfg == 4 ) ) THEN
[690]	85	IF(lwp) WRITE(numout,*)
[2528]	86	IF(lwp) WRITE(numout,*) ' ORCA R1, R2 or R4: overwrite the previous definition of ahm'
	87	IF(lwp) WRITE(numout,*) ' ================='
[690]	88	CALL ldf_dyn_c3d_orca( ld_print )
	89	ENDIF
[689]	90
[3]	91	ENDIF
[690]	92
	93	! Control print
	94	IF(lwp .AND. ld_print ) THEN
	95	WRITE(numout,*)
	96	WRITE(numout,*) ' 3D ahm1 array (k=1)'
	97	CALL prihre( ahm1(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
	98	WRITE(numout,*)
	99	WRITE(numout,*) ' 3D ahm2 array (k=1)'
	100	CALL prihre( ahm2(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
	101	ENDIF
[3]	102
[690]	103
	104	! ahm3 and ahm4 at U- and V-points (used for bilaplacian operator
	105	! ================================ whatever its orientation is)
	106	! (USER: modify ahm3 and ahm4 following your desiderata)
	107	! Here: ahm is proportional to the cube of the maximum of the gridspacing
	108	! in the to horizontal direction
	109
[3]	110	IF( ln_dynldf_bilap ) THEN
	111
[901]	112	zd_max = MAX( MAXVAL( e1u(:,:) ), MAXVAL( e2u(:,:) ) )
	113	IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain
[32]	114
[3]	115	IF(lwp) WRITE(numout,) ' bi-laplacian operator: ahm proportional to e1*3 '
[901]	116	IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
[3]	117
[2528]	118	za00 = ahm0_blp / ( zd_max * zd_max * zd_max )
[901]	119	DO jj = 1, jpj
	120	DO ji = 1, jpi
	121	zeumax = MAX( e1u(ji,jj), e2u(ji,jj) )
	122	zevmax = MAX( e1v(ji,jj), e2v(ji,jj) )
	123	ahm3(ji,jj,1) = za00 * zeumax * zeumax * zeumax
	124	ahm4(ji,jj,1) = za00 * zevmax * zevmax * zevmax
	125	END DO
	126	END DO
[3]	127
[690]	128	zh = MAX( zh, fsdept(1,1,1) ) ! at least the first reach ahm0
	129	IF( ln_zco ) THEN ! z-coordinate, same profile everywhere
	130	IF(lwp) WRITE(numout,'(36x," ahm ", 7x)')
	131	DO jk = 1, jpk
	132	IF( fsdept(1,1,jk) <= zh ) THEN
	133	zcoef(jk) = 1.e0
	134	ELSE
	135	zcoef(jk) = 1.e0 + ( zr - 1.e0 ) &
	136	& * ( 1. - EXP( ( fsdept(1,1,jk ) - zh ) / zh ) ) &
	137	& / ( 1. - EXP( ( fsdept(1,1,jpkm1) - zh ) / zh ) )
	138	ENDIF
	139	ahm3(:,:,jk) = ahm3(:,:,1) * zcoef(jk)
	140	ahm4(:,:,jk) = ahm4(:,:,1) * zcoef(jk)
	141	IF(lwp) WRITE(numout,'(34x,E7.2,8x,i3)') zcoef(jk) * ahm0, jk
	142	END DO
	143	ELSE ! partial steps or s-ccordinate
	144	zc = MAXVAL( fsdept(:,:,jpkm1) )
	145	IF( lk_mpp ) CALL mpp_max( zc ) ! max over the global domain
	146
	147	zc = 1. / ( 1. - EXP( ( zc - zh ) / zh ) )
	148	DO jk = 2, jpkm1
	149	DO jj = 1, jpj
	150	DO ji = 1, jpi
	151	IF( fsdept(ji,jj,jk) <= zh ) THEN
	152	ahm3(ji,jj,jk) = ahm3(ji,jj,1)
	153	ahm4(ji,jj,jk) = ahm4(ji,jj,1)
	154	ELSE
	155	zd = 1.e0 + ( zr - 1.e0 ) * ( 1. - EXP( ( fsdept(ji,jj,jk) - zh ) / zh ) ) * zc
	156	ahm3(ji,jj,jk) = ahm3(ji,jj,1) * zd
	157	ahm4(ji,jj,jk) = ahm4(ji,jj,1) * zd
	158	ENDIF
	159	END DO
	160	END DO
	161	END DO
	162	ahm3(:,:,jpk) = ahm3(:,:,jpkm1)
	163	ahm4(:,:,jpk) = ahm4(:,:,jpkm1)
	164	IF(lwp) WRITE(numout,'(36x," ahm ", 7x)')
	165	DO jk = 1, jpk
	166	IF(lwp) WRITE(numout,'(30x,E10.2,8x,i3)') ahm3(1,1,jk), jk
	167	END DO
	168	ENDIF
	169
[3]	170	! Control print
	171	IF( lwp .AND. ld_print ) THEN
	172	WRITE(numout,*)
[690]	173	WRITE(numout,*) 'inildf: ahm3 array at level 1'
	174	CALL prihre(ahm3(:,:,1 ),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[3]	175	WRITE(numout,*)
[690]	176	WRITE(numout,*) 'inildf: ahm4 array at level 1'
	177	CALL prihre(ahm4(:,:,jpk),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[3]	178	ENDIF
	179	ENDIF
[2715]	180	!
[3294]	181	CALL wrk_dealloc( jpk, zcoef )
[2715]	182	!
[690]	183	END SUBROUTINE ldf_dyn_c3d
[3]	184
	185
[690]	186	SUBROUTINE ldf_dyn_c3d_orca( ld_print )
[3]	187	!!----------------------------------------------------------------------
[690]	188	!! * ROUTINE ldf_dyn_c3d *
	189	!!
[2528]	190	!! ** Purpose : ORCA R1, R2 and R4 only
[3]	191	!!
[690]	192	!! ** Method : blah blah blah ....
[3]	193	!!----------------------------------------------------------------------
[2715]	194	USE ldftra_oce, ONLY: aht0
[2528]	195	!!
[2715]	196	LOGICAL, INTENT(in) :: ld_print ! If true, output arrays on numout
[2528]	197	!!
[690]	198	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[2715]	199	INTEGER :: ii0, ii1, ij0, ij1 ! local integers
	200	INTEGER :: inum, iim, ijm !
[32]	201	INTEGER :: ifreq, il1, il2, ij, ii
[2715]	202	REAL(wp) :: zahmeq, zcoff, zcoft, zmsk ! local scalars
	203	REAL(wp) :: zemax , zemin, zeref, zahmm
[3]	204	CHARACTER (len=15) :: clexp
[3294]	205	INTEGER , POINTER, DIMENSION(:,:) :: icof
	206	INTEGER , POINTER, DIMENSION(:,:) :: idata
	207	REAL(wp), POINTER, DIMENSION(: ) :: zcoef
	208	REAL(wp), POINTER, DIMENSION(:,:) :: zahm0
[3]	209	!!----------------------------------------------------------------------
[3294]	210	!
	211	CALL wrk_alloc( jpi , jpj , icof )
	212	CALL wrk_alloc( jpidta, jpjdta, idata )
	213	CALL wrk_alloc( jpk , zcoef )
	214	CALL wrk_alloc( jpi , jpj , zahm0 )
	215	!
[3]	216	IF(lwp) WRITE(numout,*)
[690]	217	IF(lwp) WRITE(numout,*) 'ldfdyn_c3d_orca : 3D eddy viscosity coefficient'
	218	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
[2715]	219	IF(lwp) WRITE(numout,*) ' orca R1, R2 or R4 configuration: reduced in the surface Eq. strip '
[3]	220
	221	! Read 2d integer array to specify western boundary increase in the
	222	! ===================== equatorial strip (20N-20S) defined at t-points
	223
[1581]	224	CALL ctl_opn( inum, 'ahmcoef', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
[32]	225	READ(inum,9101) clexp, iim, ijm
	226	READ(inum,'(/)')
[3]	227	ifreq = 40
	228	il1 = 1
	229	DO jn = 1, jpidta/ifreq+1
[32]	230	READ(inum,'(/)')
[3]	231	il2 = MIN( jpidta, il1+ifreq-1 )
[32]	232	READ(inum,9201) ( ii, ji = il1, il2, 5 )
	233	READ(inum,'(/)')
[3]	234	DO jj = jpjdta, 1, -1
[32]	235	READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 )
[3]	236	END DO
	237	il1 = il1 + ifreq
	238	END DO
	239
	240	DO jj = 1, nlcj
	241	DO ji = 1, nlci
	242	icof(ji,jj) = idata( mig(ji), mjg(jj) )
	243	END DO
	244	END DO
	245	DO jj = nlcj+1, jpj
	246	DO ji = 1, nlci
	247	icof(ji,jj) = icof(ji,nlcj)
	248	END DO
	249	END DO
	250	DO jj = 1, jpj
	251	DO ji = nlci+1, jpi
	252	icof(ji,jj) = icof(nlci,jj)
	253	END DO
	254	END DO
[690]	255
	256	9101 FORMAT(1x,a15,2i8)
	257	9201 FORMAT(3x,13(i3,12x))
	258	9202 FORMAT(i3,41i3)
	259
	260	! Set ahm1 and ahm2
[3]	261	! =================
[690]	262
[3]	263	! define ahm1 and ahm2 at the right grid point position
	264	! (USER: modify ahm1 and ahm2 following your desiderata)
[690]	265	! biharmonic : ahm1 (ahm2) defined at u- (v-) point
	266	! harmonic : ahm1 (ahm2) defined at t- (f-) point
[3]	267
[690]	268	! first level : as for 2D coefficients
[3]	269
	270	! Decrease ahm to zahmeq m2/s in the tropics
	271	! (from 90 to 20 degre: ahm = constant
	272	! from 20 to 2.5 degre: ahm = decrease in (1-cos)/2
	273	! from 2.5 to 0 degre: ahm = constant
	274	! symmetric in the south hemisphere)
[690]	275
	276	IF( jp_cfg == 4 ) THEN
	277	zahmeq = 5.0 * aht0
	278	zahmm = min( 160000.0, ahm0)
	279	zemax = MAXVAL ( e1t(:,:) * e2t(:,:), tmask(:,:,1) .GE. 0.5 )
	280	zemin = MINVAL ( e1t(:,:) * e2t(:,:), tmask(:,:,1) .GE. 0.5 )
	281	zeref = MAXVAL ( e1t(:,:) * e2t(:,:), &
	282	& tmask(:,:,1) .GE. 0.5 .AND. ABS(gphit(:,:)) .GT. 50. )
	283
	284	DO jj = 1, jpj
	285	DO ji = 1, jpi
	286	zmsk = e1t(ji,jj) * e2t(ji,jj)
	287	IF( abs(gphit(ji,jj)) .LE. 15 ) THEN
	288	zahm0(ji,jj) = ahm0
	289	ELSE
	290	IF( zmsk .GE. zeref ) THEN
	291	zahm0(ji,jj) = ahm0
	292	ELSE
	293	zahm0(ji,jj) = zahmm + (ahm0-zahmm)*(1.0 - &
	294	& cos((rpi0.5(zmsk-zemin)/(zeref-zemin))))
	295	ENDIF
	296	ENDIF
	297	END DO
	298	END DO
	299	ENDIF
[689]	300
[690]	301	IF( jp_cfg == 2 ) THEN
	302	zahmeq = aht0
	303	zahmm = ahm0
	304	zahm0(:,:) = ahm0
	305	ENDIF
	306
[2528]	307	IF( jp_cfg == 1 ) THEN
	308	zahmeq = aht0 ! reduced to aht0 on equator; set to ahm0 if no tropical reduction is required
	309	zahmm = ahm0
	310	zahm0(:,:) = ahm0
	311	ENDIF
	312
[3]	313	DO jj = 1, jpj
	314	DO ji = 1, jpi
[690]	315	IF( ABS(gphif(ji,jj)) >= 20.) THEN
	316	ahm2(ji,jj,1) = zahm0(ji,jj)
	317	ELSEIF( ABS(gphif(ji,jj)) <= 2.5) THEN
	318	ahm2(ji,jj,1) = zahmeq
[3]	319	ELSE
[690]	320	ahm2(ji,jj,1) = zahmeq + (zahm0(ji,jj)-zahmeq)/2. &
	321	& (1.-COS( rad(ABS(gphif(ji,jj))-2.5)*180./17.5 ) )
[3]	322	ENDIF
[690]	323	IF( ABS(gphit(ji,jj)) >= 20.) THEN
	324	ahm1(ji,jj,1) = zahm0(ji,jj)
	325	ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
	326	ahm1(ji,jj,1) = zahmeq
[3]	327	ELSE
[690]	328	ahm1(ji,jj,1) = zahmeq + (zahm0(ji,jj)-zahmeq)/2. &
	329	& (1.-COS( rad(ABS(gphit(ji,jj))-2.5)*180./17.5 ) )
[3]	330	ENDIF
	331	END DO
	332	END DO
[690]	333
[3]	334	! increase along western boundaries of equatorial strip
	335	! t-point
	336	DO jj = 1, jpjm1
	337	DO ji = 1, jpim1
[690]	338	zcoft = float( icof(ji,jj) ) / 100.
	339	ahm1(ji,jj,1) = zcoft * zahm0(ji,jj) + (1.-zcoft) * ahm1(ji,jj,1)
[3]	340	END DO
	341	END DO
	342	! f-point
	343	icof(:,:) = icof(:,:) * tmask(:,:,1)
	344	DO jj = 1, jpjm1
[1694]	345	DO ji = 1, jpim1 ! NO vector opt.
[3]	346	zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1)
	347	IF( zmsk == 0. ) THEN
	348	zcoff = 1.
	349	ELSE
	350	zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) ) &
	351	/ (zmsk * 100.)
	352	ENDIF
[690]	353	ahm2(ji,jj,1) = zcoff * zahm0(ji,jj) + (1.-zcoff) * ahm2(ji,jj,1)
[3]	354	END DO
	355	END DO
[690]	356
	357	! other level: re-increase the coef in the deep ocean
[2528]	358	!==================================================================
	359	! Prior to v3.3, zcoeff was hardwired according to k-index jk.
	360	!
	361	! From v3.3 onwards this has been generalised to a function of
	362	! depth so that it can be used with any number of levels.
	363	!
	364	! The function has been chosen to match the original values (shown
	365	! in the following comments) when using the standard 31 ORCA levels.
	366	! DO jk = 1, 21
	367	! zcoef(jk) = 1._wp
	368	! END DO
	369	! zcoef(22) = 2._wp
	370	! zcoef(23) = 3._wp
	371	! zcoef(24) = 5._wp
	372	! zcoef(25) = 7._wp
	373	! zcoef(26) = 9._wp
	374	! DO jk = 27, jpk
	375	! zcoef(jk) = 10._wp
	376	! END DO
	377	!==================================================================
	378
	379	IF(lwp) THEN
	380	WRITE(numout,*)
	381	WRITE(numout,*) ' 1D zcoef array '
	382	WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
	383	WRITE(numout,*)
	384	WRITE(numout,*) ' jk zcoef '
	385	ENDIF
	386
	387	DO jk=1, jpk
[4292]	388	zcoef(jk) = 1.0_wp + NINT(9.0_wp*(gdept_1d(jk)-800.0_wp)/(3000.0_wp-800.0_wp))
[2528]	389	zcoef(jk) = MIN(10.0_wp, MAX(1.0_wp, zcoef(jk)))
	390	IF(lwp) WRITE(numout,'(4x,i3,6x,f7.3)') jk,zcoef(jk)
[690]	391	END DO
	392
	393	DO jk = 2, jpk
	394	ahm1(:,:,jk) = MIN( zahm0(:,:), zcoef(jk) * ahm1(:,:,1) )
	395	ahm2(:,:,jk) = MIN( zahm0(:,:), zcoef(jk) * ahm2(:,:,1) )
	396	END DO
	397
[2528]	398	IF( jp_cfg == 4 ) THEN ! Limit AHM in Gibraltar strait
[690]	399	ij0 = 50 ; ij1 = 53
	400	ii0 = 69 ; ii1 = 71
	401	DO jk = 1, jpk
[2528]	402	ahm1(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) = MIN( zahmm, ahm1(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) )
	403	ahm2(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) = MIN( zahmm, ahm2(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),jk) )
[690]	404	END DO
	405	ENDIF
[2528]	406	CALL lbc_lnk( ahm1, 'T', 1. ) ! Lateral boundary conditions (unchanged sign)
	407	CALL lbc_lnk( ahm2, 'F', 1. )
[3]	408
[690]	409
[2528]	410	IF(lwp) THEN ! Control print
[690]	411	WRITE(numout,*)
	412	WRITE(numout,*) ' 3D ahm1 array (k=1)'
	413	CALL prihre( ahm1(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
	414	WRITE(numout,*)
	415	WRITE(numout,*) ' 3D ahm2 array (k=1)'
	416	CALL prihre( ahm2(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
	417	WRITE(numout,*)
	418	WRITE(numout,*) ' 3D ahm2 array (k=jpk)'
	419	CALL prihre( ahm2(:,:,jpk), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout )
	420	ENDIF
	421
	422
	423	! Set ahm3 and ahm4
	424	! =================
	425
	426	! define ahm3 and ahm4 at the right grid point position
	427	! initialization to a constant value
	428	! (USER: modify ahm3 and ahm4 following your desiderata)
	429	! harmonic isopycnal or geopotential:
	430	! ahm3 (ahm4) defined at u- (v-) point
	431	DO jk = 1, jpk
	432	DO jj = 2, jpj
	433	DO ji = 2, jpi
	434	ahm3(ji,jj,jk) = 0.5 * ( ahm2(ji,jj,jk) + ahm2(ji ,jj-1,jk) )
	435	ahm4(ji,jj,jk) = 0.5 * ( ahm2(ji,jj,jk) + ahm2(ji-1,jj ,jk) )
	436	END DO
	437	END DO
	438	END DO
	439	ahm3 ( :, 1, :) = ahm3 ( :, 2, :)
	440	ahm4 ( :, 1, :) = ahm4 ( :, 2, :)
	441
[2528]	442	CALL lbc_lnk( ahm3, 'U', 1. ) ! Lateral boundary conditions (unchanged sign)
	443	CALL lbc_lnk( ahm4, 'V', 1. )
[690]	444
	445	! Control print
	446
[3]	447	IF( lwp .AND. ld_print ) THEN
	448	WRITE(numout,*)
[690]	449	WRITE(numout,*) ' ahm3 array level 1'
	450	CALL prihre(ahm3(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[3]	451	WRITE(numout,*)
[690]	452	WRITE(numout,*) ' ahm4 array level 1'
	453	CALL prihre(ahm4(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[3]	454	ENDIF
[2715]	455	!
[3294]	456	CALL wrk_dealloc( jpi , jpj , icof )
	457	CALL wrk_dealloc( jpidta, jpjdta, idata )
	458	CALL wrk_dealloc( jpk , zcoef )
	459	CALL wrk_dealloc( jpi , jpj , zahm0 )
	460	!
[690]	461	END SUBROUTINE ldf_dyn_c3d_orca

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

Download in other formats:

Original Format