Context Navigation

lib_fortran.F90 @ 15725

Last change on this file since 15725 was 15376, checked in by clem, 3 years ago
add functions glob_min_vec and glob_max_vec. They would need to be merged into the generic module lib_fortran_generic.h90 at some point, along with glob_sum_vec
Property svn:keywords set to `Id`
File size: 32.6 KB

Rev	Line
[2003]	1	MODULE lib_fortran
	2	!!======================================================================
	3	!! * MODULE lib_fortran *
	4	!! Fortran utilities: includes some low levels fortran functionality
	5	!!======================================================================
[2307]	6	!! History : 3.2 ! 2010-05 (M. Dunphy, R. Benshila) Original code
[4161]	7	!! 3.4 ! 2013-06 (C. Rousset) add glob_min, glob_max
	8	!! + 3d dim. of input is fexible (jpk, jpl...)
[7646]	9	!! 4.0 ! 2016-06 (T. Lovato) double precision global sum by default
[2003]	10	!!----------------------------------------------------------------------
[2307]	11
[2003]	12	!!----------------------------------------------------------------------
[3764]	13	!! glob_sum : generic interface for global masked summation over
[2307]	14	!! the interior domain for 1 or 2 2D or 3D arrays
[3764]	15	!! it works only for T points
[2307]	16	!! SIGN : generic interface for SIGN to overwrite f95 behaviour
	17	!! of intrinsinc sign function
	18	!!----------------------------------------------------------------------
[3632]	19	USE par_oce ! Ocean parameter
	20	USE dom_oce ! ocean domain
	21	USE in_out_manager ! I/O manager
	22	USE lib_mpp ! distributed memory computing
[10425]	23	USE lbclnk ! ocean lateral boundary conditions
[2003]	24
	25	IMPLICIT NONE
	26	PRIVATE
	27
[15145]	28	PUBLIC glob_sum ! used in many places (masked with tmask_i = ssmask * tmask_h)
	29	PUBLIC glob_sum_full ! used in many places (masked with tmask_h, excluding all duplicated points halos+periodicity)
[10425]	30	PUBLIC local_sum ! used in trcrad, local operation before glob_sum_delay
	31	PUBLIC sum3x3 ! used in trcrad, do a sum over 3x3 boxes
[6140]	32	PUBLIC DDPDD ! also used in closea module
[4161]	33	PUBLIC glob_min, glob_max
[15048]	34	PUBLIC glob_sum_vec
	35	PUBLIC glob_sum_full_vec
[15376]	36	PUBLIC glob_min_vec, glob_max_vec
[2341]	37	#if defined key_nosignedzero
[2003]	38	PUBLIC SIGN
	39	#endif
	40
	41	INTERFACE glob_sum
[10425]	42	MODULE PROCEDURE glob_sum_1d, glob_sum_2d, glob_sum_3d
[2003]	43	END INTERFACE
[6140]	44	INTERFACE glob_sum_full
	45	MODULE PROCEDURE glob_sum_full_2d, glob_sum_full_3d
	46	END INTERFACE
[10425]	47	INTERFACE local_sum
	48	MODULE PROCEDURE local_sum_2d, local_sum_3d
	49	END INTERFACE
	50	INTERFACE sum3x3
	51	MODULE PROCEDURE sum3x3_2d, sum3x3_3d
	52	END INTERFACE
[4161]	53	INTERFACE glob_min
[10425]	54	MODULE PROCEDURE glob_min_2d, glob_min_3d
[4161]	55	END INTERFACE
	56	INTERFACE glob_max
[10425]	57	MODULE PROCEDURE glob_max_2d, glob_max_3d
[4161]	58	END INTERFACE
[15048]	59	INTERFACE glob_sum_vec
	60	MODULE PROCEDURE glob_sum_vec_3d, glob_sum_vec_4d
	61	END INTERFACE
	62	INTERFACE glob_sum_full_vec
	63	MODULE PROCEDURE glob_sum_full_vec_3d, glob_sum_full_vec_4d
	64	END INTERFACE
[15376]	65	INTERFACE glob_min_vec
	66	MODULE PROCEDURE glob_min_vec_3d, glob_min_vec_4d
	67	END INTERFACE
	68	INTERFACE glob_max_vec
	69	MODULE PROCEDURE glob_max_vec_3d, glob_max_vec_4d
	70	END INTERFACE
[2003]	71
[3764]	72	#if defined key_nosignedzero
[2003]	73	INTERFACE SIGN
[2307]	74	MODULE PROCEDURE SIGN_SCALAR, SIGN_ARRAY_1D, SIGN_ARRAY_2D, SIGN_ARRAY_3D, &
[3764]	75	& SIGN_ARRAY_1D_A, SIGN_ARRAY_2D_A, SIGN_ARRAY_3D_A, &
	76	& SIGN_ARRAY_1D_B, SIGN_ARRAY_2D_B, SIGN_ARRAY_3D_B
[2003]	77	END INTERFACE
	78	#endif
	79
[12377]	80	!! * Substitutions
	81	# include "do_loop_substitute.h90"
[2307]	82	!!----------------------------------------------------------------------
[9598]	83	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[3764]	84	!! $Id$
[10068]	85	!! Software governed by the CeCILL license (see ./LICENSE)
[2307]	86	!!----------------------------------------------------------------------
[3764]	87	CONTAINS
[2003]	88
[10425]	89	# define GLOBSUM_CODE
[3764]	90
[10425]	91	# define DIM_1d
	92	# define FUNCTION_GLOBSUM glob_sum_1d
	93	# include "lib_fortran_generic.h90"
	94	# undef FUNCTION_GLOBSUM
	95	# undef DIM_1d
[2307]	96
[10425]	97	# define DIM_2d
	98	# define OPERATION_GLOBSUM
	99	# define FUNCTION_GLOBSUM glob_sum_2d
	100	# include "lib_fortran_generic.h90"
	101	# undef FUNCTION_GLOBSUM
	102	# undef OPERATION_GLOBSUM
	103	# define OPERATION_FULL_GLOBSUM
	104	# define FUNCTION_GLOBSUM glob_sum_full_2d
	105	# include "lib_fortran_generic.h90"
	106	# undef FUNCTION_GLOBSUM
	107	# undef OPERATION_FULL_GLOBSUM
	108	# undef DIM_2d
[2307]	109
[10425]	110	# define DIM_3d
	111	# define OPERATION_GLOBSUM
	112	# define FUNCTION_GLOBSUM glob_sum_3d
	113	# include "lib_fortran_generic.h90"
	114	# undef FUNCTION_GLOBSUM
	115	# undef OPERATION_GLOBSUM
	116	# define OPERATION_FULL_GLOBSUM
	117	# define FUNCTION_GLOBSUM glob_sum_full_3d
	118	# include "lib_fortran_generic.h90"
	119	# undef FUNCTION_GLOBSUM
	120	# undef OPERATION_FULL_GLOBSUM
	121	# undef DIM_3d
[2307]	122
[10425]	123	# undef GLOBSUM_CODE
[2307]	124
[10425]	125
	126	# define GLOBMINMAX_CODE
	127
	128	# define DIM_2d
	129	# define OPERATION_GLOBMIN
	130	# define FUNCTION_GLOBMINMAX glob_min_2d
	131	# include "lib_fortran_generic.h90"
	132	# undef FUNCTION_GLOBMINMAX
	133	# undef OPERATION_GLOBMIN
	134	# define OPERATION_GLOBMAX
	135	# define FUNCTION_GLOBMINMAX glob_max_2d
	136	# include "lib_fortran_generic.h90"
	137	# undef FUNCTION_GLOBMINMAX
	138	# undef OPERATION_GLOBMAX
	139	# undef DIM_2d
	140
	141	# define DIM_3d
	142	# define OPERATION_GLOBMIN
	143	# define FUNCTION_GLOBMINMAX glob_min_3d
	144	# include "lib_fortran_generic.h90"
	145	# undef FUNCTION_GLOBMINMAX
	146	# undef OPERATION_GLOBMIN
	147	# define OPERATION_GLOBMAX
	148	# define FUNCTION_GLOBMINMAX glob_max_3d
	149	# include "lib_fortran_generic.h90"
	150	# undef FUNCTION_GLOBMINMAX
	151	# undef OPERATION_GLOBMAX
	152	# undef DIM_3d
	153	# undef GLOBMINMAX_CODE
	154
	155	! ! FUNCTION local_sum !
	156
	157	FUNCTION local_sum_2d( ptab )
[2307]	158	!!----------------------------------------------------------------------
[10425]	159	REAL(wp), INTENT(in ) :: ptab(:,:) ! array on which operation is applied
[13226]	160	COMPLEX(dp) :: local_sum_2d
[10425]	161	!
[2003]	162	!!-----------------------------------------------------------------------
[2307]	163	!
[13226]	164	COMPLEX(dp):: ctmp
[2307]	165	REAL(wp) :: ztmp
[10425]	166	INTEGER :: ji, jj ! dummy loop indices
	167	INTEGER :: ipi, ipj ! dimensions
[2307]	168	!!-----------------------------------------------------------------------
	169	!
[10425]	170	ipi = SIZE(ptab,1) ! 1st dimension
	171	ipj = SIZE(ptab,2) ! 2nd dimension
[4161]	172	!
[10425]	173	ctmp = CMPLX( 0.e0, 0.e0, wp ) ! warning ctmp is cumulated
	174
	175	DO jj = 1, ipj
	176	DO ji = 1, ipi
	177	ztmp = ptab(ji,jj) * tmask_i(ji,jj)
[13226]	178	CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp )
[10425]	179	END DO
[2307]	180	END DO
	181	!
[10425]	182	local_sum_2d = ctmp
	183
	184	END FUNCTION local_sum_2d
[2307]	185
[10425]	186	FUNCTION local_sum_3d( ptab )
[6140]	187	!!----------------------------------------------------------------------
[10425]	188	REAL(wp), INTENT(in ) :: ptab(:,:,:) ! array on which operation is applied
[13226]	189	COMPLEX(dp) :: local_sum_3d
[10425]	190	!
[6140]	191	!!-----------------------------------------------------------------------
	192	!
[13226]	193	COMPLEX(dp):: ctmp
[6140]	194	REAL(wp) :: ztmp
	195	INTEGER :: ji, jj, jk ! dummy loop indices
[10425]	196	INTEGER :: ipi, ipj, ipk ! dimensions
[6140]	197	!!-----------------------------------------------------------------------
	198	!
[10425]	199	ipi = SIZE(ptab,1) ! 1st dimension
	200	ipj = SIZE(ptab,2) ! 2nd dimension
	201	ipk = SIZE(ptab,3) ! 3rd dimension
[6140]	202	!
[10425]	203	ctmp = CMPLX( 0.e0, 0.e0, wp ) ! warning ctmp is cumulated
[6140]	204
[10425]	205	DO jk = 1, ipk
	206	DO jj = 1, ipj
	207	DO ji = 1, ipi
	208	ztmp = ptab(ji,jj,jk) * tmask_i(ji,jj)
[13226]	209	CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp )
[10425]	210	END DO
	211	END DO
[4161]	212	END DO
	213	!
[10425]	214	local_sum_3d = ctmp
	215
	216	END FUNCTION local_sum_3d
[4161]	217
[10425]	218	! ! FUNCTION sum3x3 !
[4161]	219
[10425]	220	SUBROUTINE sum3x3_2d( p2d )
[4161]	221	!!-----------------------------------------------------------------------
[10425]	222	!! * routine sum3x3_2d *
[4161]	223	!!
[10425]	224	!! ** Purpose : sum over 3x3 boxes
	225	!!----------------------------------------------------------------------
	226	REAL(wp), DIMENSION (:,:), INTENT(inout) :: p2d
[4161]	227	!
[10425]	228	INTEGER :: ji, ji2, jj, jj2 ! dummy loop indices
	229	!!----------------------------------------------------------------------
[4161]	230	!
[10425]	231	IF( SIZE(p2d,1) /= jpi ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_2d, the first dimension is not equal to jpi' )
	232	IF( SIZE(p2d,2) /= jpj ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_2d, the second dimension is not equal to jpj' )
[4161]	233	!
[13324]	234	! work over the whole domain (guarantees all internal cells are set when nn_hls=2)
	235	!
	236	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
[14433]	237	IF( MOD(mig(ji), 3) == MOD(nn_hls, 3) .AND. & ! 1st bottom left corner always at (Nis0-1, Njs0-1)
[13327]	238	& MOD(mjg(jj), 3) == MOD(nn_hls, 3) ) THEN ! bottom left corner of a 3x3 box
[12377]	239	ji2 = MIN(mig(ji)+2, jpiglo) - nimpp + 1 ! right position of the box
	240	jj2 = MIN(mjg(jj)+2, jpjglo) - njmpp + 1 ! upper position of the box
	241	IF( ji2 <= jpi .AND. jj2 <= jpj ) THEN ! the box is fully included in the local mpi domain
	242	p2d(ji:ji2,jj:jj2) = SUM(p2d(ji:ji2,jj:jj2))
[10425]	243	ENDIF
[12377]	244	ENDIF
	245	END_2D
[13226]	246	CALL lbc_lnk( 'lib_fortran', p2d, 'T', 1.0_wp )
[14433]	247	! no need for 2nd exchange when nn_hls > 1
	248	IF( nn_hls == 1 ) THEN
	249	IF( mpiRnei(nn_hls,jpwe) > -1 ) THEN ! 1st column was changed during the previous call to lbc_lnk
	250	IF( MOD(mig( 1), 3) == 1 ) & ! 1st box start at i=1 -> column 1 to 3 correctly computed locally
	251	p2d( 1,:) = p2d( 2,:) ! previous lbc_lnk corrupted column 1 -> put it back using column 2
	252	IF( MOD(mig( 1), 3) == 2 ) & ! 1st box start at i=3 -> column 1 and 2 correctly computed on west neighbourh
	253	p2d( 2,:) = p2d( 1,:) ! previous lbc_lnk fix column 1 -> copy it to column 2
[13324]	254	ENDIF
[14433]	255	IF( mpiRnei(nn_hls,jpea) > -1 ) THEN
[13324]	256	IF( MOD(mig(jpi-2), 3) == 1 ) p2d( jpi,:) = p2d(jpi-1,:)
	257	IF( MOD(mig(jpi-2), 3) == 0 ) p2d(jpi-1,:) = p2d( jpi,:)
	258	ENDIF
[14433]	259	IF( mpiRnei(nn_hls,jpso) > -1 ) THEN
[13324]	260	IF( MOD(mjg( 1), 3) == 1 ) p2d(:, 1) = p2d(:, 2)
	261	IF( MOD(mjg( 1), 3) == 2 ) p2d(:, 2) = p2d(:, 1)
	262	ENDIF
[14433]	263	IF( mpiRnei(nn_hls,jpno) > -1 ) THEN
[13324]	264	IF( MOD(mjg(jpj-2), 3) == 1 ) p2d(:, jpj) = p2d(:,jpj-1)
	265	IF( MOD(mjg(jpj-2), 3) == 0 ) p2d(:,jpj-1) = p2d(:, jpj)
	266	ENDIF
	267	CALL lbc_lnk( 'lib_fortran', p2d, 'T', 1.0_wp )
[10425]	268	ENDIF
[4161]	269
[10425]	270	END SUBROUTINE sum3x3_2d
	271
	272	SUBROUTINE sum3x3_3d( p3d )
[4161]	273	!!-----------------------------------------------------------------------
[10425]	274	!! * routine sum3x3_3d *
[4161]	275	!!
[10425]	276	!! ** Purpose : sum over 3x3 boxes
	277	!!----------------------------------------------------------------------
	278	REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: p3d
[4161]	279	!
[10425]	280	INTEGER :: ji, ji2, jj, jj2, jn ! dummy loop indices
	281	INTEGER :: ipn ! Third dimension size
	282	!!----------------------------------------------------------------------
[4161]	283	!
[10425]	284	IF( SIZE(p3d,1) /= jpi ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_3d, the first dimension is not equal to jpi' )
	285	IF( SIZE(p3d,2) /= jpj ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_3d, the second dimension is not equal to jpj' )
	286	ipn = SIZE(p3d,3)
[4161]	287	!
[10425]	288	DO jn = 1, ipn
[13324]	289	!
	290	! work over the whole domain (guarantees all internal cells are set when nn_hls=2)
	291	!
	292	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
[14433]	293	IF( MOD(mig(ji), 3) == MOD(nn_hls, 3) .AND. & ! 1st bottom left corner always at (Nis0-1, Njs0-1)
[13327]	294	& MOD(mjg(jj), 3) == MOD(nn_hls, 3) ) THEN ! bottom left corner of a 3x3 box
[12377]	295	ji2 = MIN(mig(ji)+2, jpiglo) - nimpp + 1 ! right position of the box
	296	jj2 = MIN(mjg(jj)+2, jpjglo) - njmpp + 1 ! upper position of the box
	297	IF( ji2 <= jpi .AND. jj2 <= jpj ) THEN ! the box is fully included in the local mpi domain
	298	p3d(ji:ji2,jj:jj2,jn) = SUM(p3d(ji:ji2,jj:jj2,jn))
[10425]	299	ENDIF
[12377]	300	ENDIF
	301	END_2D
[4161]	302	END DO
[13226]	303	CALL lbc_lnk( 'lib_fortran', p3d, 'T', 1.0_wp )
[14433]	304	! no need for 2nd exchange when nn_hls > 1
	305	IF( nn_hls == 1 ) THEN
	306	IF( mpiRnei(nn_hls,jpwe) > -1 ) THEN ! 1st column was changed during the previous call to lbc_lnk
	307	IF( MOD(mig( 1), 3) == 1 ) & ! 1st box start at i=1 -> column 1 to 3 correctly computed locally
	308	p3d( 1,:,:) = p3d( 2,:,:) ! previous lbc_lnk corrupted column 1 -> put it back using column 2
	309	IF( MOD(mig( 1), 3) == 2 ) & ! 1st box start at i=3 -> column 1 and 2 correctly computed on west neighbourh
	310	p3d( 2,:,:) = p3d( 1,:,:) ! previous lbc_lnk fix column 1 -> copy it to column 2
[13324]	311	ENDIF
[14433]	312	IF( mpiRnei(nn_hls,jpea) > -1 ) THEN
[13324]	313	IF( MOD(mig(jpi-2), 3) == 1 ) p3d( jpi,:,:) = p3d(jpi-1,:,:)
	314	IF( MOD(mig(jpi-2), 3) == 0 ) p3d(jpi-1,:,:) = p3d( jpi,:,:)
	315	ENDIF
[14433]	316	IF( mpiRnei(nn_hls,jpso) > -1 ) THEN
[13324]	317	IF( MOD(mjg( 1), 3) == 1 ) p3d(:, 1,:) = p3d(:, 2,:)
	318	IF( MOD(mjg( 1), 3) == 2 ) p3d(:, 2,:) = p3d(:, 1,:)
	319	ENDIF
[14433]	320	IF( mpiRnei(nn_hls,jpno) > -1 ) THEN
[13324]	321	IF( MOD(mjg(jpj-2), 3) == 1 ) p3d(:, jpj,:) = p3d(:,jpj-1,:)
	322	IF( MOD(mjg(jpj-2), 3) == 0 ) p3d(:,jpj-1,:) = p3d(:, jpj,:)
	323	ENDIF
	324	CALL lbc_lnk( 'lib_fortran', p3d, 'T', 1.0_wp )
[10425]	325	ENDIF
[4161]	326
[10425]	327	END SUBROUTINE sum3x3_3d
[4161]	328
	329
[15048]	330	FUNCTION glob_sum_vec_3d( cdname, ptab ) RESULT( ptmp )
	331	!!----------------------------------------------------------------------
	332	CHARACTER(len=*), INTENT(in) :: cdname ! name of the calling subroutine
	333	REAL(wp), INTENT(in) :: ptab(:,:,:) ! array on which operation is applied
	334	REAL(wp), DIMENSION(SIZE(ptab,3)) :: ptmp
	335	!
	336	COMPLEX(dp), DIMENSION(:), ALLOCATABLE :: ctmp
	337	REAL(wp) :: ztmp
	338	INTEGER :: ji , jj , jk ! dummy loop indices
	339	INTEGER :: ipi, ipj, ipk ! dimensions
[15145]	340	INTEGER :: iis, iie, ijs, ije ! loop start and end
[15048]	341	!!-----------------------------------------------------------------------
	342	!
	343	ipi = SIZE(ptab,1) ! 1st dimension
	344	ipj = SIZE(ptab,2) ! 2nd dimension
	345	ipk = SIZE(ptab,3) ! 3rd dimension
	346	!
[15145]	347	IF( ipi == jpi .AND. ipj == jpj ) THEN ! do 2D loop only over the inner domain (-> avoid to use undefined values)
	348	iis = Nis0 ; iie = Nie0
	349	ijs = Njs0 ; ije = Nje0
	350	ELSE ! I think we are never in this case...
	351	iis = 1 ; iie = jpi
	352	ijs = 1 ; ije = jpj
	353	ENDIF
	354	!
[15048]	355	ALLOCATE( ctmp(ipk) )
	356	!
	357	DO jk = 1, ipk
	358	ctmp(jk) = CMPLX( 0.e0, 0.e0, dp ) ! warning ctmp is cumulated
[15145]	359	DO jj = ijs, ije
	360	DO ji = iis, iie
[15048]	361	ztmp = ptab(ji,jj,jk) * tmask_i(ji,jj)
	362	CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jk) )
	363	END DO
	364	END DO
	365	END DO
	366	CALL mpp_sum( cdname, ctmp(:) ) ! sum over the global domain
	367	!
	368	ptmp = REAL( ctmp(:), wp )
	369	!
	370	DEALLOCATE( ctmp )
	371	!
	372	END FUNCTION glob_sum_vec_3d
	373
	374	FUNCTION glob_sum_vec_4d( cdname, ptab ) RESULT( ptmp )
	375	!!----------------------------------------------------------------------
	376	CHARACTER(len=*), INTENT(in) :: cdname ! name of the calling subroutine
	377	REAL(wp), INTENT(in) :: ptab(:,:,:,:) ! array on which operation is applied
	378	REAL(wp), DIMENSION(SIZE(ptab,4)) :: ptmp
	379	!
	380	COMPLEX(dp), DIMENSION(:), ALLOCATABLE :: ctmp
	381	REAL(wp) :: ztmp
	382	INTEGER :: ji , jj , jk , jl ! dummy loop indices
	383	INTEGER :: ipi, ipj, ipk, ipl ! dimensions
[15145]	384	INTEGER :: iis, iie, ijs, ije ! loop start and end
[15048]	385	!!-----------------------------------------------------------------------
	386	!
	387	ipi = SIZE(ptab,1) ! 1st dimension
	388	ipj = SIZE(ptab,2) ! 2nd dimension
	389	ipk = SIZE(ptab,3) ! 3rd dimension
	390	ipl = SIZE(ptab,4) ! 4th dimension
	391	!
[15145]	392	IF( ipi == jpi .AND. ipj == jpj ) THEN ! do 2D loop only over the inner domain (-> avoid to use undefined values)
	393	iis = Nis0 ; iie = Nie0
	394	ijs = Njs0 ; ije = Nje0
	395	ELSE ! I think we are never in this case...
	396	iis = 1 ; iie = jpi
	397	ijs = 1 ; ije = jpj
	398	ENDIF
	399	!
[15048]	400	ALLOCATE( ctmp(ipl) )
	401	!
	402	DO jl = 1, ipl
	403	ctmp(jl) = CMPLX( 0.e0, 0.e0, dp ) ! warning ctmp is cumulated
	404	DO jk = 1, ipk
[15145]	405	DO jj = ijs, ije
	406	DO ji = iis, iie
[15048]	407	ztmp = ptab(ji,jj,jk,jl) * tmask_i(ji,jj)
	408	CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jl) )
	409	END DO
	410	END DO
	411	END DO
	412	END DO
	413	CALL mpp_sum( cdname, ctmp(:) ) ! sum over the global domain
	414	!
	415	ptmp = REAL( ctmp(:), wp )
	416	!
	417	DEALLOCATE( ctmp )
	418	!
	419	END FUNCTION glob_sum_vec_4d
	420
	421	FUNCTION glob_sum_full_vec_3d( cdname, ptab ) RESULT( ptmp )
	422	!!----------------------------------------------------------------------
	423	CHARACTER(len=*), INTENT(in) :: cdname ! name of the calling subroutine
	424	REAL(wp), INTENT(in) :: ptab(:,:,:) ! array on which operation is applied
	425	REAL(wp), DIMENSION(SIZE(ptab,3)) :: ptmp
	426	!
	427	COMPLEX(dp), DIMENSION(:), ALLOCATABLE :: ctmp
	428	REAL(wp) :: ztmp
	429	INTEGER :: ji , jj , jk ! dummy loop indices
	430	INTEGER :: ipi, ipj, ipk ! dimensions
[15145]	431	INTEGER :: iis, iie, ijs, ije ! loop start and end
[15048]	432	!!-----------------------------------------------------------------------
	433	!
	434	ipi = SIZE(ptab,1) ! 1st dimension
	435	ipj = SIZE(ptab,2) ! 2nd dimension
	436	ipk = SIZE(ptab,3) ! 3rd dimension
	437	!
[15145]	438	IF( ipi == jpi .AND. ipj == jpj ) THEN ! do 2D loop only over the inner domain (-> avoid to use undefined values)
	439	iis = Nis0 ; iie = Nie0
	440	ijs = Njs0 ; ije = Nje0
	441	ELSE ! I think we are never in this case...
	442	iis = 1 ; iie = jpi
	443	ijs = 1 ; ije = jpj
	444	ENDIF
	445	!
[15048]	446	ALLOCATE( ctmp(ipk) )
	447	!
	448	DO jk = 1, ipk
	449	ctmp(jk) = CMPLX( 0.e0, 0.e0, dp ) ! warning ctmp is cumulated
[15145]	450	DO jj = ijs, ije
	451	DO ji = iis, iie
[15048]	452	ztmp = ptab(ji,jj,jk) * tmask_h(ji,jj)
	453	CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jk) )
	454	END DO
	455	END DO
	456	END DO
	457	CALL mpp_sum( cdname, ctmp(:) ) ! sum over the global domain
	458	!
	459	ptmp = REAL( ctmp(:), wp )
	460	!
	461	DEALLOCATE( ctmp )
	462	!
	463	END FUNCTION glob_sum_full_vec_3d
	464
	465	FUNCTION glob_sum_full_vec_4d( cdname, ptab ) RESULT( ptmp )
	466	!!----------------------------------------------------------------------
	467	CHARACTER(len=*), INTENT(in) :: cdname ! name of the calling subroutine
	468	REAL(wp), INTENT(in) :: ptab(:,:,:,:) ! array on which operation is applied
	469	REAL(wp), DIMENSION(SIZE(ptab,4)) :: ptmp
	470	!
	471	COMPLEX(dp), DIMENSION(:), ALLOCATABLE :: ctmp
	472	REAL(wp) :: ztmp
	473	INTEGER :: ji , jj , jk , jl ! dummy loop indices
	474	INTEGER :: ipi, ipj, ipk, ipl ! dimensions
[15145]	475	INTEGER :: iis, iie, ijs, ije ! loop start and end
[15048]	476	!!-----------------------------------------------------------------------
	477	!
	478	ipi = SIZE(ptab,1) ! 1st dimension
	479	ipj = SIZE(ptab,2) ! 2nd dimension
	480	ipk = SIZE(ptab,3) ! 3rd dimension
	481	ipl = SIZE(ptab,4) ! 4th dimension
	482	!
[15145]	483	IF( ipi == jpi .AND. ipj == jpj ) THEN ! do 2D loop only over the inner domain (-> avoid to use undefined values)
	484	iis = Nis0 ; iie = Nie0
	485	ijs = Njs0 ; ije = Nje0
	486	ELSE ! I think we are never in this case...
	487	iis = 1 ; iie = jpi
	488	ijs = 1 ; ije = jpj
	489	ENDIF
	490	!
[15048]	491	ALLOCATE( ctmp(ipl) )
	492	!
	493	DO jl = 1, ipl
	494	ctmp(jl) = CMPLX( 0.e0, 0.e0, dp ) ! warning ctmp is cumulated
	495	DO jk = 1, ipk
[15145]	496	DO jj = ijs, ije
	497	DO ji = iis, iie
[15048]	498	ztmp = ptab(ji,jj,jk,jl) * tmask_h(ji,jj)
	499	CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jl) )
	500	END DO
	501	END DO
	502	END DO
	503	END DO
	504	CALL mpp_sum( cdname, ctmp(:) ) ! sum over the global domain
	505	!
	506	ptmp = REAL( ctmp(:), wp )
	507	!
	508	DEALLOCATE( ctmp )
	509	!
	510	END FUNCTION glob_sum_full_vec_4d
	511
