source: CONFIG_DEVT/IPSLCM6.5_work_ENSEMBLES/modeles/NEMO/tools/OBSTOOLS/src/fbgenerate_coords.F90 @ 5501

Last change on this file since 5501 was 5501, checked in by aclsce, 4 years ago

First import of IPSLCM6.5_work_ENSEMBLES working configuration

File size: 15.0 KB
Line 
1MODULE fbgenerate_coords
2USE obs_fbm
3USE date_utils
4
5CONTAINS
6
7   REAL(KIND=fbdp) FUNCTION fb_dates(timein,datein) 
8      IMPLICIT NONE   
9      INTEGER, INTENT(IN) :: timein     ! Format: HHMM
10                INTEGER, INTENT(IN) :: datein   ! Format: YYYYMMDD
11      INTEGER :: iyea,imon,iday
12      INTEGER :: ihr,imin,isec
13 
14      iyea=datein/10000
15      imon=datein/100-iyea*100
16      iday=datein-iyea*10000-imon*100
17               
18                ihr = timein/100
19                imin = timein - ihr*100
20                isec = 0
21               
22      CALL greg2jul(isec,imin,ihr,iday,imon,iyea,fb_dates)
23   
24   END FUNCTION fb_dates
25
26
27
28        SUBROUTINE set_spatial_coords(lats,lons,n,FillVal)
29        IMPLICIT NONE   
30        INTEGER :: i, j, k, p, nlats, nlons, nlats_in_list, nlons_in_list
31        INTEGER, INTENT(IN) :: n
32        REAL(KIND=fbdp), INTENT(INOUT) :: lats(:), lons(:)
33        REAL(KIND=fbdp) :: FillVal 
34       
35        ! A single non-FillVal value should be replicated n times
36        ! Three non-FillVal values should be used as the start, end, step in conjunction with n
37        !       unless n is three, in which case treat as list, not bounds.
38        ! A full list of n non-FillVals should be left unaltered.
39       
40        !Find number of lats and lons
41        ! by finding number of non-FillVal entries in arrays
42        ! and expanding bounds if necessary
43       
44        nlats_in_list = last_element(lats,FillVal)
45        nlons_in_list = last_element(lons,FillVal)
46
47   !Expand bounds only if 3 values given and number obs>3
48        IF (nlats_in_list == 1) THEN
49                lats(1:nlats) = lats(1)
50        ELSE IF ((nlats_in_list==3).AND.((n>3))) THEN ! Treat as list of three, not bounds
51                CALL expand_bounds(lats, FillVal)
52        END IF
53
54        IF (nlons_in_list == 1) THEN
55                lons(1:nlons) = lons(1)
56        ELSE IF ((nlons_in_list==3).AND.((n>3))) THEN ! Treat as list of three, not bounds
57                CALL expand_bounds(lons, FillVal)
58        END IF
59
60        nlats = last_element(lats,FillVal)
61        nlons = last_element(lons,FillVal)
62       
63        IF ((n /= nlats) .AND. (n /= nlons)) THEN
64                WRITE(*,*) "ERROR: Number of lat/lons not equal to nobs", nlats, nlons, n
65                CALL abort     
66        END IF
67       
68        END SUBROUTINE set_spatial_coords
69
70
71        SUBROUTINE set_spatial_coords_grid(lats,lons,n,nlats,nlons,FillVal)
72        IMPLICIT NONE   
73        INTEGER :: i, j, k, p, nlats_in_list, nlons_in_list
74        INTEGER, INTENT(IN) :: n, nlats, nlons
75        REAL(KIND=fbdp), INTENT(INOUT) :: lats(:), lons(:)
76        REAL(KIND=fbdp), ALLOCATABLE :: tmp_lats(:), tmp_lons(:)
77        REAL(KIND=fbdp) :: FillVal 
78       
79        ! A single non-FillVal value should be replicated n times
80        ! Three non-FillVal values should be used as the start, end, step in conjunction with n
81        !       unless n is three, in which case treat as list, not bounds.
82        ! A full list of n non-FillVals should be left unaltered.
83       
84        !Find number of lats and lons in list
85        ! by finding number of non-FillVal entries in arrays
86        ! and then expand bounds if necessary
87       
88        nlats_in_list = last_element(lats,FillVal)
89        nlons_in_list = last_element(lons,FillVal)
90       
91   !Expand bounds only if 3 values given and number lats (or lons) not equal to 3
92        IF (nlats_in_list == 1) THEN
93                lats(1:nlats) = lats(1)
94        ELSE IF ((nlats_in_list==3).AND.((nlats>3))) THEN ! Treat as bounds
95                CALL expand_bounds(lats, FillVal)
96        END IF
97
98        IF (nlons_in_list == 1) THEN
99                lons(1:nlons) = lons(1)
100        ELSE IF ((nlons_in_list==3).AND.((nlons>3))) THEN ! Treat as bounds
101                CALL expand_bounds(lons, FillVal)
102        END IF
103
104        nlats_in_list = last_element(lats,FillVal)
105        nlons_in_list = last_element(lons,FillVal)
106
107
108        IF ((nlats_in_list * nlons_in_list) == n) THEN
109                ALLOCATE(tmp_lons(n))
110                ALLOCATE(tmp_lats(n))
111                k = 0
112                p = 0
113                DO i = 1, nlats
114                   p = p + 1
115                        DO j = 1, nlons
116                                k = k + 1
117                                tmp_lons(k) = lons(j)   
118                                tmp_lats(k) = lats(p)   
119                        END DO
120                END DO
121                lats(:) = tmp_lats(:)
122                lons(:) = tmp_lons(:)
123                DEALLOCATE(tmp_lons)
124                DEALLOCATE(tmp_lats)
125        ELSE
126                WRITE(*,*) 
127                WRITE(*,*) "ERROR: Number of lat * lon values not equal to nobs", nlats_in_list, nlons_in_list, n
128                CALL abort     
129        END IF
130       
131        END SUBROUTINE set_spatial_coords_grid
132
133
134        SUBROUTINE set_spatial_coords_physpace(lats,lons,n,FillVal,phys_spacing)
135        IMPLICIT NONE   
136        INTEGER :: i, j, k, nlats, nlons
137        INTEGER, INTENT(IN) :: n
138        REAL(KIND=fbdp), INTENT(INOUT) :: lats(:), lons(:)
139        REAL(KIND=fbdp), ALLOCATABLE :: tmp_lats(:)
140        REAL(KIND=fbdp) :: FillVal 
141   REAL(KIND=fbdp), PARAMETER :: earth_radius = 6371.0_fbdp ! km
142   REAL(KIND=fbdp), PARAMETER :: pi = 3.141592654_fbdp
143        REAL(KIND=fbdp), INTENT(IN) :: phys_spacing
144   REAL(KIND=fbdp) :: lat_step
145       
146        ! Use physical spacing to set lats and lons
147
148   nlats = INT(pi * earth_radius / phys_spacing)
149   IF (MOD(nlats,2)==0) nlats=nlats-1   ! Ensure nlats is odd to have even number around equator. 
150     
151   ALLOCATE(tmp_lats(nlats))
152
153   lat_step = 180.0_fbdp / (nlats-1)
154
155   tmp_lats(1) = 0.0_fbdp
156   DO i=1,(nlats-1)/2
157      tmp_lats(i+1) = tmp_lats(i) + lat_step
158      tmp_lats(((nlats-1)/2)+i+1) = tmp_lats(i+1) * (-1.0_fbdp)
159   END DO
160
161   k = 0
162   DO i=1,nlats
163      nlons = INT(2.0_fbdp * pi * earth_radius * cos(tmp_lats(i)* (pi/180.0_fbdp)) / phys_spacing)
164      IF (nlons>0) THEN
165         DO j = 1, nlons
166            k = k + 1
167            lons(k) = REAL(j-1) * 360.0_fbdp / nlons
168            lats(k) = tmp_lats(i)
169         END DO
170      END IF
171   END DO
172   DEALLOCATE(tmp_lats)
173       
174        END SUBROUTINE set_spatial_coords_physpace
175
176
177        SUBROUTINE expand_bounds(array, fillval)
178        !
179        ! Checks if there are three entries and expands bounds.
180        !
181        IMPLICIT NONE
182        INTEGER :: i, nentries
183        REAL(KIND=fbdp), INTENT(INOUT) :: array(:)
184        REAL(KIND=fbdp), INTENT(IN) :: fillval
185        REAL(KIND=fbdp) :: nsteps, start, step
186
187        ! Find number of elements given in list
188        nentries = last_element(array, fillval)
189
190        IF (nentries==3) THEN                   ! Bounds given 
191                nsteps = (array(2) - array(1)) /  array(3) 
192           IF ((array(2) < array(1)) .AND. (array(3) >= 0.0)) THEN
193                   WRITE(*,*) "Error: if upper bound is less than lower bound, step must be negative.", array
194                        CALL abort
195                ELSE IF ( NINT(nsteps)+1 > SIZE(array) ) THEN
196                   WRITE(*,*) "Error: bounds not compatible with length of array.", array
197                   WRITE(*,*) "Error: bounds not compatible with length of array.", nsteps
198                        CALL abort
199                END IF
200                start = array(1) 
201                step = array(3)
202                DO i=1,NINT(nsteps)+1
203                        array(i) = start + (i-1)*step           
204                END DO
205        END IF
206       
207        END SUBROUTINE expand_bounds
208
209
210
211        INTEGER FUNCTION last_element(array,fillval)
212        !
213        ! Returns index of the last non-fill value element in a list
214        ! Will return 1, even in first element is the FillValue
215        !
216        IMPLICIT NONE
217        REAL(KIND=fbdp), INTENT(IN) :: array(:), fillval
218        INTEGER :: i
219
220        last_element = 1
221        IF (SIZE(array) > 1) THEN
222                DO i=2, SIZE(array)
223                        IF (array(i) /= fillval) THEN
224                                last_element = i
225                        END IF
226                END DO
227        END IF
228               
229        END FUNCTION last_element
230
231
232
233        SUBROUTINE set_depths(array,n,FillVal)
234        IMPLICIT NONE   
235        INTEGER :: i
236        INTEGER, INTENT(IN) :: n
237        REAL(KIND=fbdp), INTENT(INOUT) :: array(:)
238        REAL(KIND=fbdp) :: FillVal
239        REAL(KIND=fbdp) :: start, step, nstep
240        ! Top bound not necessarily inclusive - will use start, step and number of obs
241       
242        ! A single non-FillVal value should be replicated n times
243        ! Three non-FillVal values should be used as the start, end, step in conjunction with n
244        !       unless n is three, in which case treat as list, not bounds.
245        ! A full list of n non-FillVals should be left unaltered.
246        IF (n==1) THEN
247                array(:) = array(1)
248        ELSE IF ((n > 1) .AND. (array(2)==FillVal)) THEN
249                array(:) = array(1)
250        ELSE IF ((n==2) .AND. (array(2)/=FillVal)) THEN
251                array(:) = array(:)
252        ELSE IF (n==3) THEN ! Treat as list of three, not bounds
253                array(:) = array(:)
254        ELSE IF ((n>=4) .AND. (array(4)==FillVal)) THEN ! Assume start, stop, step
255                nstep = (array(2) - array(1)) /  array(3) 
256           IF ((array(2) < array(1)) .AND. (array(3) >= 0.0)) THEN
257                   WRITE(*,*) "Error: if upper bound is less than lower bound, step must be negative.", array
258                        CALL abort
259                ELSE IF ( ( nstep < (0.99 * real(n-1)) ) .OR. &
260                          ( nstep > (1.01 * real(n-1)) ) ) THEN
261                   WRITE(*,*) "Error: depth bounds not compatible.", array
262                   WRITE(*,*) "Error: depth bounds not compatible.", nstep
263                        CALL abort
264                END IF
265                start = array(1) 
266                step = array(3)
267                DO i=1,n
268                        array(i) = start + (i-1)*step           
269                END DO
270        END IF
271       
272        END SUBROUTINE set_depths
273
274
275
276
277
278        SUBROUTINE set_datetime(date,time,julian_date,n,FillVal)
279        !
280        ! Transform input array from m values describing n dates to a list of n dates.
281        !
282        IMPLICIT NONE   
283        INTEGER :: i
284        INTEGER, INTENT(IN) :: n
285        INTEGER, INTENT(INOUT) :: date(:), time(:)
286        INTEGER :: FillVal
287        REAL(KIND=fbdp), INTENT(INOUT) :: julian_date(:)
288        REAL(KIND=fbdp) :: start, step, finish
289
290   ! Check if user has supplied start and end times to be split amongst number of obs
291   ! (i.e. more than 2 obs, but only 2 dates and 2 times given.)
292   IF ((n>2) .AND. (date(2)/=FillVal) .AND. (date(3)==FillVal) &
293      &      .AND. (time(2)/=FillVal) .AND. (time(3)==FillVal)) THEN
294     
295      !work out JD of each
296      start = fb_dates(time(1),date(1)) 
297      finish = fb_dates(time(2),date(2)) 
298
299      ! use info to calc step
300      step = ( finish - start ) / n
301      ! calc all JDs and put output in dates
302      DO i=1,n 
303              ! Could make this an elemental function if dateutils were pure funcs
304         julian_date(i)= start + (i-1)*step
305      END DO
306
307   ELSE
308      CALL set_date(date,n,FillVal)
309      CALL set_time(time,n,FillVal)
310      ! Replace dates with JD including time info
311      DO i=1,n 
312              ! Could make this an elemental function if dateutils were pure funcs
313         julian_date(i)= fb_dates(time(i),date(i))
314      END DO
315   END IF
316
317
318        END SUBROUTINE set_datetime
319
320
321        SUBROUTINE set_date(array,n,FillVal)
322        !
323        ! Transform input array from m values describing n dates to a list of n dates.
324        !
325        IMPLICIT NONE   
326        INTEGER :: i
327        INTEGER, INTENT(IN) :: n
328        INTEGER, INTENT(INOUT) :: array(:)
329        INTEGER :: FillVal
330        INTEGER :: start, step, diff
331        ! Top bound not necessarily inclusive - will use start, step and number of obs
332       
333        ! A single non-FillVal value should be replicated n times
334        ! Three non-FillVal values should be used as the start, end, step in conjunction with n
335        !       unless n is three, in which case treat as list, not bounds.
336        ! A full list of n non-FillVals should be left unaltered.
337        IF (n==1) THEN
338                array(:) = array(1)
339        ELSE IF ((n > 1) .AND. (array(2)==FillVal)) THEN
340                array(:) = array(1)
341        ELSE IF ((n==2) .AND. (array(2)/=FillVal)) THEN
342                array(:) = array(:)
343        ELSE IF (n==3) THEN ! Treat as list of three, not bounds
344                array(:) = array(:)
345        ELSE IF ((n>=4) .AND. (array(4)==FillVal)) THEN ! Assume start, stop, step
346        diff = diffdate(array(2),array(1))      ! in number of days
347           IF ((array(2) < array(1)) .AND. (array(3) >= 0)) THEN
348                   WRITE(*,*) "Error: if upper bound is less than lower bound, step must be negative.", array
349                        CALL abort
350                ELSE IF ( ( diff / ABS(array(3)) ) /= (n-1) ) THEN
351                   WRITE(*,*) "Error: date bounds not compatible.", array
352                        CALL abort
353                END IF
354                start = array(1) 
355                step = array(3)
356                DO i=1,n
357                        CALL add_days_to_date(start,(i-1)*step,array(i))               
358                END DO
359        END IF
360       
361        END SUBROUTINE set_date
362
363
364
365        SUBROUTINE set_time(array,n,FillVal)
366        !
367        ! Transform input array from m (m<=n) values describing n times to a list of n times.
368        !
369        IMPLICIT NONE   
370        INTEGER :: i
371        INTEGER, INTENT(IN) :: n
372        INTEGER, INTENT(INOUT) :: array(:)
373        INTEGER :: FillVal
374        INTEGER :: start, step, nstep
375        ! Top bound not necessarily inclusive - will use start, step and number of obs
376       
377        ! A single non-FillVal value should be replicated n times
378        ! Three non-FillVal values should be used as the start, end, step in conjunction with n
379        !       unless n is three, in which case treat as list, not bounds.
380        ! A full list of n non-FillVals should be left unaltered.
381        IF (n==1) THEN
382                array(:) = array(1)
383        ELSE IF ((n > 1) .AND. (array(2)==FillVal)) THEN
384                array(:) = array(1)
385        ELSE IF ((n==2) .AND. (array(2)/=FillVal)) THEN
386                array(:) = array(:)
387        ELSE IF (n==3) THEN ! Treat as list of three, not bounds
388                array(:) = array(:)
389        ELSE IF ((n>=4) .AND. (array(4)==FillVal)) THEN ! Assume start, stop, step
390                IF (array(3)<0) nstep = difftime(array(2),array(1)) / ABS(array(3)) 
391                IF (array(3)>=0) nstep = difftime(array(1),array(2)) / ABS(array(3)) 
392                IF ( nstep /= (n-1) ) THEN
393                   WRITE(*,*) "Error: time bounds not compatible.", array
394                        CALL abort
395                END IF
396                start = array(1) 
397                step = array(3)
398                DO i=1,n
399                        array(i) = add_mins_to_time(start,(i-1)*step)           
400                END DO
401        END IF
402       
403        END SUBROUTINE set_time
404       
405       
406
407        SUBROUTINE set_obs_values(array,p,n,m,FillVal)
408        !
409        ! n profiles at m depths
410        !
411        IMPLICIT NONE   
412        INTEGER :: i, j, k
413        INTEGER, INTENT(IN) :: n, m, p
414        REAL(KIND=fbdp), INTENT(INOUT) :: array(:,:,:)                  ! (nlevels, nprofiles, nvars) = (m,n,p)
415        REAL(KIND=fbdp) :: FillVal 
416       
417       
418        DO k=1,p
419           IF ((k > 1) .AND. (array(1,1,k) == FillVal)) THEN  ! set to same values as first variable
420                        array(:,:,k) = array(:,:,1)
421                ELSE
422                        IF ((n==1).AND.(m==1)) THEN
423                                array(:,:,k) = array(1,1,k)
424
425                        ! If mult depths, but not specified, set all to one value       
426                        ELSE IF ((m > 1) .AND. (array(2,1,k) == FillVal)) THEN
427                                array(:,:,k) = array(1,1,k)
428
429                        ! If mult profiles and mult depths and only one value set in first profile     
430                        ELSE IF ((n > 1) .AND. (m > 1) .AND. (array(2,1,k) == FillVal)) THEN
431                                array(:,:,k) = array(1,1,k)
432                               
433                        ! If mult profiles and mult depths and only one first profile set       
434                        ELSE IF ((n > 1) .AND. (m > 1) .AND. (array(1,2,k) == FillVal)) THEN
435                                DO j=1,m
436                                        array(j,:,k) = array(j,1,k)
437                                END DO
438                               
439                        ELSE
440                        array(:,:,k) = array(:,:,k)
441                        END IF
442                END IF
443        END DO 
444       
445        END SUBROUTINE set_obs_values
446
447
448   ! Unbiased shuffle of array
449   SUBROUTINE shuffle(a)
450   REAL(KIND=fbdp), INTENT(INOUT) :: a(:)
451   INTEGER :: i, randpos
452   REAL(KIND=fbdp) :: r, temp
453
454   CALL random_seed()
455   DO i = SIZE(a), 2, -1
456      CALL random_number(r)
457      randpos = int(r * i) + 1
458      temp = a(randpos)
459      a(randpos) = a(i)
460      a(i) = temp
461   END DO
462   
463   END SUBROUTINE shuffle
464   
465   
466   ! Add a random perturbation to the lats and lons
467   ! Perturbation is sampled from a uniform distribution +/-perturb_limit
468   SUBROUTINE perturb_positions(lats,lons,perturb_limit)
469   INTEGER :: i
470   REAL(KIND=fbdp), INTENT(INOUT) :: lats(:), lons(:)
471   REAL(KIND=fbdp), INTENT(IN) :: perturb_limit
472   REAL(KIND=fbdp) :: randpos, lat_perturb, lon_perturb
473   REAL(KIND=fbdp), PARAMETER :: earth_radius = 6371.0_fbdp
474   REAL(KIND=fbdp), PARAMETER :: pi = 3.141592654_fbdp
475   
476        IF ( SIZE(lats) /= SIZE(lons) ) THEN
477                WRITE(*,*) "Error: different number of lat and lon elements", SIZE(lats), SIZE(lons)
478                CALL abort
479        END IF
480
481   CALL random_seed()
482   
483   ! Convert physical sep into a latidue sep in degrees
484   lat_perturb = 360.0_fbdp * perturb_limit / (2.0_fbdp * pi * earth_radius)
485   
486   DO i=1,SIZE(lats)
487   
488      ! Perturb lats first, as lon conversion uses lat
489      CALL random_number(randpos)
490      lats(i) = lats(i) + randpos*lat_perturb
491      IF (lats(i) > 90.0_fbdp) lats(i) = 180.0_fbdp - lats(i) 
492      IF (lats(i) < -90.0_fbdp) lats(i) = (lats(i) + 180.0_fbdp) * (-1.0_fbdp)
493     
494      ! Use lat to convert physical size to delta_longitude   
495      IF (ABS(lats(i)) == 90.0_fbdp) THEN
496         lon_perturb = 0.0_fbdp
497      ELSE
498         lon_perturb = 360.0_fbdp * perturb_limit / (2.0_fbdp * pi * earth_radius * cos(lats(i)*(pi/180.0_fbdp)))
499      END IF
500      CALL random_number(randpos)
501      lons(i) = lons(i) + randpos*lon_perturb
502      IF (lons(i) >= 360.0_fbdp) lons(i) = lons(i) - 360.0_fbdp 
503      IF (lons(i) < 0.0_fbdp) lons(i) = lons(i) + 360.0_fbdp 
504   
505   END DO
506   
507   END SUBROUTINE perturb_positions
508   
509END MODULE fbgenerate_coords
Note: See TracBrowser for help on using the repository browser.