trunk/dyn3d/conf_dat2d.f

module conf_dat2d_m

  ! This module is clean: no C preprocessor directive, no include line.
  ! From conf_dat2d.F, version 1.2 2006/01/27 15:14:22

  IMPLICIT NONE

contains

  SUBROUTINE conf_dat2d(xd, yd, xf, yf, champd, interbar)

    ! Auteur : P. Le Van

    ! Ce sous-programme configure le champ de données 2D 'champd' et
    ! les longitudes et latitudes de telle façon qu'on ait - pi à pi
    ! en longitude et pi/2 à - pi/2 en latitude.
    
    ! This procedure receives a 2D field, with the corresponding
    ! coordinate variables: longitude and latitude.
    ! The procedure converts longitude and latitude to radians, if the
    ! input values are in degrees.
    ! If the input longitudes are between 0 and 2 pi, the procedure
    ! computes the congruent longitudes between -pi and pi, and permutes
    ! them so they stay in increasing order.
    ! If the input latitudes are from south pole to north pole, the
    ! procedure permutes them so they become from north to south.
    ! Any change on longitudes or latitudes induces a change on the 2D field.
    ! If required, the longitudes and latitudes are finally replaced
    ! by their mid-values.

    use numer_rec, only: assert_eq
    use comconst, only: pi

    REAL, intent(in):: xd(:)
    ! (longitudes, in degrees or radians, in increasing order, from 0°
    ! to 360° or -180° to 180°)

    REAL, intent(in):: yd(:)
    ! (latitudes, in degrees or radians, in increasing or decreasing
    ! order, from pole to pole)

    LOGICAL, intent(in), optional:: interbar
    REAL, intent(out):: xf(:), yf(:) ! longitudes and latitudes, in rad
    REAL, intent(inout):: champd(:, :)

    ! Variables locales:

    INTEGER  lons, lats
    LOGICAL radianlon ! "xd" is in degrees
    logical invlon ! "xd" contains longitudes between 0 and 2 pi
    logical radianlat ! "yd" is in rad
    REAL rlatmin, rlatmax, old_xf_1
    INTEGER i, j
    logical mid_values

    !------------------------------

    lons = assert_eq(size(xd), size(xf), size(champd, 1), "conf_dat2d lons")
    lats = assert_eq(size(yd), size(yf), size(champd, 2), "conf_dat2d lats")

    IF (xd(1) >= - pi -0.5 .AND. xd(lons) <= pi + 0.5) THEN
       radianlon = .TRUE.
       invlon    = .FALSE.
    ELSE IF (xd(1) >= -0.5 .AND. xd(lons) <= 2 * pi+0.5) THEN
       radianlon = .TRUE.
       invlon    = .TRUE.
    ELSE IF (xd(1) >= -180.5 .AND. xd(lons) <= 180.5) THEN
       radianlon = .FALSE.
       invlon    = .FALSE.
    ELSE IF (xd(1) >= -0.5 .AND. xd(lons) <= 360.5) THEN
       radianlon = .FALSE.
       invlon    = .TRUE.
    ELSE
       stop '"conf_dat2d": problem with longitudes'
    ENDIF

    rlatmin = MIN(yd(1), yd(lats))
    rlatmax = MAX(yd(1), yd(lats))

    IF (rlatmin >= -pi / 2 - 0.5 .AND. rlatmax <= pi / 2 + 0.5)THEN
       radianlat = .TRUE.
    ELSE IF (rlatmin >= -90. - 0.5 .AND. rlatmax <= 90. + 0.5) THEN
       radianlat = .FALSE.
    ELSE
       stop '"conf_dat2d": problem with latitudes'
    ENDIF

    IF (radianlon)  THEN
       xf(:) = xd(:)
    else
       xf(:) = xd(:) * pi / 180. ! convert to rad
    ENDIF

    IF (radianlat)  THEN
       yf(:) = yd(:)
    else
       yf(:) = yd(:) * pi / 180. ! convert to rad
    ENDIF

    IF (invlon) THEN
       ! On tourne les longitudes pour avoir - pi à + pi :

       ! Get the index of the first longitude > pi:
       i = 1
       do while (xf(i) <= pi)
          i = i + 1
       end do

       xf(i:) = xf(i:) - 2 * pi
       xf(:) = cshift(xf, shift=i - 1)
       champd(:, :) = cshift(champd, shift=i - 1)
    ENDIF

    IF (yd(1) < yd(lats)) THEN
       ! "yd" contains latitudes from south pole to north pole,
       ! reverse their order in "yf":
       yf(lats:1:-1) = yf(:)
       champd(:, lats:1:-1) = champd(:, :)
    ENDIF

    if (present(interbar)) then
       mid_values = interbar
    else
       mid_values = .true. ! default
    end if
    if (mid_values) then
       ! Replace longitudes and latitudes by their mid-values:
       old_xf_1 = xf(1)
       forall (i = 1: lons - 1) xf(i) = 0.5 * (xf(i) + xf(i+1))
       xf(lons) = 0.5 * (xf(lons) + old_xf_1 + 2 * pi)
          
       forall (j = 1: lats - 1) yf(j) = 0.5 * (yf(j) + yf(j+1))
    end if

  END SUBROUTINE conf_dat2d