[15376]	512	FUNCTION glob_min_vec_3d( cdname, ptab ) RESULT( ptmp )
	513	!!----------------------------------------------------------------------
	514	CHARACTER(len=*), INTENT(in) :: cdname ! name of the calling subroutine
	515	REAL(wp), INTENT(in) :: ptab(:,:,:) ! array on which operation is applied
	516	REAL(wp), DIMENSION(SIZE(ptab,3)) :: ptmp
	517	!
	518	INTEGER :: jk ! dummy loop indice & dimension
	519	INTEGER :: ipk ! dimension
	520	!!-----------------------------------------------------------------------
	521	!
	522	ipk = SIZE(ptab,3)
	523	DO jk = 1, ipk
	524	ptmp(jk) = MINVAL( ptab(:,:,jk) * tmask_i(:,:) )
	525	ENDDO
	526	!
	527	CALL mpp_min( cdname, ptmp (:) )
	528	!
	529	END FUNCTION glob_min_vec_3d
	530
	531	FUNCTION glob_min_vec_4d( cdname, ptab ) RESULT( ptmp )
	532	!!----------------------------------------------------------------------
	533	CHARACTER(len=*), INTENT(in) :: cdname ! name of the calling subroutine
	534	REAL(wp), INTENT(in) :: ptab(:,:,:,:) ! array on which operation is applied
	535	REAL(wp), DIMENSION(SIZE(ptab,4)) :: ptmp
	536	!
	537	INTEGER :: jk , jl ! dummy loop indice & dimension
	538	INTEGER :: ipk, ipl ! dimension
	539	!!-----------------------------------------------------------------------
	540	!
	541	ipk = SIZE(ptab,3)
	542	ipl = SIZE(ptab,4)
	543	DO jl = 1, ipl
	544	ptmp(jl) = MINVAL( ptab(:,:,1,jl) * tmask_i(:,:) )
	545	DO jk = 2, ipk
	546	ptmp(jl) = MIN( ptmp(jl), MINVAL( ptab(:,:,jk,jl) * tmask_i(:,:) ) )
	547	ENDDO
	548	ENDDO
	549	!
	550	CALL mpp_min( cdname, ptmp (:) )
	551	!
	552	END FUNCTION glob_min_vec_4d
	553
	554	FUNCTION glob_max_vec_3d( cdname, ptab ) RESULT( ptmp )
	555	!!----------------------------------------------------------------------
	556	CHARACTER(len=*), INTENT(in) :: cdname ! name of the calling subroutine
	557	REAL(wp), INTENT(in) :: ptab(:,:,:) ! array on which operation is applied
	558	REAL(wp), DIMENSION(SIZE(ptab,3)) :: ptmp
	559	!
	560	INTEGER :: jk ! dummy loop indice & dimension
	561	INTEGER :: ipk ! dimension
	562	!!-----------------------------------------------------------------------
	563	!
	564	ipk = SIZE(ptab,3)
	565	DO jk = 1, ipk
	566	ptmp(jk) = MAXVAL( ptab(:,:,jk) * tmask_i(:,:) )
	567	ENDDO
	568	!
	569	CALL mpp_max( cdname, ptmp (:) )
	570	!
	571	END FUNCTION glob_max_vec_3d
	572
	573	FUNCTION glob_max_vec_4d( cdname, ptab ) RESULT( ptmp )
	574	!!----------------------------------------------------------------------
	575	CHARACTER(len=*), INTENT(in) :: cdname ! name of the calling subroutine
	576	REAL(wp), INTENT(in) :: ptab(:,:,:,:) ! array on which operation is applied
	577	REAL(wp), DIMENSION(SIZE(ptab,4)) :: ptmp
	578	!
	579	INTEGER :: jk , jl ! dummy loop indice & dimension
	580	INTEGER :: ipk, ipl ! dimension
	581	!!-----------------------------------------------------------------------
	582	!
	583	ipk = SIZE(ptab,3)
	584	ipl = SIZE(ptab,4)
	585	DO jl = 1, ipl
	586	ptmp(jl) = MAXVAL( ptab(:,:,1,jl) * tmask_i(:,:) )
	587	DO jk = 2, ipk
	588	ptmp(jl) = MAX( ptmp(jl), MAXVAL( ptab(:,:,jk,jl) * tmask_i(:,:) ) )
	589	ENDDO
	590	ENDDO
	591	!
	592	CALL mpp_max( cdname, ptmp (:) )
	593	!
	594	END FUNCTION glob_max_vec_4d
	595
[2003]	596	SUBROUTINE DDPDD( ydda, yddb )
	597	!!----------------------------------------------------------------------
	598	!! * ROUTINE DDPDD *
[3764]	599	!!
[2003]	600	!! ** Purpose : Add a scalar element to a sum
	601	!!
[3764]	602	!!
	603	!! ** Method : The code uses the compensated summation with doublet
[2003]	604	!! (sum,error) emulated useing complex numbers. ydda is the
[3764]	605	!! scalar to add to the summ yddb
[2003]	606	!!
[3764]	607	!! ** Action : This does only work for MPI.
	608	!!
[2003]	609	!! References : Using Acurate Arithmetics to Improve Numerical
	610	!! Reproducibility and Sability in Parallel Applications
[3764]	611	!! Yun HE and Chris H. Q. DING, Journal of Supercomputing 18, 259-277, 2001
[2003]	612	!!----------------------------------------------------------------------
[13226]	613	COMPLEX(dp), INTENT(in ) :: ydda
	614	COMPLEX(dp), INTENT(inout) :: yddb
[2307]	615	!
[13226]	616	REAL(dp) :: zerr, zt1, zt2 ! local work variables
[2307]	617	!!-----------------------------------------------------------------------
	618	!
[2003]	619	! Compute ydda + yddb using Knuth's trick.
[2307]	620	zt1 = REAL(ydda) + REAL(yddb)
	621	zerr = zt1 - REAL(ydda)
	622	zt2 = ( (REAL(yddb) - zerr) + (REAL(ydda) - (zt1 - zerr)) ) &
	623	& + AIMAG(ydda) + AIMAG(yddb)
	624	!
[2003]	625	! The result is t1 + t2, after normalization.
[2307]	626	yddb = CMPLX( zt1 + zt2, zt2 - ((zt1 + zt2) - zt1), wp )
	627	!
[2003]	628	END SUBROUTINE DDPDD
	629
	630	#if defined key_nosignedzero
[2307]	631	!!----------------------------------------------------------------------
	632	!! 'key_nosignedzero' F90 SIGN
	633	!!----------------------------------------------------------------------
[3764]	634
[2307]	635	FUNCTION SIGN_SCALAR( pa, pb )
[2003]	636	!!-----------------------------------------------------------------------
	637	!! * FUNCTION SIGN_SCALAR *
	638	!!
	639	!! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
	640	!!-----------------------------------------------------------------------
	641	REAL(wp) :: pa,pb ! input
[2307]	642	REAL(wp) :: SIGN_SCALAR ! result
	643	!!-----------------------------------------------------------------------
	644	IF ( pb >= 0.e0) THEN ; SIGN_SCALAR = ABS(pa)
	645	ELSE ; SIGN_SCALAR =-ABS(pa)
[2003]	646	ENDIF
	647	END FUNCTION SIGN_SCALAR
	648
[2307]	649
[3764]	650	FUNCTION SIGN_ARRAY_1D( pa, pb )
[2003]	651	!!-----------------------------------------------------------------------
	652	!! * FUNCTION SIGN_ARRAY_1D *
	653	!!
	654	!! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
	655	!!-----------------------------------------------------------------------
[2307]	656	REAL(wp) :: pa,pb(:) ! input
[2003]	657	REAL(wp) :: SIGN_ARRAY_1D(SIZE(pb,1)) ! result
[2307]	658	!!-----------------------------------------------------------------------
	659	WHERE ( pb >= 0.e0 ) ; SIGN_ARRAY_1D = ABS(pa)
	660	ELSEWHERE ; SIGN_ARRAY_1D =-ABS(pa)
[2003]	661	END WHERE
	662	END FUNCTION SIGN_ARRAY_1D
	663
[2307]	664
[3764]	665	FUNCTION SIGN_ARRAY_2D(pa,pb)
[2003]	666	!!-----------------------------------------------------------------------
	667	!! * FUNCTION SIGN_ARRAY_2D *
	668	!!
	669	!! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
	670	!!-----------------------------------------------------------------------
	671	REAL(wp) :: pa,pb(:,:) ! input
	672	REAL(wp) :: SIGN_ARRAY_2D(SIZE(pb,1),SIZE(pb,2)) ! result
[2307]	673	!!-----------------------------------------------------------------------
	674	WHERE ( pb >= 0.e0 ) ; SIGN_ARRAY_2D = ABS(pa)
	675	ELSEWHERE ; SIGN_ARRAY_2D =-ABS(pa)
[2003]	676	END WHERE
	677	END FUNCTION SIGN_ARRAY_2D
	678
[3764]	679	FUNCTION SIGN_ARRAY_3D(pa,pb)
[2003]	680	!!-----------------------------------------------------------------------
	681	!! * FUNCTION SIGN_ARRAY_3D *
	682	!!
	683	!! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
	684	!!-----------------------------------------------------------------------
	685	REAL(wp) :: pa,pb(:,:,:) ! input
	686	REAL(wp) :: SIGN_ARRAY_3D(SIZE(pb,1),SIZE(pb,2),SIZE(pb,3)) ! result
[2307]	687	!!-----------------------------------------------------------------------
	688	WHERE ( pb >= 0.e0 ) ; SIGN_ARRAY_3D = ABS(pa)
	689	ELSEWHERE ; SIGN_ARRAY_3D =-ABS(pa)
