source: trunk/MESOSCALE_DEV/SRC/ARWpost/src/module_interp.f90 @ 1068

Last change on this file since 1068 was 207, checked in by aslmd, 14 years ago

MESOSCALE: A GENERAL CLEAN-UP FOLLOWING UPDATING THE USER MANUAL. EVERYTHING ESSENTIAL IS IN MESOSCALE (much lighter than before). EVERYTHING FOR DEVELOPPERS OR EXPERTS IS IN MESOSCALE_DEV.
File size: 4.6 KB
RevLine 
[207]1MODULE module_interp
2
3  CONTAINS
4!--------------------------------------------------------
5
6  SUBROUTINE interp( data_in, nx, ny, nz, &
7                     data_out, nxout, nyout, nzout, &
8                     z_data, z_levs, number_of_zlevs )
9
10  USE module_model_basics
11
12  implicit none
13
14 ! Arguments
15  integer                                                       :: nx, ny, nz, number_of_zlevs
16  real, dimension(west_east_dim,south_north_dim,bottom_top_dim) :: z_data
17  real, dimension(nx,ny,nz)                                     :: data_in
18  real, pointer, dimension(:,:,:)                               :: data_out
19  real, dimension(number_of_zlevs)                              :: z_levs
20
21  ! Local variables
22  integer                                                       :: nxout, nyout, nzout
23  real, allocatable, dimension(:,:,:)                           :: SCR2
24  real, dimension(bottom_top_dim)                               :: data_in_1d, z_data_1d
25  real, dimension(number_of_zlevs)                              :: data_out_1d
26  integer                                                       :: i,j,k
27
28
29    IF ( ALLOCATED(SCR2) ) DEALLOCATE(SCR2)
30    IF ( ASSOCIATED(data_out) ) DEALLOCATE(data_out)
31    nxout = nx
32    nyout = ny
33    nzout = nz
34
35    !! We may be dealing with a staggered field
36    IF ( nx .gt. west_east_dim ) THEN
37      ALLOCATE(SCR2(west_east_dim,south_north_dim,bottom_top_dim))
38      SCR2 = 0.5*(data_in(1:west_east_dim,:,:)+data_in(2:west_east_dim+1,:,:))
39      nxout = west_east_dim
40    ELSE IF ( ny .gt. south_north_dim ) THEN
41      ALLOCATE(SCR2(west_east_dim,south_north_dim,bottom_top_dim))
42      SCR2 = 0.5*(data_in(:,1:south_north_dim,:)+data_in(:,2:south_north_dim+1,:))
43      nyout = south_north_dim
44    ELSE IF ( nz .gt. bottom_top_dim ) THEN 
45      ALLOCATE(SCR2(west_east_dim,south_north_dim,bottom_top_dim))
46      SCR2 = 0.5*(data_in(:,:,1:bottom_top_dim)+data_in(:,:,2:bottom_top_dim+1))
47      nzout = bottom_top_dim
48    ELSE
49      ALLOCATE(SCR2(nx,ny,nz))
50      SCR2 = data_in
51    ENDIF
52
53
54    IF ( iprogram .ge. 6 .AND. nzout .gt. 10 .AND.                  &
55         (vertical_type == 'p' .or. vertical_type == 'z') ) THEN
56
57      ALLOCATE(data_out(west_east_dim,south_north_dim,number_of_zlevs))
58      DO i=1,west_east_dim
59      DO j=1,south_north_dim
60 
61        DO k=1,bottom_top_dim
62          data_in_1d(k) = SCR2(i,j,k)
63          z_data_1d(k)  = z_data(i,j,k)
64        ENDDO
65 
66        CALL interp_1d( data_in_1d, z_data_1d, bottom_top_dim, &
67                        data_out_1d, z_levs, number_of_zlevs,  &
68                        vertical_type)
69 
70        DO k=1,number_of_zlevs
71          data_out(i,j,k) = data_out_1d(k)
72        ENDDO
73
74      ENDDO
75      ENDDO
76     
77      nzout = number_of_zlevs
78
79    ELSE
80
81      ALLOCATE(data_out(nxout,nyout,nzout))
82      data_out = SCR2
83      DEALLOCATE(SCR2)
84
85    ENDIF
86
87  END SUBROUTINE interp
88
89!----------------------------------------------
90
91  SUBROUTINE interp_1d( a, xa, na, b, xb, nb, vertical_type)
92
93  implicit none
94
95 ! Arguments
96  integer, intent(in)              :: na, nb
97  real, intent(in), dimension(na)  :: a, xa
98  real, intent(in), dimension(nb)  :: xb
99  real, intent(out), dimension(nb) :: b
100  character (len=1)                :: vertical_type
101
102  ! Local variables
103  real                             :: missing_value
104  integer                          :: n_in, n_out
105  real                             :: w1, w2
106  logical                          :: interp
107
108!  parameter (MISSING_VALUE=1.0E37)
109!!!!
110!!!!AYM AYM
111!!!!
112  parameter (MISSING_VALUE=-9999.)
113
114
115
116  IF ( vertical_type == 'p' ) THEN
117
118    DO n_out = 1, nb
119
120      b(n_out) = missing_value
121      interp = .false.
122      n_in = 1
123
124      DO WHILE ( (.not.interp) .and. (n_in < na) )
125        IF( (xa(n_in)   >= xb(n_out)) .and. &
126            (xa(n_in+1) <= xb(n_out))        ) THEN
127          interp = .true.
128          w1 = (xa(n_in+1)-xb(n_out))/(xa(n_in+1)-xa(n_in))
129          w2 = 1. - w1
130          b(n_out) = w1*a(n_in) + w2*a(n_in+1)
131        END IF
132        n_in = n_in +1
133      ENDDO
134
135    ENDDO
136
137  ELSE
138
139    DO n_out = 1, nb
140 
141      b(n_out) = missing_value
142      interp = .false.
143      n_in = 1
144 
145      DO WHILE ( (.not.interp) .and. (n_in < na) )
146        IF( (xa(n_in)   <= xb(n_out)) .and. &
147            (xa(n_in+1) >= xb(n_out))        ) THEN
148          interp = .true.
149          w1 = (xa(n_in+1)-xb(n_out))/(xa(n_in+1)-xa(n_in))
150          w2 = 1. - w1
151          b(n_out) = w1*a(n_in) + w2*a(n_in+1)
152        END IF
153        n_in = n_in +1
154      ENDDO
155
156    ENDDO
157
158  END IF
159
160  END SUBROUTINE interp_1d
161
162!-------------------------------------------------------------------------
163
164END MODULE module_interp
Note: See TracBrowser for help on using the repository browser.