end module conf_dat2d_m
1	module conf_dat2d_m
2
3	! This module is clean: no C preprocessor directive, no include line.
4	! From conf_dat2d.F, version 1.2 2006/01/27 15:14:22
5
6	IMPLICIT NONE
7
8	contains
9
10	SUBROUTINE conf_dat2d(xd, yd, xf, yf, champd, interbar)
11
12	! Auteur : P. Le Van
13
14	! Ce sous-programme configure le champ de données 2D 'champd' et
15	! les longitudes et latitudes de telle façon qu'on ait - pi à pi
16	! en longitude et pi/2 à - pi/2 en latitude.
17
18	! This procedure receives a 2D field, with the corresponding
19	! coordinate variables: longitude and latitude.
20	! The procedure converts longitude and latitude to radians, if the
21	! input values are in degrees.
22	! If the input longitudes are between 0 and 2 pi, the procedure
23	! computes the congruent longitudes between -pi and pi, and permutes
24	! them so they stay in increasing order.
25	! If the input latitudes are from south pole to north pole, the
26	! procedure permutes them so they become from north to south.
27	! Any change on longitudes or latitudes induces a change on the 2D field.
28	! If required, the longitudes and latitudes are finally replaced
29	! by their mid-values.
30
31	use numer_rec, only: assert_eq
32	use comconst, only: pi
33
34	REAL, intent(in):: xd(:)
35	! (longitudes, in degrees or radians, in increasing order, from 0°
36	! to 360° or -180° to 180°)
37
38	REAL, intent(in):: yd(:)
39	! (latitudes, in degrees or radians, in increasing or decreasing
40	! order, from pole to pole)
41
42	LOGICAL, intent(in), optional:: interbar
43	REAL, intent(out):: xf(:), yf(:) ! longitudes and latitudes, in rad
44	REAL, intent(inout):: champd(:, :)
45
46	! Variables locales:
47
48	INTEGER lons, lats
49	LOGICAL radianlon ! "xd" is in degrees
50	logical invlon ! "xd" contains longitudes between 0 and 2 pi
51	logical radianlat ! "yd" is in rad
52	REAL rlatmin, rlatmax, old_xf_1
53	INTEGER i, j
54	logical mid_values
55
56	!------------------------------
57
58	lons = assert_eq(size(xd), size(xf), size(champd, 1), "conf_dat2d lons")
59	lats = assert_eq(size(yd), size(yf), size(champd, 2), "conf_dat2d lats")
60
61	IF (xd(1) >= - pi -0.5 .AND. xd(lons) <= pi + 0.5) THEN
62	radianlon = .TRUE.
63	invlon = .FALSE.
64	ELSE IF (xd(1) >= -0.5 .AND. xd(lons) <= 2 * pi+0.5) THEN
65	radianlon = .TRUE.
66	invlon = .TRUE.
67	ELSE IF (xd(1) >= -180.5 .AND. xd(lons) <= 180.5) THEN
68	radianlon = .FALSE.
69	invlon = .FALSE.
70	ELSE IF (xd(1) >= -0.5 .AND. xd(lons) <= 360.5) THEN
71	radianlon = .FALSE.
72	invlon = .TRUE.
73	ELSE
74	stop '"conf_dat2d": problem with longitudes'
75	ENDIF
76
77	rlatmin = MIN(yd(1), yd(lats))
78	rlatmax = MAX(yd(1), yd(lats))
79
80	IF (rlatmin >= -pi / 2 - 0.5 .AND. rlatmax <= pi / 2 + 0.5)THEN
81	radianlat = .TRUE.
82	ELSE IF (rlatmin >= -90. - 0.5 .AND. rlatmax <= 90. + 0.5) THEN
83	radianlat = .FALSE.
84	ELSE
85	stop '"conf_dat2d": problem with latitudes'
86	ENDIF
87
88	IF (radianlon) THEN
89	xf(:) = xd(:)
90	else
91	xf(:) = xd(:) * pi / 180. ! convert to rad
92	ENDIF
93
94	IF (radianlat) THEN
95	yf(:) = yd(:)
96	else
97	yf(:) = yd(:) * pi / 180. ! convert to rad
98	ENDIF
99
100	IF (invlon) THEN
101	! On tourne les longitudes pour avoir - pi à + pi :
102
103	! Get the index of the first longitude > pi:
104	i = 1
105	do while (xf(i) <= pi)
106	i = i + 1
107	end do
108
109	xf(i:) = xf(i:) - 2 * pi
110	xf(:) = cshift(xf, shift=i - 1)
111	champd(:, :) = cshift(champd, shift=i - 1)
112	ENDIF
113
114	IF (yd(1) < yd(lats)) THEN
115	! "yd" contains latitudes from south pole to north pole,
116	! reverse their order in "yf":
117	yf(lats:1:-1) = yf(:)
118	champd(:, lats:1:-1) = champd(:, :)
119	ENDIF
120
121	if (present(interbar)) then
122	mid_values = interbar
123	else
124	mid_values = .true. ! default
125	end if
126	if (mid_values) then
127	! Replace longitudes and latitudes by their mid-values:
128	old_xf_1 = xf(1)
129	forall (i = 1: lons - 1) xf(i) = 0.5 * (xf(i) + xf(i+1))
130	xf(lons) = 0.5 * (xf(lons) + old_xf_1 + 2 * pi)
131
132	forall (j = 1: lats - 1) yf(j) = 0.5 * (yf(j) + yf(j+1))
133	end if
134
135	END SUBROUTINE conf_dat2d
136
137	end module conf_dat2d_m