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

Last change on this file since 12080 was 12080, checked in by jpaul, 10 months 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!----------------------------------------------------------------------
94MODULE 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
141CONTAINS
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_dlat*0.5)**2 + &
1391         &           COS(dl_lat1)*COS(dl_lat2)*SIN(dl_dlon*0.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   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1450END MODULE math
Note: See TracBrowser for help on using the repository browser.