[2003]	690	END WHERE
	691	END FUNCTION SIGN_ARRAY_3D
	692
[2307]	693
[3764]	694	FUNCTION SIGN_ARRAY_1D_A(pa,pb)
[2003]	695	!!-----------------------------------------------------------------------
	696	!! * FUNCTION SIGN_ARRAY_1D_A *
	697	!!
	698	!! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
	699	!!-----------------------------------------------------------------------
	700	REAL(wp) :: pa(:),pb(:) ! input
[2307]	701	REAL(wp) :: SIGN_ARRAY_1D_A(SIZE(pb,1)) ! result
	702	!!-----------------------------------------------------------------------
	703	WHERE ( pb >= 0.e0 ) ; SIGN_ARRAY_1D_A = ABS(pa)
	704	ELSEWHERE ; SIGN_ARRAY_1D_A =-ABS(pa)
[2003]	705	END WHERE
	706	END FUNCTION SIGN_ARRAY_1D_A
	707
[2307]	708
[3764]	709	FUNCTION SIGN_ARRAY_2D_A(pa,pb)
[2003]	710	!!-----------------------------------------------------------------------
	711	!! * FUNCTION SIGN_ARRAY_2D_A *
	712	!!
	713	!! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
	714	!!-----------------------------------------------------------------------
	715	REAL(wp) :: pa(:,:),pb(:,:) ! input
	716	REAL(wp) :: SIGN_ARRAY_2D_A(SIZE(pb,1),SIZE(pb,2)) ! result
[2307]	717	!!-----------------------------------------------------------------------
	718	WHERE ( pb >= 0.e0 ) ; SIGN_ARRAY_2D_A = ABS(pa)
	719	ELSEWHERE ; SIGN_ARRAY_2D_A =-ABS(pa)
[2003]	720	END WHERE
	721	END FUNCTION SIGN_ARRAY_2D_A
	722
[2307]	723
[3764]	724	FUNCTION SIGN_ARRAY_3D_A(pa,pb)
[2003]	725	!!-----------------------------------------------------------------------
	726	!! * FUNCTION SIGN_ARRAY_3D_A *
	727	!!
	728	!! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
	729	!!-----------------------------------------------------------------------
	730	REAL(wp) :: pa(:,:,:),pb(:,:,:) ! input
	731	REAL(wp) :: SIGN_ARRAY_3D_A(SIZE(pb,1),SIZE(pb,2),SIZE(pb,3)) ! result
[2307]	732	!!-----------------------------------------------------------------------
	733	WHERE ( pb >= 0.e0 ) ; SIGN_ARRAY_3D_A = ABS(pa)
	734	ELSEWHERE ; SIGN_ARRAY_3D_A =-ABS(pa)
[2003]	735	END WHERE
	736	END FUNCTION SIGN_ARRAY_3D_A
	737
[2307]	738
[3764]	739	FUNCTION SIGN_ARRAY_1D_B(pa,pb)
[2003]	740	!!-----------------------------------------------------------------------
	741	!! * FUNCTION SIGN_ARRAY_1D_B *
	742	!!
	743	!! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
	744	!!-----------------------------------------------------------------------
	745	REAL(wp) :: pa(:),pb ! input
	746	REAL(wp) :: SIGN_ARRAY_1D_B(SIZE(pa,1)) ! result
[2307]	747	!!-----------------------------------------------------------------------
	748	IF( pb >= 0.e0 ) THEN ; SIGN_ARRAY_1D_B = ABS(pa)
	749	ELSE ; SIGN_ARRAY_1D_B =-ABS(pa)
[2003]	750	ENDIF
	751	END FUNCTION SIGN_ARRAY_1D_B
	752
[2307]	753
[3764]	754	FUNCTION SIGN_ARRAY_2D_B(pa,pb)
[2003]	755	!!-----------------------------------------------------------------------
	756	!! * FUNCTION SIGN_ARRAY_2D_B *
	757	!!
	758	!! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
	759	!!-----------------------------------------------------------------------
	760	REAL(wp) :: pa(:,:),pb ! input
	761	REAL(wp) :: SIGN_ARRAY_2D_B(SIZE(pa,1),SIZE(pa,2)) ! result
[2307]	762	!!-----------------------------------------------------------------------
	763	IF( pb >= 0.e0 ) THEN ; SIGN_ARRAY_2D_B = ABS(pa)
	764	ELSE ; SIGN_ARRAY_2D_B =-ABS(pa)
[2003]	765	ENDIF
	766	END FUNCTION SIGN_ARRAY_2D_B
	767
[2307]	768
[3764]	769	FUNCTION SIGN_ARRAY_3D_B(pa,pb)
[2003]	770	!!-----------------------------------------------------------------------
	771	!! * FUNCTION SIGN_ARRAY_3D_B *
	772	!!
	773	!! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
	774	!!-----------------------------------------------------------------------
	775	REAL(wp) :: pa(:,:,:),pb ! input
	776	REAL(wp) :: SIGN_ARRAY_3D_B(SIZE(pa,1),SIZE(pa,2),SIZE(pa,3)) ! result
[2307]	777	!!-----------------------------------------------------------------------
	778	IF( pb >= 0.e0 ) THEN ; SIGN_ARRAY_3D_B = ABS(pa)
	779	ELSE ; SIGN_ARRAY_3D_B =-ABS(pa)
[2003]	780	ENDIF
	781	END FUNCTION SIGN_ARRAY_3D_B
	782	#endif
	783
[2307]	784	!!======================================================================
[2003]	785	END MODULE lib_fortran

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

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

Context Navigation

source: NEMO/trunk/src/OCE/lib_fortran.F90 @ 15725

Download in other formats: