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 branches/UKMO/dev_r5518_v3.6_asm_nemovar_community/NEMOGCM/TOOLS/SIREN/src – NEMO

source: branches/UKMO/dev_r5518_v3.6_asm_nemovar_community/NEMOGCM/TOOLS/SIREN/src/math.f90 @ 7731

Last change on this file since 7731 was 7731, checked in by dford, 7 years ago

Merge in revisions 6625:7726 of dev_r5518_v3.4_asm_nemovar_community, so this branch will be identical to revison 7726 of dev_r5518_v3.6_asm_nemovar_community.

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