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.
fbgenerate_coords.F90 in branches/UKMO/r6232_tracer_advection/NEMOGCM/TOOLS/OBSTOOLS/src – NEMO

source: branches/UKMO/r6232_tracer_advection/NEMOGCM/TOOLS/OBSTOOLS/src/fbgenerate_coords.F90 @ 9295

Last change on this file since 9295 was 4990, checked in by timgraham, 9 years ago

Merged branches/2014/dev_MERGE_2014 back onto the trunk as follows:

In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
1 conflict in LIM_SRC_3/limdiahsb.F90
Resolved by keeping the version from dev_MERGE_2014 branch
and commited at r4989

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2014/dev_MERGE_2014
to merge the branch into the trunk - no conflicts at this stage.

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.