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.

math.f90 in utils/tools/SIREN/src – NEMO

Context Navigation

source: utils/tools/SIREN/src/math.f90 @ 12080

Last change on this file since 12080 was 12080, checked in by jpaul, 4 years ago
update nemo trunk
File size: 47.9 KB

Line
1	!----------------------------------------------------------------------
2	! NEMO system team, System and Interface for oceanic RElocable Nesting
3	!----------------------------------------------------------------------
4	!
5	! DESCRIPTION:
6	!> @brief
7	!> This module groups some useful mathematical function.
8	!>
9	!> @details
10	!>
11	!> to compute the mean of an array:<br/>
12	!> @code
13	!> dl_value=math_mean( dl_value, dd_fill )
14	!> @endcode
15	!> - dl_value is 1D or 2D array
16	!> - dd_fill is FillValue
17	!>
18	!> to compute the median of an array:<br/>
19	!> @code
20	!> dl_value=math_median( dl_value, dd_fill )
21	!> @endcode
22	!> - dl_value is 1D or 2D array
23	!> - dd_fill is FillValue
24	!>
25	!> to compute the mean without extremum of an array:<br/>
26	!> @code
27	!> dl_value=math_mwe( dl_value, id_next, dd_fill )
28	!> @endcode
29	!> - dl_value is 1D or 2D array
30	!> - id_next is the number of extremum to be removed
31	!> - dd_fill is FillValue
32	!>
33	!> to sort an 1D array:<br/>
34	!> @code
35	!> CALL math_QsortC(dl_value)
36	!> @endcode
37	!> - dl_value is 1D array
38	!>
39	!> to correct phase angles to produce smoother phase:<br/>
40	!> @code
41	!> CALL math_unwrap(dl_value, [dl_discont])
42	!> @endcode
43	!> - dl_value is 1D array
44	!> - dl_discont maximum discontinuity between values, default pi
45	!>
46	!> to compute simple operation
47	!> @code
48	!> dl_res=math_compute(cl_var)
49	!> @endcode
50	!> - cl_var operation to compute (string of character)
51	!> - dl_res result of the operation, real(dp)
52	!>
53	!> to compute first derivative of 1D array:<br/>
54	!> @code
55	!> dl_value(:)=math_deriv_1D( dd_value(:), dd_fill, [ld_discont] )
56	!> @endcode
57	!> - dd_value is 1D array of variable
58	!> - dd_fill is FillValue of variable
59	!> - ld_discont is logical to take into account longitudinal
60	!> East-West discontinuity [optional]
61	!>
62	!> to compute first derivative of 2D array:<br/>
63	!> @code
64	!> dl_value(:,:)=math_deriv_2D( dd_value(:,:), dd_fill, cd_dim,
65	!> [ld_discont] )
66	!> @endcode
67	!> - dd_value is 2D array of variable
68	!> - dd_fill is FillValue of variable
69	!> - cd_dim is character to compute derivative on first (I) or
70	!> second (J) dimension
71	!> - ld_discont is logical to take into account longitudinal
72	!> East-West discontinuity [optional]
73	!>
74	!> to compute first derivative of 3D array:<br/>
75	!> @code
76	!> dl_value(:,:,:)=math_deriv_3D( dd_value(:,:,:), dd_fill, cd_dim,
77	!> [ld_discont] )
78	!> @endcode
79	!> - dd_value is 3D array of variable
80	!> - dd_fill is FillValue of variable
81	!> - cd_dim is character to compute derivative on first (I), second (J),
82	!> or third (K) dimension
83	!> - ld_discont is logical to take into account longitudinal East-West
84	!> discontinuity [optional]
85	!>
86	!>
87	!> @author
88	!> J.Paul
89	!>
90	!> @date January, 2015 - Initial version
91	!>
92	!> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
93	!----------------------------------------------------------------------
94	MODULE math
95
96	USE kind ! F90 kind parameter
97	USE global ! global variable
98	USE phycst ! physical constant
99	USE fct ! basic useful function
100	IMPLICIT NONE
101	! NOTE_avoid_public_variables_if_possible
102
103	! function and subroutine
104	PUBLIC :: math_mean !< return mean of an array
105	PUBLIC :: math_median !< return median of an array
106	PUBLIC :: math_mwe !< return mean without extremum of an array
107	PUBLIC :: math_QsortC !< sort an 1D array
108	PUBLIC :: math_unwrap !< correct phase angles to produce smoother phase
109	PUBLIC :: math_compute !< compute simple operation
110	PUBLIC :: math_deriv_1D !< compute first derivative of 1D array
111	PUBLIC :: math_deriv_2D !< compute first derivative of 2D array
112	PUBLIC :: math_deriv_3D !< compute first derivative of 3D array
113	PUBLIC :: math_ortho !< compute orthodome distance
114	PUBLIC :: math_euclid !< compute euclidian distance
115
116	PRIVATE :: math__Partition
117	PRIVATE :: math__mean_1d
118	PRIVATE :: math__mean_2d
119	PRIVATE :: math__median_1d
120	PRIVATE :: math__median_2d
121	PRIVATE :: math__mwe_1d
122	PRIVATE :: math__mwe_2d
123	PRIVATE :: math__parentheses
124
125
126	INTERFACE math_mean
127	MODULE PROCEDURE math__mean_1d ! return mean of an array 1D
128	MODULE PROCEDURE math__mean_2d ! return mean of an array 2D
129	END INTERFACE math_mean
130
131	INTERFACE math_median
132	MODULE PROCEDURE math__median_1d ! return median of an array 1D
133	MODULE PROCEDURE math__median_2d ! return median of an array 2D
134	END INTERFACE math_median
135
136	INTERFACE math_mwe
137	MODULE PROCEDURE math__mwe_1d ! return mean without extremum of an array 1D
138	MODULE PROCEDURE math__mwe_2d ! return mean without extremum of an array 2D
139	END INTERFACE math_mwe
140
141	CONTAINS
142	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
143	PURE FUNCTION math__mean_1d(dd_array, dd_fill) &
144	& RESULT (df_mean)
145	!-------------------------------------------------------------------
146	!> @brief This function compute the mean of a 1D array.
147	!>
148	!> @author J.Paul
149	!> @date January, 2015 - Initial Version
150	!>
151	!> @param[in] dd_array 1D array
152	!> @param[in] dd_fill fillValue
153	!> @return mean value, real(dp)
154	!-------------------------------------------------------------------
155
156	IMPLICIT NONE
157
158	! Argument
159	REAL(dp), DIMENSION(:), INTENT(IN) :: dd_array
160	REAL(dp), INTENT(IN), OPTIONAL :: dd_fill
161
162	! function
163	REAL(dp) :: df_mean
164
165	! local variable
166	INTEGER(i4) :: il_count
167	REAL(dp) :: dl_sum
168	REAL(dp) :: dl_count
169	!----------------------------------------------------------------
170
171	IF( PRESENT(dd_fill) )THEN
172	il_count=COUNT(dd_array(:)/=dd_fill)
173	IF( il_count > 0 )THEN
174	dl_sum =SUM ( dd_array(:), dd_array(:)/= dd_fill )
175	dl_count=REAL( il_count, dp)
176
177	df_mean=dl_sum/dl_count
178	ELSE
179	df_mean=dd_fill
180	ENDIF
181	ELSE
182	il_count=SIZE(dd_array(:))
183	IF( il_count > 0 )THEN
184	dl_sum =SUM ( dd_array(:) )
185	dl_count=REAL( il_count, dp)
186
187	df_mean=dl_sum/dl_count
188	ELSE
189	df_mean=0
190	ENDIF
191	ENDIF
192
193	END FUNCTION math__mean_1d
194	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
195	PURE FUNCTION math__mean_2d(dd_array, dd_fill) &
196	& RESULT (df_mean)
197	!-------------------------------------------------------------------
198	!> @brief This function compute the mean of a 2D array.
199	!>
200	!> @author J.Paul
201	!> @date January, 2015 - Initial Version
202	!>
203	!> @param[in] dd_array 2D array
204	!> @param[in] dd_fill fillValue
205	!> @return mean value, real(dp)
206	!-------------------------------------------------------------------
207
208	IMPLICIT NONE
209
210	! Argument
211	REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_array
212	REAL(dp), INTENT(IN), OPTIONAL :: dd_fill
213
214	! function
215	REAL(dp) :: df_mean
216
217	! local variable
218	INTEGER(i4) :: il_count
219
220	REAL(dp), DIMENSION(:), ALLOCATABLE :: dl_list
221	!----------------------------------------------------------------
222
223	IF( PRESENT(dd_fill) )THEN
224	il_count=COUNT(dd_array(:,:)/=dd_fill)
225	IF( il_count > 0 )THEN
226	ALLOCATE( dl_list(il_count) )
227	dl_list(:)=PACK(dd_array(:,:),dd_array(:,:)/=dd_fill)
228	ELSE
229	df_mean=dd_fill
230	ENDIF
231	ELSE
232	il_count=SIZE(dd_array)
233	IF( il_count > 0 )THEN
234	ALLOCATE( dl_list(il_count) )
235	dl_list(:)=PACK(dd_array(:,:), MASK=.TRUE.)
236	ELSE
237	df_mean=0
238	ENDIF
239	ENDIF
240
241	IF( ALLOCATED(dl_list) )THEN
242	df_mean=math_mean(dl_list(:))
243	DEALLOCATE( dl_list )
244	ENDIF
245
246	END FUNCTION math__mean_2d
247	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
248	PURE FUNCTION math__median_1d(dd_array, dd_fill) &
249	& RESULT (df_median)
250	!-------------------------------------------------------------------
251	!> @brief This function compute the median of a 1D array.
252	!>
253	!> @author J.Paul
254	!> @date January, 2015 - Initial Version
255	!>
256	!> @param[in] dd_array 1D array
257	!> @param[in] dd_fill fillValue
258	!> @return median value, real(dp)
259	!-------------------------------------------------------------------
260
261	IMPLICIT NONE
262
263	! Argument
264	REAL(dp), DIMENSION(:), INTENT(IN) :: dd_array
265	REAL(dp), INTENT(IN), OPTIONAL :: dd_fill
266
267	! function
268	REAL(dp) :: df_median
269
270	! local variable
271	INTEGER(i4) :: il_count
272
273	REAL(dp), DIMENSION(:), ALLOCATABLE :: dl_list
274	!----------------------------------------------------------------
275
276	IF( PRESENT(dd_fill) )THEN
277	il_count=COUNT(dd_array(:)/=dd_fill)
278	IF( il_count > 0 )THEN
279	ALLOCATE( dl_list(il_count) )
280	dl_list(:)=PACK(dd_array(:),dd_array(:)/=dd_fill)
281	ELSE
282	df_median=dd_fill
283	ENDIF
284	ELSE
285	il_count=SIZE(dd_array(:))
286	IF( il_count > 0 )THEN
287	ALLOCATE( dl_list(il_count) )
288	dl_list(:)=dd_array(:)
289	ELSE
290	df_median=0
291	ENDIF
292	ENDIF
293
294	IF( ALLOCATED(dl_list) )THEN
295	CALL math_QsortC(dl_list(:))
296
297	IF( MOD(il_count,2) == 0 )THEN
298	df_median=(dl_list(il_count/2+1)+dl_list(il_count/2))/2_dp
299	ELSE
300	df_median=dl_list(il_count/2+1)
301	ENDIF
302
303	DEALLOCATE(dl_list)
304	ENDIF
305
306	END FUNCTION math__median_1d
307	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
308	PURE FUNCTION math__median_2d(dd_array, dd_fill) &
309	& RESULT (df_median)
310	!-------------------------------------------------------------------
311	!> @brief This function compute the median of a 2D array.
312	!>
313	!> @author J.Paul
314	!> @date January, 2015 - Initial Version
315	!>
316	!> @param[in] dd_array 2D array
317	!> @param[in] dd_fill fillValue
318	!> @return median value, real(dp)
319	!-------------------------------------------------------------------
320
321	IMPLICIT NONE
322
323	! Argument
324	REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_array
325	REAL(dp), INTENT(IN), OPTIONAL :: dd_fill
326
327	! funtion
328	REAL(dp) :: df_median
329
330	! local variable
331	INTEGER(i4) :: il_count
332
333	REAL(dp), DIMENSION(:), ALLOCATABLE :: dl_list
334	!----------------------------------------------------------------
335
336	IF( PRESENT(dd_fill) )THEN
337	il_count=COUNT(dd_array(:,:)/=dd_fill)
338	IF( il_count > 0 )THEN
339	ALLOCATE( dl_list(il_count) )
340	dl_list(:)=PACK(dd_array(:,:),dd_array(:,:)/=dd_fill)
341	ELSE
342	df_median=dd_fill
343	ENDIF
344	ELSE
345	il_count=SIZE(dd_array(:,:))
346	IF( il_count > 0 )THEN
347	ALLOCATE( dl_list(il_count) )
348	dl_list(:)=PACK(dd_array(:,:), MASK=.TRUE.)
349	ELSE
350	df_median=0
351	ENDIF
352	ENDIF
353
354	IF( ALLOCATED(dl_list) )THEN
355	df_median=math_median(dl_list(:))
356	DEALLOCATE(dl_list)
357	ENDIF
358
359	END FUNCTION math__median_2d
360	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
361	PURE FUNCTION math__mwe_1d(dd_array, id_next, dd_fill) &
362	& RESULT (df_mwe)
363	!-------------------------------------------------------------------
364	!> @brief This function compute the mean without extremum of a 1D array.
365	!>
366	!> @author J.Paul
367	!> @date January, 2015 - Initial Version
368	!>
369	!> @param[in] dd_array 1D array
370	!> @param[in] id_next number of extremum to be removed
371	!> @param[in] dd_fill fillValue
372	!> @return median value, real(dp)
373	!-------------------------------------------------------------------
374
375	IMPLICIT NONE
376
377	! Argument
378	REAL(dp), DIMENSION(:), INTENT(IN) :: dd_array
379	INTEGER(i4) , INTENT(IN), OPTIONAL :: id_next
380	REAL(dp), INTENT(IN), OPTIONAL :: dd_fill
381
382	! function
383	REAL(dp) :: df_mwe
384
385	! local variable
386	INTEGER(i4) :: il_next
387	INTEGER(i4) :: il_count
388	INTEGER(i4) :: il_size
389
390	REAL(dp), DIMENSION(:), ALLOCATABLE :: dl_list
391	!----------------------------------------------------------------
392
393	il_next=2
394	IF( PRESENT(id_next) ) il_next=id_next
395
396	il_size=SIZE(dd_array(:))
397	IF( PRESENT(dd_fill) )THEN
398	il_count=COUNT(dd_array(:)/=dd_fill)
399	IF( il_count > 0 )THEN
400	ALLOCATE( dl_list(il_count) )
401	dl_list(:)=PACK(dd_array(:),dd_array(:)/=dd_fill)
402	ELSE
403	df_mwe=dd_fill
404	ENDIF
405	ELSE
406	il_count=SIZE(dd_array(:))
407	IF( il_count > 0 )THEN
408	ALLOCATE( dl_list(il_count) )
409	dl_list(:)=dd_array(:)
410	ELSE
411	df_mwe=0
412	ENDIF
413	ENDIF
414
415	IF( ALLOCATED(dl_list) )THEN
416	CALL math_QsortC(dl_list(:))
417
418	IF( il_count == il_size )THEN
419	! no fillValue
420	df_mwe=math_mean(dl_list(il_next+1:il_size-il_next))
421	ELSEIF( il_count > il_size-2*il_next )THEN
422	! remove one extremum each side
423	df_mwe=math_mean(dl_list(2:il_size-1))
424	ELSE ! il_count <= il_size-2*il_next
425	! more than 2*il_next fillValue
426	! compute mean only
427	df_mwe=math_mean(dl_list(:))
428	ENDIF
429
430	DEALLOCATE(dl_list)
431	ENDIF
432
433	END FUNCTION math__mwe_1d
434	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
435	PURE FUNCTION math__mwe_2d(dd_array, id_next, dd_fill) &
436	& RESULT (df_mwe)
437	!-------------------------------------------------------------------
438	!> @brief This function compute the mean without extremum of a 2D array.
439	!>
440	!> @author J.Paul
441	!> @date January, 2015 - Initial Version
442	!>
443	!> @param[in] dd_array 2D array
444	!> @param[in] id_next number of extremum to be removed
445	!> @param[in] dd_fill fillValue
446	!> @return median value, real(dp)
447	!-------------------------------------------------------------------
448
449	IMPLICIT NONE
450
451	! Argument
452	REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_array
453	INTEGER(i4) , INTENT(IN), OPTIONAL :: id_next
454	REAL(dp), INTENT(IN), OPTIONAL :: dd_fill
455
456	! function
457	REAL(dp) :: df_mwe
458
459	! local variable
460	INTEGER(i4) :: il_count
461
462	REAL(dp), DIMENSION(:), ALLOCATABLE :: dl_list
463	!----------------------------------------------------------------
464
465	IF( PRESENT(dd_fill) )THEN
466	il_count=COUNT(dd_array(:,:)/=dd_fill)
467	IF( il_count > 0 )THEN
468	ALLOCATE( dl_list(il_count) )
469	dl_list(:)=PACK(dd_array(:,:),dd_array(:,:)/=dd_fill)
470
471	df_mwe=math_mwe(dl_list(:), id_next)
472	ELSE
473	df_mwe=dd_fill
474	ENDIF
475	ELSE
476	il_count=SIZE(dd_array(:,:))
477	IF( il_count > 0 )THEN
478	ALLOCATE( dl_list(il_count) )
479	dl_list(:)=PACK(dd_array(:,:), MASK=.TRUE.)
480
481	df_mwe=math_mwe(dl_list(:), id_next)
482	ELSE
483	df_mwe=0
484	ENDIF
485	ENDIF
486
487	IF( ALLOCATED(dl_list) ) DEALLOCATE(dl_list)
488
489	END FUNCTION math__mwe_2d
490	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
491	PURE RECURSIVE SUBROUTINE math_QsortC(dd_array)
492	!-------------------------------------------------------------------
493	!> @brief This subroutine sort a 1D array.
494	!>
495	!> @details
496	!> Recursive Fortran 95 quicksort routine
497	!> sorts real numbers into ascending numerical order
498	!> Author: Juli Rew, SCD Consulting (juliana@ucar.edu), 9/03
499	!> Based on algorithm from Cormen et al., Introduction to Algorithms,
500	!> 1997 printing
501	!>
502	!> @author J.Paul
503	!> @date January, 2015 - Rewrite with SIREN coding rules
504	!>
505	!> @param[inout] dd_array 1D array
506	!-------------------------------------------------------------------
507
508	IMPLICIT NONE
509
510	! Argument
511	REAL(dp), DIMENSION(:), INTENT(INOUT) :: dd_array
512
513	! local variable
514	INTEGER(i4) :: il_iq
515	!----------------------------------------------------------------
516
517	IF( SIZE(dd_array(:)) > 1 )THEN
518	CALL math__Partition(dd_array, il_iq)
519	CALL math_QsortC(dd_array(:il_iq-1))
520	CALL math_QsortC(dd_array(il_iq:))
521	ENDIF
522
523	END SUBROUTINE math_QsortC
524	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
525	PURE SUBROUTINE math__Partition(dd_array, id_marker)
526	!-------------------------------------------------------------------
527	!> @brief This subroutine partition a 1D array.
528	!>
529	!> @details
530	!> Author: Juli Rew, SCD Consulting (juliana@ucar.edu), 9/03
531	!> Based on algorithm from Cormen et al., Introduction to Algorithms,
532	!> 1997 printing
533	!>
534	!> @author J.Paul
535	!> @date January, 2015 - Rewrite with SIREN coding rules
536	!> @date November, 2017
537	!> - use the correct loop index to look for element bigger than pivot point.
538	!>
539	!> @param[inout] dd_array 1D array
540	!> @param[in] id_marker
541	!-------------------------------------------------------------------
542
543	IMPLICIT NONE
544
545	! Argument
546	REAL(dp) , DIMENSION(:), INTENT(INOUT) :: dd_array
547	INTEGER(i4), INTENT( OUT) :: id_marker
548
549	! local variable
550	REAL(dp) :: dl_temp
551	REAL(dp) :: dl_x ! pivot point
552
553	! loop indices
554	INTEGER(i4) :: ji
555	INTEGER(i4) :: jj
556	!----------------------------------------------------------------
557
558	dl_x = dd_array(1)
559	ji= 0
560	jj= SIZE(dd_array(:)) + 1
561
562	DO
563	jj=jj-1
564	DO
565	IF( dd_array(jj) <= dl_x ) EXIT
566	jj=jj-1
567	ENDDO
568	ji=ji+1
569	DO
570	IF( dd_array(ji) >= dl_x ) EXIT
571	ji=ji+1
572	ENDDO
573	IF( ji < jj )THEN
574	! exchange dd_array(ji) and dd_array(jj)
575	dl_temp= dd_array(ji)
576	dd_array(ji) = dd_array(jj)
577	dd_array(jj) = dl_temp
578	ELSEIF( ji==jj )THEN
579	id_marker=ji+1
580	RETURN
581	ELSE
582	id_marker=ji
583	RETURN
584	ENDIF
585	ENDDO
586
587	END SUBROUTINE math__Partition
588	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
589	PURE SUBROUTINE math_unwrap(dd_array, dd_discont)
590	!-------------------------------------------------------------------
591	!> @brief This subroutine correct phase angles to produce smoother
592	!> phase plots.
593	!>
594	!> @details
595	!> This code is based on numpy unwrap function
596	!>
597	!> Unwrap by changing deltas between values to 2*pi complement.
598	!>
599	!> Unwrap radian phase `dd_array` by changing absolute jumps greater than
600	!> `dd_discont` to their 2*pi complement.
601	!>
602	!> @note If the discontinuity in `dd_array` is smaller than ``pi``,
603	!> but larger than `dd_discont`, no unwrapping is done because taking
604	!> the 2*pi complement would only make the discontinuity larger.
605	!>
606	!> @author J.Paul
607	!> @date Marsh, 2015 - Rewrite in fortran, with SIREN coding rules
608	!>
609	!> @param[inout] dd_array 1D array
610	!> @param[in] dd_discont maximum discontinuity between values, default pi
611	!-------------------------------------------------------------------
612
613	IMPLICIT NONE
614
615	! Argument
616	REAL(dp) , DIMENSION(:), INTENT(INOUT) :: dd_array
617	REAL(dp) , INTENT(IN ), OPTIONAL :: dd_discont
618
619	! local variable
620	INTEGER(i4) :: il_size
621
622	REAL(dp) :: dl_discont
623
624	REAL(dp), DIMENSION(:), ALLOCATABLE :: dl_diff
625	REAL(dp), DIMENSION(:), ALLOCATABLE :: dl_mod
626	REAL(dp), DIMENSION(:), ALLOCATABLE :: dl_correct
627	REAL(dp), DIMENSION(:), ALLOCATABLE :: dl_tmp
628
629	! loop indices
630	INTEGER(i4) :: ji
631	!----------------------------------------------------------------
632	dl_discont=dp_pi
633	IF( PRESENT(dd_discont) ) dl_discont=dd_discont
634
635	il_size=SIZE(dd_array)
636	ALLOCATE(dl_diff(il_size-1))
637	DO ji=1,il_size-1
638	dl_diff(ji)=dd_array(ji+1)-dd_array(ji)
639	ENDDO
640
641	ALLOCATE(dl_mod(il_size-1))
642	DO ji=1,il_size-1
643	dl_mod(ji) = MOD(dl_diff(ji) + dp_pi, 2*dp_pi) - dp_pi
644	ENDDO
645
646	WHERE( (dl_mod(:) == -dp_pi) .AND. (dl_diff(:) > 0._dp ) )
647	dl_mod(:)=dp_pi
648	END WHERE
649
650	ALLOCATE(dl_correct(il_size-1))
651	dl_correct(:)=dl_mod(:)-dl_diff(:)
652
653	DEALLOCATE(dl_mod)
654
655	WHERE( ABS(dl_diff(:)) < dl_discont )
656	dl_correct(:)=0._dp
657	END WHERE
658
659	DEALLOCATE(dl_diff)
660
661	ALLOCATE(dl_tmp(il_size))
662	dl_tmp(:)=dd_array(:)
663
664	DO ji=1,il_size-1
665	dd_array(ji+1)=dl_tmp(ji+1)+SUM(dl_correct(1:ji))
666	ENDDO
667
668	DEALLOCATE(dl_correct)
669	DEALLOCATE(dl_tmp)
670
671	END SUBROUTINE math_unwrap
672	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
673	RECURSIVE FUNCTION math_compute(cd_var) &
674	& RESULT(df_res)
675	!-------------------------------------------------------------------
676	!> @brief This function compute simple operation
677	!>
678	!> @details
679	!> - operation should be write as a string of character.
680	!> - operators allowed are : +,-,*,/
681	!> - to ordered operation you should use parentheses
682	!>
683	!> exemples: '1e6/(16/122)', '(3/2)*(2+1)'
684	!>
685	!> @author J.Paul
686	!> @date June, 2015 - initial version
687	!>
688	!> @param[in] cd_var operation to compute (string of character)
689	!> @return result of the operation, real(dp)
690	!-------------------------------------------------------------------
691
692	IMPLICIT NONE
693
694	! Argument
695	CHARACTER(LEN=*), INTENT(IN) :: cd_var
696
697	! fucntion
698	REAL(dp) :: df_res
699
700	! local variables
701	CHARACTER(LEN=lc) :: cl_var
702	CHARACTER(LEN=lc) :: cl_str1
703	CHARACTER(LEN=lc) :: cl_str2
704
705	INTEGER(i4) :: il_ind
706	! loop indices
707	!----------------------------------------------------------------
708
709
710	IF(fct_is_real(cd_var))THEN
711	READ(cd_var,*) df_res
712	ELSE
713
714
715	CALL math__parentheses(cd_var, cl_var)
716
717	IF(fct_is_real(cl_var))THEN
718	READ(cl_var,*) df_res
719	ELSE
720	il_ind=SCAN(TRIM(cl_var),'*')
721	IF( il_ind /= 0 )THEN
722	cl_str1=cl_var(1:il_ind-1)
723	cl_str2=cl_var(il_ind+1:)
724	df_res=math_compute(cl_str1)*math_compute(cl_str2)
725	ELSE
726	il_ind=SCAN(TRIM(cl_var),'/')
727	IF( il_ind /= 0 )THEN
728	cl_str1=cl_var(1:il_ind-1)
729	cl_str2=cl_var(il_ind+1:)
730	df_res=math_compute(cl_str1)/math_compute(cl_str2)
731	ELSE
732	il_ind=SCAN(TRIM(cl_var),'+')
733	IF( il_ind /= 0 )THEN
734	cl_str1=cl_var(1:il_ind-1)
735	cl_str2=cl_var(il_ind+1:)
736	df_res=math_compute(cl_str1)+math_compute(cl_str2)
737	ELSE
738	il_ind=SCAN(TRIM(cl_var),'-')
739	IF( il_ind /= 0 )THEN
740	cl_str1=cl_var(1:il_ind-1)
741	cl_str2=cl_var(il_ind+1:)
742	df_res=math_compute(cl_str1)-math_compute(cl_str2)
743	ELSE
744	df_res=dp_fill
745	ENDIF
746	ENDIF
747	ENDIF
748	ENDIF
749	ENDIF
750	ENDIF
751
752	END FUNCTION math_compute
753	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
754	SUBROUTINE math__parentheses(cd_varin, cd_varout)
755	!-------------------------------------------------------------------
756	!> @brief This subroutine replace sub string inside parentheses
757	!> by the value of the operation inside.
758	!>
759	!> @details
760	!> exemple :
761	!> - '2.6+(3/2)' => '2.6+1.5000'
762	!>
763	!> @author J.Paul
764	!> @date June, 2015 - initial version
765	!>
766	!> @param[in] cd_varin string of character with operation inside
767	!> parentheses
768	!> @param[out] cd_varout string of character with result of
769	!> operation inside parentheses
770	!-------------------------------------------------------------------
771
772	IMPLICIT NONE
773
774	! Argument
775	CHARACTER(LEN=*) , INTENT(IN) :: cd_varin
776	CHARACTER(LEN=lc), INTENT(OUT) :: cd_varout
777
778	! local variables
779	CHARACTER(LEN=lc) :: cl_cpt
780	INTEGER(i4) :: il_ind
781	INTEGER(i4) :: il_count
782
783	! loop indices
784	INTEGER(i4) :: ji
785	!----------------------------------------------------------------
786	il_ind=INDEX(cd_varin,'(')
787	IF( il_ind /= 0 )THEN
788	il_count=0
789	DO ji=il_ind+1,LEN(cd_varin)
790	IF( cd_varin(ji:ji) == '(' )THEN
791	il_count=il_count+1
792	ELSEIF( cd_varin(ji:ji) == ')' )THEN
793	IF( il_count == 0 )THEN
794	WRITE(cl_cpt,*) math_compute(cd_varin(il_ind+1:ji-1))
795	cd_varout=TRIM(cd_varin(1:il_ind-1))//TRIM(ADJUSTL(cl_cpt))//&
796	& TRIM(cd_varin(ji+1:))
797	EXIT
798	ELSE
799	il_count=il_count-1
800	ENDIF
801	ENDIF
802	ENDDO
803	ELSE
804	cd_varout=TRIM(cd_varin)
805	ENDIF
806
807	END SUBROUTINE math__parentheses
808	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
809	PURE FUNCTION math_deriv_1D(dd_value, dd_fill, ld_discont) &
810	& RESULT (df_deriv)
811	!-------------------------------------------------------------------
812	!> @brief
813	!> This function compute derivative of 1D array.
814	!>
815	!> @details
816	!> optionaly you could specify to take into account east west discontinuity
817	!> (-180° 180° or 0° 360° for longitude variable)
818	!>
819	!> @author J.Paul
820	!> @date November, 2013 - Initial Version
821	!>
822	!> @param[in] dd_value 1D array of variable to be extrapolated
823	!> @param[in] dd_fill FillValue of variable
824	!> @param[in] ld_discont logical to take into account east west discontinuity
825	!-------------------------------------------------------------------
826
827	IMPLICIT NONE
828
829	! Argument
830	REAL(dp) , DIMENSION(:), INTENT(IN) :: dd_value
831	REAL(dp) , INTENT(IN) :: dd_fill
832	LOGICAL , INTENT(IN), OPTIONAL :: ld_discont
833
834	! function
835	REAL(dp), DIMENSION(SIZE(dd_value,DIM=1) ) :: df_deriv
836
837	! local variable
838	INTEGER(i4) :: il_imin
839	INTEGER(i4) :: il_imax
840	INTEGER(i4), DIMENSION(1) :: il_shape
841
842	REAL(dp) :: dl_min
843	REAL(dp) :: dl_max
844	REAL(dp) , DIMENSION(:), ALLOCATABLE :: dl_value
845
846	LOGICAL :: ll_discont
847
848	! loop indices
849	INTEGER(i4) :: ji
850
851	INTEGER(i4) :: i1
852	INTEGER(i4) :: i2
853	!----------------------------------------------------------------
854	! init
855	df_deriv(:)=dd_fill
856
857	ll_discont=.FALSE.
858	IF( PRESENT(ld_discont) ) ll_discont=ld_discont
859
860	il_shape(:)=SHAPE(dd_value(:))
861
862	ALLOCATE( dl_value(3))
863
864	! compute derivative in i-direction
865	DO ji=1,il_shape(1)
866
867	il_imin=MAX(ji-1,1)
868	il_imax=MIN(ji+1,il_shape(1))
869
870	IF( il_imin==ji-1 .AND. il_imax==ji+1 )THEN
871	i1=1 ; i2=3
872	ELSEIF( il_imin==ji .AND. il_imax==ji+1 )THEN
873	i1=1 ; i2=2
874	ELSEIF( il_imin==ji-1 .AND. il_imax==ji )THEN
875	i1=2 ; i2=3
876	ENDIF
877
878	dl_value(i1:i2)=dd_value(il_imin:il_imax)
879	IF( il_imin == 1 )THEN
880	dl_value(:)=EOSHIFT( dl_value(:), &
881	& DIM=1, &
882	& SHIFT=-1, &
883	& BOUNDARY=dl_value(1) )
884	ENDIF
885	IF( il_imax == il_shape(1) )THEN
886	dl_value(:)=EOSHIFT( dl_value(:), &
887	& DIM=1, &
888	& SHIFT=1, &
889	& BOUNDARY=dl_value(3))
890	ENDIF
891
892	IF( ll_discont )THEN
893	dl_min=MINVAL( dl_value(:), dl_value(:)/=dd_fill )
894	dl_max=MAXVAL( dl_value(:), dl_value(:)/=dd_fill )
895	IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
896	WHERE( dl_value(:) < 0._dp )
897	dl_value(:) = dl_value(:)+360._dp
898	END WHERE
899	ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
900	WHERE( dl_value(:) > 180._dp )
901	dl_value(:) = dl_value(:)-180._dp
902	END WHERE
903	ENDIF
904	ENDIF
905
906	IF( dl_value( 2) /= dd_fill .AND. & ! ji
907	& dl_value( 3) /= dd_fill .AND. & ! ji+1
908	& dl_value( 1) /= dd_fill )THEN ! ji-1
909
910	df_deriv(ji)= (dl_value(3) - dl_value(1)) / REAL(il_imax-il_imin,dp)
911
912	ENDIF
913
914	ENDDO
915
916	DEALLOCATE( dl_value )
917
918	END FUNCTION math_deriv_1D
919	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
920	FUNCTION math_deriv_2D(dd_value, dd_fill, cd_dim, ld_discont) &
921	& RESULT (df_deriv)
922	!-------------------------------------------------------------------
923	!> @brief
924	!> This function compute derivative of 2D array.
925	!> you have to specify in which direction derivative have to be computed:
926	!> first (I) or second (J) dimension.
927	!>
928	!> @details
929	!> optionaly you could specify to take into account east west discontinuity
930	!> (-180° 180° or 0° 360° for longitude variable)
931	!>
932	!> @author J.Paul
933	!> @date November, 2013 - Initial Version
934	!>
935	!> @param[in] dd_value 2D array of variable to be extrapolated
936	!> @param[in] dd_fill FillValue of variable
937	!> @param[in] cd_dim compute derivative on first (I) or second (J) dimension
938	!> @param[in] ld_discont logical to take into account east west discontinuity
939	!-------------------------------------------------------------------
940
941	IMPLICIT NONE
942
943	! Argument
944	REAL(dp) , DIMENSION(:,:), INTENT(IN) :: dd_value
945	REAL(dp) , INTENT(IN) :: dd_fill
946	CHARACTER(LEN=*) , INTENT(IN) :: cd_dim
947	LOGICAL , INTENT(IN), OPTIONAL :: ld_discont
948
949	! function
950	REAL(dp), DIMENSION(SIZE(dd_value,DIM=1), &
951	& SIZE(dd_value,DIM=2) ) :: df_deriv
952
953	! local variable
954	INTEGER(i4) :: il_imin
955	INTEGER(i4) :: il_imax
956	INTEGER(i4) :: il_jmin
957	INTEGER(i4) :: il_jmax
958	INTEGER(i4), DIMENSION(2) :: il_shape
959
960	REAL(dp) :: dl_min
961	REAL(dp) :: dl_max
962	REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_value
963
964	LOGICAL :: ll_discont
965
966	! loop indices
967	INTEGER(i4) :: ji
968	INTEGER(i4) :: jj
969
970	INTEGER(i4) :: i1
971	INTEGER(i4) :: i2
972
973	INTEGER(i4) :: j1
974	INTEGER(i4) :: j2
975	!----------------------------------------------------------------
976	! init
977	df_deriv(:,:)=dd_fill
978
979	ll_discont=.FALSE.
980	IF( PRESENT(ld_discont) ) ll_discont=ld_discont
981
982	il_shape(:)=SHAPE(dd_value(:,:))
983
984	SELECT CASE(TRIM(fct_upper(cd_dim)))
985
986	CASE('I')
987
988	ALLOCATE( dl_value(3,il_shape(2)) )
989	! compute derivative in i-direction
990	DO ji=1,il_shape(1)
991
992	! init
993	dl_value(:,:)=dd_fill
994
995	il_imin=MAX(ji-1,1)
996	il_imax=MIN(ji+1,il_shape(1))
997
998	IF( il_imin==ji-1 .AND. il_imax==ji+1 )THEN
999	i1=1 ; i2=3
1000	ELSEIF( il_imin==ji .AND. il_imax==ji+1 )THEN
1001	i1=1 ; i2=2
1002	ELSEIF( il_imin==ji-1 .AND. il_imax==ji )THEN
1003	i1=2 ; i2=3
1004	ENDIF
1005
1006	dl_value(i1:i2,:)=dd_value(il_imin:il_imax,:)
1007	IF( il_imin == 1 )THEN
1008	dl_value(:,:)=EOSHIFT( dl_value(:,:), &
1009	& DIM=1, &
1010	& SHIFT=-1, &
1011	& BOUNDARY=dl_value(1,:) )
1012	ENDIF
1013	IF( il_imax == il_shape(1) )THEN
1014	dl_value(:,:)=EOSHIFT( dl_value(:,:), &
1015	& DIM=1, &
1016	& SHIFT=1, &
1017	& BOUNDARY=dl_value(3,:))
1018	ENDIF
1019
1020	IF( ll_discont )THEN
1021	dl_min=MINVAL( dl_value(:,:), dl_value(:,:)/=dd_fill )
1022	dl_max=MAXVAL( dl_value(:,:), dl_value(:,:)/=dd_fill )
1023	IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
1024	WHERE( dl_value(:,:) < 0_dp )
1025	dl_value(:,:) = dl_value(:,:)+360._dp
1026	END WHERE
1027	ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
1028	WHERE( dl_value(:,:) > 180 )
1029	dl_value(:,:) = dl_value(:,:)-180._dp
1030	END WHERE
1031	ENDIF
1032	ENDIF
1033
1034	WHERE( dl_value(2,:) /= dd_fill .AND. & ! ji
1035	& dl_value(3,:) /= dd_fill .AND. & ! ji+1
1036	& dl_value(1,:) /= dd_fill ) ! ji-1
1037
1038	df_deriv(ji,:)= (dl_value(3,:) - dl_value(1,:)) / REAL(il_imax-il_imin,dp)
1039
1040	END WHERE
1041
1042	ENDDO
1043
1044	CASE('J')
1045
1046	ALLOCATE( dl_value(il_shape(1),3) )
1047	! compute derivative in j-direction
1048	DO jj=1,il_shape(2)
1049
1050	il_jmin=MAX(jj-1,1)
1051	il_jmax=MIN(jj+1,il_shape(2))
1052
1053	IF( il_jmin==jj-1 .AND. il_jmax==jj+1 )THEN
1054	j1=1 ; j2=3
1055	ELSEIF( il_jmin==jj .AND. il_jmax==jj+1 )THEN
1056	j1=1 ; j2=2
1057	ELSEIF( il_jmin==jj-1 .AND. il_jmax==jj )THEN
1058	j1=2 ; j2=3
1059	ENDIF
1060
1061	dl_value(:,j1:j2)=dd_value(:,il_jmin:il_jmax)
1062	IF( il_jmin == 1 )THEN
1063	dl_value(:,:)=EOSHIFT( dl_value(:,:), &
1064	& DIM=2, &
1065	& SHIFT=-1, &
1066	& BOUNDARY=dl_value(:,1))
1067	ENDIF
1068	IF( il_jmax == il_shape(2) )THEN
1069	dl_value(:,:)=EOSHIFT( dl_value(:,:), &
1070	& DIM=2, &
1071	& SHIFT=1, &
1072	& BOUNDARY=dl_value(:,3))
1073	ENDIF
1074
1075	IF( ll_discont )THEN
1076	dl_min=MINVAL( dl_value(:,:), dl_value(:,:)/=dd_fill )
1077	dl_max=MAXVAL( dl_value(:,:), dl_value(:,:)/=dd_fill )
1078	IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
1079	WHERE( dl_value(:,:) < 0_dp )
1080	dl_value(:,:) = dl_value(:,:)+360._dp
1081	END WHERE
1082	ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
1083	WHERE( dl_value(:,:) > 180 )
1084	dl_value(:,:) = dl_value(:,:)-180._dp
1085	END WHERE
1086	ENDIF
1087	ENDIF
1088
1089	WHERE( dl_value(:, 2) /= dd_fill .AND. & ! jj
1090	& dl_value(:, 3) /= dd_fill .AND. & ! jj+1
1091	& dl_value(:, 1) /= dd_fill ) ! jj-1
1092
1093	df_deriv(:,jj)= (dl_value(:,3) - dl_value(:,1)) / REAL(il_jmax-il_jmin,dp)
1094
1095	END WHERE
1096
1097	ENDDO
1098
1099	END SELECT
1100
1101	DEALLOCATE( dl_value )
1102
1103	END FUNCTION math_deriv_2D
1104	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1105	PURE FUNCTION math_deriv_3D(dd_value, dd_fill, cd_dim, ld_discont) &
1106	& RESULT (df_deriv)
1107	!-------------------------------------------------------------------
1108	!> @brief
1109	!> This function compute derivative of 3D array.
1110	!> you have to specify in which direction derivative have to be computed:
1111	!> first (I), second (J) or third (K) dimension.
1112	!>
1113	!> @details
1114	!> optionaly you could specify to take into account east west discontinuity
1115	!> (-180° 180° or 0° 360° for longitude variable)
1116	!>
1117	!> @author J.Paul
1118	!> @date November, 2013 - Initial Version
1119	!>
1120	!> @param[inout] dd_value 3D array of variable to be extrapolated
1121	!> @param[in] dd_fill FillValue of variable
1122	!> @param[in] cd_dim compute derivative on first (I) second (J) or third (K) dimension
1123	!> @param[in] ld_discont logical to take into account east west discontinuity
1124	!-------------------------------------------------------------------
1125
1126	IMPLICIT NONE
1127
1128	! Argument
1129	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value
1130	REAL(dp) , INTENT(IN) :: dd_fill
1131	CHARACTER(LEN=*) , INTENT(IN) :: cd_dim
1132	LOGICAL , INTENT(IN), OPTIONAL :: ld_discont
1133
1134	! function
1135	REAL(dp), DIMENSION(SIZE(dd_value,DIM=1), &
1136	& SIZE(dd_value,DIM=2), &
1137	& SIZE(dd_value,DIM=3)) :: df_deriv
1138
1139	! local variable
1140	INTEGER(i4) :: il_imin
1141	INTEGER(i4) :: il_imax
1142	INTEGER(i4) :: il_jmin
1143	INTEGER(i4) :: il_jmax
1144	INTEGER(i4) :: il_kmin
1145	INTEGER(i4) :: il_kmax
1146	INTEGER(i4), DIMENSION(3) :: il_shape
1147
1148	REAL(dp) :: dl_min
1149	REAL(dp) :: dl_max
1150	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_value
1151
1152	LOGICAL :: ll_discont
1153
1154	! loop indices
1155	INTEGER(i4) :: ji
1156	INTEGER(i4) :: jj
1157	INTEGER(i4) :: jk
1158
1159	INTEGER(i4) :: i1
1160	INTEGER(i4) :: i2
1161
1162	INTEGER(i4) :: j1
1163	INTEGER(i4) :: j2
1164
1165	INTEGER(i4) :: k1
1166	INTEGER(i4) :: k2
1167	!----------------------------------------------------------------
1168	! init
1169	df_deriv(:,:,:)=dd_fill
1170
1171	ll_discont=.FALSE.
1172	IF( PRESENT(ld_discont) ) ll_discont=ld_discont
1173
1174	il_shape(:)=SHAPE(dd_value(:,:,:))
1175
1176
1177	SELECT CASE(TRIM(fct_upper(cd_dim)))
1178
1179	CASE('I')
1180
1181	ALLOCATE( dl_value(3,il_shape(2),il_shape(3)) )
1182	! compute derivative in i-direction
1183	DO ji=1,il_shape(1)
1184
1185	il_imin=MAX(ji-1,1)
1186	il_imax=MIN(ji+1,il_shape(1))
1187
1188	IF( il_imin==ji-1 .AND. il_imax==ji+1 )THEN
1189	i1=1 ; i2=3
1190	ELSEIF( il_imin==ji .AND. il_imax==ji+1 )THEN
1191	i1=1 ; i2=2
1192	ELSEIF( il_imin==ji-1 .AND. il_imax==ji )THEN
1193	i1=2 ; i2=3
1194	ENDIF
1195
1196	dl_value(i1:i2,:,:)=dd_value(il_imin:il_imax,:,:)
1197	IF( il_imin == 1 )THEN
1198	dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), &
1199	& DIM=1, &
1200	& SHIFT=-1, &
1201	& BOUNDARY=dl_value(1,:,:) )
1202	ENDIF
1203	IF( il_imax == il_shape(1) )THEN
1204	dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), &
1205	& DIM=1, &
1206	& SHIFT=1, &
1207	& BOUNDARY=dl_value(3,:,:))
1208	ENDIF
1209
1210	IF( ll_discont )THEN
1211	dl_min=MINVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill )
1212	dl_max=MAXVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill )
1213	IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
1214	WHERE( dl_value(:,:,:) < 0_dp )
1215	dl_value(:,:,:) = dl_value(:,:,:)+360._dp
1216	END WHERE
1217	ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
1218	WHERE( dl_value(:,:,:) > 180 )
1219	dl_value(:,:,:) = dl_value(:,:,:)-180._dp
1220	END WHERE
1221	ENDIF
1222	ENDIF
1223
1224	WHERE( dl_value(2,:,:) /= dd_fill .AND. & ! ji
1225	& dl_value(3,:,:) /= dd_fill .AND. & !ji+1
1226	& dl_value(1,:,:) /= dd_fill ) !ji-1
1227
1228	df_deriv(ji,:,:)= (dl_value(3,:,:) - dl_value(1,:,:)) / REAL(il_imax-il_imin,dp)
1229
1230	END WHERE
1231
1232	ENDDO
1233
1234	CASE('J')
1235
1236	ALLOCATE( dl_value(il_shape(1),3,il_shape(3)) )
1237	! compute derivative in j-direction
1238	DO jj=1,il_shape(2)
1239
1240	il_jmin=MAX(jj-1,1)
1241	il_jmax=MIN(jj+1,il_shape(2))
1242
1243	IF( il_jmin==jj-1 .AND. il_jmax==jj+1 )THEN
1244	j1=1 ; j2=3
1245	ELSEIF( il_jmin==jj .AND. il_jmax==jj+1 )THEN
1246	j1=1 ; j2=2
1247	ELSEIF( il_jmin==jj-1 .AND. il_jmax==jj )THEN
1248	j1=2 ; j2=3
1249	ENDIF
1250
1251	dl_value(:,j1:j2,:)=dd_value(:,il_jmin:il_jmax,:)
1252	IF( il_jmin == 1 )THEN
1253	dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), &
1254	& DIM=2, &
1255	& SHIFT=-1, &
1256	& BOUNDARY=dl_value(:,1,:) )
1257	ENDIF
1258	IF( il_jmax == il_shape(2) )THEN
1259	dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), &
1260	& DIM=2, &
1261	& SHIFT=1, &
1262	& BOUNDARY=dl_value(:,3,:))
1263	ENDIF
1264
1265	IF( ll_discont )THEN
1266	dl_min=MINVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill )
1267	dl_max=MAXVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill )
1268	IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
1269	WHERE( dl_value(:,:,:) < 0_dp )
1270	dl_value(:,:,:) = dl_value(:,:,:)+360._dp
1271	END WHERE
1272	ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
1273	WHERE( dl_value(:,:,:) > 180 )
1274	dl_value(:,:,:) = dl_value(:,:,:)-180._dp
1275	END WHERE
1276	ENDIF
1277	ENDIF
1278
1279	WHERE( dl_value(:, 2,:) /= dd_fill .AND. & ! jj
1280	& dl_value(:, 3,:) /= dd_fill .AND. & ! jj+1
1281	& dl_value(:, 1,:) /= dd_fill ) ! jj-1
1282
1283	df_deriv(:,jj,:)= (dl_value(:,3,:) - dl_value(:,1,:)) / REAL(il_jmax - il_jmin,dp)
1284
1285	END WHERE
1286
1287	ENDDO
1288
1289	CASE('K')
1290
1291	ALLOCATE( dl_value(il_shape(1),il_shape(2),3) )
1292	! compute derivative in k-direction
1293	DO jk=1,il_shape(3)
1294
1295	il_kmin=MAX(jk-1,1)
1296	il_kmax=MIN(jk+1,il_shape(3))
1297
1298	IF( il_kmin==jk-1 .AND. il_kmax==jk+1 )THEN
1299	k1=1 ; k2=3
1300	ELSEIF( il_kmin==jk .AND. il_kmax==jk+1 )THEN
1301	k1=1 ; k2=2
1302	ELSEIF( il_kmin==jk-1 .AND. il_kmax==jk )THEN
1303	k1=2 ; k2=3
1304	ENDIF
1305
1306	dl_value(:,:,k1:k2)=dd_value(:,:,il_kmin:il_kmax)
1307	IF( il_kmin == 1 )THEN
1308	dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), &
1309	& DIM=3, &
1310	& SHIFT=-1, &
1311	& BOUNDARY=dl_value(:,:,1) )
1312	ENDIF
1313	IF( il_kmax == il_shape(3) )THEN
1314	dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), &
1315	& DIM=3, &
1316	& SHIFT=1, &
1317	& BOUNDARY=dl_value(:,:,3))
1318	ENDIF
1319
1320	IF( ll_discont )THEN
1321	dl_min=MINVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill )
1322	dl_max=MAXVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill )
1323	IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
1324	WHERE( dl_value(:,:,:) < 0_dp )
1325	dl_value(:,:,:) = dl_value(:,:,:)+360._dp
1326	END WHERE
1327	ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
1328	WHERE( dl_value(:,:,:) > 180 )
1329	dl_value(:,:,:) = dl_value(:,:,:)-180._dp
1330	END WHERE
1331	ENDIF
1332	ENDIF
1333
1334	WHERE( dl_value(:,:,2) /= dd_fill .AND. & ! jk
1335	& dl_value(:,:,3) /= dd_fill .AND. & ! jk+1
1336	& dl_value(:,:,1) /= dd_fill ) ! jk-1
1337
1338	df_deriv(:,:,jk)= (dl_value(:,:,3) - dl_value(:,:,1)) / REAL(il_kmax-il_kmin,dp)
1339
1340	END WHERE
1341
1342	ENDDO
1343
1344	END SELECT
1345
1346	DEALLOCATE( dl_value )
1347
1348	END FUNCTION math_deriv_3D
1349	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1350	FUNCTION math_ortho(dd_latm) &
1351	& RESULT(df_ortho)
1352	!-------------------------------------------------------------------
1353	!> @brief
1354	!> This function compute orthodome distance between opposite point of a cell
1355	!> of one degree.
1356	!>
1357	!> @details
1358	!>
1359	!> @author J.Paul
1360	!> @date April, 2017 - Initial Version
1361	!>
1362	!> @param[in] dd_latm mean latitude of the cell
1363	!> @return orthodome distance
1364	!-------------------------------------------------------------------
1365
1366	IMPLICIT NONE
1367
1368	! Argument
1369	REAL(dp), TARGET :: dd_latm
1370
1371	! function
1372	REAL(dp) :: df_ortho
1373
1374	! local
1375	REAL(dp) :: dl_dlat
1376	REAL(dp) :: dl_dlon
1377	REAL(dp) :: dl_lat1
1378	REAL(dp) :: dl_lat2
1379	REAL(dp) :: dl_tmp
1380	!----------------------------------------------------------------
1381
1382	! one degree cell
1383	dl_dlat= 1._dp * dp_deg2rad
1384	dl_dlon= 1._dp * dp_deg2rad
1385
1386	!
1387	dl_lat1 = (dd_latm - 0.5_dp) * dp_deg2rad
1388	dl_lat2 = (dd_latm + 0.5_dp) * dp_deg2rad
1389
1390	dl_tmp = SQRT( SIN(dl_dlat0.5)*2 + &
1391	& COS(dl_lat1)COS(dl_lat2)SIN(dl_dlon0.5)*2 )
1392
1393	df_ortho= 2* dp_rearth * ASIN( dl_tmp )
1394
1395	END FUNCTION math_ortho
1396	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1397	FUNCTION math_euclid(dd_lonm,dd_latm) &
1398	& RESULT(df_euclid)
1399	!-------------------------------------------------------------------
1400	!> @brief
1401	!> This function compute euclidian distance between opposite point of a cell
1402	!> of one degree, center on (lonm,latm).
1403	!>
1404	!> @details
1405	!>
1406	!> @author J.Paul
1407	!> @date April, 2017 - Initial Version
1408	!>
1409	!> @param[in] dd_lonm mean longitude of the cell
1410	!> @param[in] dd_latm mean latitude of the cell
1411	!> @return orthodome distance
1412	!-------------------------------------------------------------------
1413
1414	IMPLICIT NONE
1415
1416	! Argument
1417	REAL(dp), TARGET :: dd_lonm
1418	REAL(dp), TARGET :: dd_latm
1419
1420	! function
1421	REAL(dp) :: df_euclid
1422
1423	! local
1424	REAL(dp) :: dl_lata
1425	REAL(dp) :: dl_lona
1426	REAL(dp) :: dl_latb
1427	REAL(dp) :: dl_lonb
1428	REAL(dp) :: xa,ya,za
1429	REAL(dp) :: xb,yb,zb
1430	!----------------------------------------------------------------
1431
1432	dl_lata=(dd_latm-0.5)*dp_deg2rad
1433	dl_lona=(dd_lonm-0.5)*dp_deg2rad
1434
1435	xa = dp_rearth * COS(dl_lata) * COS(dl_lona)
1436	ya = dp_rearth * COS(dl_lata) * SIN(dl_lona)
1437	za = dp_rearth * SIN(dl_lata)
1438
1439	dl_latb=(dd_latm+0.5)*dp_deg2rad
1440	dl_lonb=(dd_lonm+0.5)*dp_deg2rad
1441
1442	xb = dp_rearth * COS(dl_latb) * COS(dl_lonb)
1443	yb = dp_rearth * COS(dl_latb) * SIN(dl_lonb)
1444	zb = dp_rearth * SIN(dl_latb)
1445
1446	df_euclid = ((xb-xa)2 + (yb-ya)2 + (zb-za)**2)
1447
1448	END FUNCTION math_euclid
1449	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1450	END MODULE math

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

Download in other formats: