-
Notifications
You must be signed in to change notification settings - Fork 0
/
mod_floats.F
executable file
·2781 lines (2753 loc) · 89.9 KB
/
mod_floats.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#if defined(ROW_LAND)
#define SEA_P .true.
#define SEA_U .true.
#define SEA_V .true.
#elif defined(ROW_ALLSEA)
#define SEA_P allip(j).or.ip(i,j).ne.0
#define SEA_U alliu(j).or.iu(i,j).ne.0
#define SEA_V alliv(j).or.iv(i,j).ne.0
#else
#define SEA_P ip(i,j).ne.0
#define SEA_U iu(i,j).ne.0
#define SEA_V iv(i,j).ne.0
#endif
module mod_floats
use mod_xc ! HYCOM communication interface
use mod_pipe ! HYCOM debugging interface
#if defined(STOKES)
use mod_stokes ! Stokes Drift Velocity Module
#endif
c
implicit none
c
c --- HYCOM synthetic floats, drifters and moorings
c --- See subroutine floats (below) for more information
c
integer, parameter, public :: nfldim=400 !maximum number of synthetic floats
real, allocatable, dimension(:,:,:), save, public ::
& wveli, ! interface vertical velocity
& uold2,
& vold2,
& wold2u,
& wold2d
real, allocatable, dimension(:,:), save, public ::
& dlondx,
& dlondy,
& dlatdx,
& dlatdy
real, allocatable, dimension(:,:,:), save, private ::
& wvelup,
& wveldn
real, allocatable, dimension(:,:), save, private ::
& dpdxup,
& dpdyup,
& dpdxdn,
& dpdydn
real, save :: flt(nfldim,13),
& deltfl,fldepm,tbvar,tdecri,dtturb,uturb0
integer, save :: kfloat(nfldim),iflnum(nfldim),ifltyp(nfldim),
& nflsam,nfldta,fltflg,nfladv,intpfl,iturbv,ismpfl,
& nflout,nfl,nflt,nstepfl
logical, save :: synflt,turbvel,samplfl,wvelfl,hadvfl,nonlatlon
character(48), save :: flnmflti,flnmflto,flnmfltio
contains
subroutine floats_init(m,n,time0)
use mod_cb_arrays ! HYCOM saved arrays
implicit none
integer m,n
real time0
c
c --- read initial float data
c --- initialize parameters and arrays required for floats
c
integer i,j,k,l,margin
integer nbcday,nsmday,ityp
c
include 'stmt_fns.h'
c
c --- initialize flags
c
c --- wvelfl: calculate vertical velocity?
wvelfl=.false.
c
c --- hadvfl: horizontally advect the float?
hadvfl=.false.
c
c --- nonlatlon: does the model grid cross latitude/longitude lines anywhere
c --- in the domain?
nonlatlon=.false.
c
c --- initialize file names
flnmflti = 'floats.input'
flnmflto = 'floats_out'
flnmfltio= 'floats.input_out'
c
c -------------------------
c --- set timing parameters
c -------------------------
c
c --- must have an integer number of baroclinic time steps per day
nbcday=nint(86400.0/baclin)
if (real(nbcday).ne.86400.0/baclin) then
if (mnproc.eq.1) then
write(lp,'(/ a /)')
& 'error - need integer no. of bacl. time steps/day'
call flush(lp)
endif !1st tile
call xcstop('(floats_init)')
stop '(floats_init)'
endif
c
c --- parameter nfladv is entered in blkdat.input as the number of
c --- baroclinic time steps between updates of float position. at model
c --- time nstep when the float is to be advected, velocity fields saved
c --- at times nstep, nstep-nfladv/2, and nstep-nfladv are used to perform
c --- the runga-kutta time interpolation.
c
c --- nfladv must be set so the float is advected every 2-4 hours, but
c --- also must be no smaller than 4. This ensures that the runga-
c --- kutta time interpolation is performed with a delta-time of 1-2 hours.
c
c --- if all floats are stationary (synthetic moorings), then nfladv
c --- must still be no smaller than 4, but the above time restrection
c --- is not necessary. this allows the user to sample at very high
c --- frequency. for example, if baclin is 360 seconds, then setting
c --- nfladv to 10 will allow water properties to be sampled once per hour
c
c --- set variable nfldta to nfladv/2 so that it equals the number of time
c --- steps separating the velocity fields input into the runga-kutta
c --- interpolation
c
nfldta=nfladv/2
c
c --- make sure that the float will be advected at an integer number of
c --- temporal points per day
nsmday=nbcday/nfladv
if (int(real(nbcday)/real(2*nfldta)).ne.nsmday) then
if (mnproc.eq.1) then
write(lp,'(/ a /)')
& 'error - need integer no. of adv. times per day'
call flush(lp)
endif !1st tile
call xcstop('(floats_init)')
stop '(floats_init)'
endif
c
c --- test to make sure that the advection time interval is between 2
c --- and 4 hours. only stop the program later if there are floats to
c --- be advected
if (nfladv.ne.4 .and. (nsmday.gt.12 .or. nsmday.lt.6)) then
nfladv=-nfladv
endif
c
c --- initialize number of time steps between float output
if (nflsam.gt.0) then
nflsam=nflsam*nbcday
elseif (nflsam.eq.0) then
nflsam=abs(nfladv)
endif
c
c --- calculate the time interval for the runga-kutta interpolation
c --- as 1/2 of the advection time interval
deltfl=nfldta*baclin
c
c --- set minimum float depth (m)
fldepm=1.0
c
c --- -------------------------------------------------------------
c --- initialize synthetic drifters, floats, and moored instruments
c --- -------------------------------------------------------------
c
open(unit=uoff+99,file=trim(flnminp)//flnmflti,status='old') !on all nodes
c
c --- first line in float input file is the number of floats
read(uoff+99,*) nflt
if (nflt.gt.nfldim) then
if (mnproc.eq.1) then
write(lp,'(/ a /)')
& 'error - increase nfldim in mod_floats.F'
call flush(lp)
endif !1st tile
call xcstop('(floats_init)')
stop '(floats_init)'
endif
if (nflt.eq.0) then
synflt=.false.
endif
c
c --- second line in float input file is the initial float time step
c --- this number is set to zero for the first model run segment
read(uoff+99,*) nstepfl
c
do nfl=1,nflt
c
c --- read input file containing the following information for each float:
c ---
c --- column 1 - initial sequential float number
c --- column 2 - float type
c --- 1 = 3-d lagrangian (vertically advected by diagnosed w)
c --- 2 = isopycnic
c --- 3 = isobaric (surface drifter when released in sfc. layer)
c --- 4 = stationary (synthetic moored instrument)
c --- column 3 - deployment time (days from model start, 0.0 = immediate)
c --- column 4 - termination time (days from model start, 0.0 = forever)
c --- column 5 - initial longitude (must be between minimum and maximum
c --- longitudes defined in regional.grid.b)
c --- column 6 - initial latitude (must be between minimum and maximum
c --- latitudes defined in regional.grid.b)
c --- column 7 - initial depth (or reference sigma for isopycnic floats)
c
read(uoff+99,*) iflnum(nfl),ifltyp(nfl),flt(nfl,8),flt(nfl,9),
& flt(nfl,1),flt(nfl,2),flt(nfl,3)
if (ifltyp(nfl).eq.1 .or. ifltyp(nfl).eq.4) then
wvelfl=.true.
endif
if (ifltyp(nfl).ne.4) then
hadvfl=.true.
c
c --- since this float is to be advected, stop if the advection time interval
c --- is not between 2 and 4 hr
if (nfladv.lt.0) then
if (mnproc.eq.1) then
write(lp,'(/ a /)')
& 'error - set nfladv to advect floats every 2-4 hr'
call flush(lp)
endif !1st tile
call xcstop('(floats_init)')
stop '(floats_init)'
endif
c
endif
c
flt(nfl,4)=0.0
flt(nfl,5)=0.0
flt(nfl,6)=0.0
flt(nfl,7)=0.0
if(ifltyp(nfl).eq.2) then
flt(nfl,7)=flt(nfl,3)
flt(nfl,3)=0.0
endif
flt(nfl,10)=0.0
flt(nfl,11)=0.0
flt(nfl,12)=0.0
flt(nfl,13)=0.0
c
kfloat(nfl)=-1
c
enddo
c
nfladv=abs(nfladv)
c
close(unit=uoff+99) !file='float.input'
c
if (mnproc.eq.1) then
write(lp,955) nflt
write(lp,956) nfladv
write(lp,957) nflsam
write(lp,958) min(nflt,10)
955 format(' read initial data for',i6,' floats, time step',i9)
956 format(' float advected every',i7,' baroclinic time steps')
957 format(' float sampled every',i8,' baroclinic time steps'/)
958 format(' input data for first',i3,' floats')
do nfl=1,min(nflt,10)
write(lp,959) nfl,(flt(nfl,i),i=1,9)
959 format(3x,i2,1p,5e13.5/5x,4e13.5/)
enddo
call flush(lp)
endif !1st tile
c
c ---------------------------------------------------------------------------
c --- allocate arrays for floats
c ---------------------------------------------------------------------------
c
allocate( wvelup(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kdm),
& wveldn(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kdm),
& dpdxup(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy),
& dpdyup(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy),
& dpdxdn(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy),
& dpdydn(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) )
allocate( wveli( 1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy, kdm+1),
& uold2( 1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,2*kdm),
& vold2( 1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,2*kdm),
& wold2u(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,2*kdm),
& wold2d(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,2*kdm),
& dlondx(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy),
& dlondy(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy),
& dlatdx(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy),
& dlatdy(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) )
call mem_stat_add( 11*(idm+2*nbdy)*(jdm+2*nbdy)*kdm )
call mem_stat_add( 9*(idm+2*nbdy)*(jdm+2*nbdy) )
c
if (synflt) then
!$OMP PARALLEL DO PRIVATE(j,i,k)
!$OMP& SCHEDULE(STATIC,jblk)
do j=1-nbdy,jj+nbdy
do i=1-nbdy,ii+nbdy
do k=1,kk
uold2( i,j,k )=hugel
uold2( i,j,k+kk)=hugel
vold2( i,j,k )=hugel
vold2( i,j,k+kk)=hugel
wold2u(i,j,k )=hugel
wold2u(i,j,k+kk)=hugel
wold2d(i,j,k )=hugel
wold2d(i,j,k+kk)=hugel
wveli( i,j,k) =hugel
enddo !k
wveli(i,j,kk+1)=hugel
enddo
enddo
!$OMP END PARALLEL DO
endif !synflt
c
c ---------------------------------------------------------------------------
c --- initialize old horizontal velocities for runga-kutta time interpolation
c --- calculate dlondx, dlondy, dlatdx, dlatdy required to convert u, v to
c --- longitude/time and latitude/time for float advection
c ---------------------------------------------------------------------------
c
margin = nbdy - 1
c
do j=1-margin,jj+margin
do i=1-margin,ii+margin
if (SEA_U) then
do k=1,kk
#if defined(STOKES)
uold2(i,j,k )=u(i,j,k,n)+usd(i,j,k)+ubavg(i,j,n)
uold2(i,j,k+kk)=u(i,j,k,n)+usd(i,j,k)+ubavg(i,j,n)
#else
uold2(i,j,k )=u(i,j,k,n)+ubavg(i,j,n)
uold2(i,j,k+kk)=u(i,j,k,n)+ubavg(i,j,n)
#endif
enddo !k
dlondx(i,j)=(plon(i ,j)*scpy(i ,j)-
& plon(i-1,j)*scpy(i-1,j))/scu2(i,j)
dlatdx(i,j)=(plat(i ,j)*scpy(i ,j)-
& plat(i-1,j)*scpy(i-1,j))/scu2(i,j)
if (plat(i,j)-plat(i-1,j).ne.0.0) then
nonlatlon=.true.
endif
endif !iu
if (SEA_V) then
do k=1,kk
#if defined(STOKES)
vold2(i,j,k )=v(i,j,k,n)+vsd(i,j,k)+vbavg(i,j,n)+vsdbavg(i,j)
vold2(i,j,k+kk)=v(i,j,k,n)+vsd(i,j,k)+vbavg(i,j,n)+vsdbavg(i,j)
#else
vold2(i,j,k )=v(i,j,k,n)+vbavg(i,j,n)
vold2(i,j,k+kk)=v(i,j,k,n)+vbavg(i,j,n)
#endif
enddo !k
dlondy(i,j)=(plon(i,j )*scpx(i,j )-
& plon(i,j-1)*scpx(i,j-1))/scv2(i,j)
dlatdy(i,j)=(plat(i,j )*scpx(i,j )-
& plat(i,j-1)*scpx(i,j-1))/scv2(i,j)
if (plon(i,j)-plon(i,j-1).ne.0.0) then
nonlatlon=.true.
endif
endif !iv
if (SEA_P) then
do k=1,kk
wold2u(i,j,k )=0.0
wold2d(i,j,k )=0.0
wold2u(i,j,k+kk)=0.0
wold2d(i,j,k+kk)=0.0
wveli (i,j,k )=0.0
enddo !k
wveli(i,j,kk+1)=0.0
endif !ip
enddo !i
enddo !j
c
c ----------------------------------------
c --- turbulent horizontal velocity option
c ----------------------------------------
c
c --- initialize turbulence constants
c
c --- tdecri is in inverse days and dtturb is the turbulent time step in days
c
c --- dtturb is set to the runga-kutta interpolation time step, which must
c --- be within 1 and 2 hours
c
if (turbvel) then
dtturb=deltfl/86400.0
uturb0=sqrt(2.*tbvar*tdecri*dtturb)
endif
c
return
end subroutine floats_init
subroutine floats_restart
implicit none
c
c -----------------------------
c --- output float restart file
c -----------------------------
c
integer i,k,l
c
if (mnproc.eq.1) then
c
c --- first remove floats that have run aground or left the domain
i=0
do l=1,nflt
if (flt(l,1).gt.-999.0) then
i=i+1
endif
enddo !1
c
write(lp,*) 'writing new float input file for restart'
call flush(lp)
c
open(unit=uoff+802,file=flnmfltio,
& form='formatted',status='new')
write(uoff+802,970) i
write(uoff+802,970) nstepfl
do l=1,nflt
if (flt(l,1).gt.-999.0) then
if (ifltyp(l).ne.2) then
write(uoff+802,971) iflnum(l),ifltyp(l),flt(l,8),flt(l,9),
& (flt(l,k),k=1,3)
else
write(uoff+802,971) iflnum(l),ifltyp(l),flt(l,8),flt(l,9),
& (flt(l,k),k=1,2),flt(l,7)
endif
endif
enddo !l
close(uoff+802)
endif !1st tile
970 format(i9)
971 format(i6,i2,2f9.2,3f15.9)
return
end subroutine floats_restart
subroutine floats(m,n,timefl,ioflag)
use mod_cb_arrays ! HYCOM saved arrays
implicit none
c
integer, parameter :: nfl_debug = -1 !no debugging
c integer, parameter :: nfl_debug = 6 !debug float nfl_debug
c
integer m,n,ioflag
real timefl
c
c --- local variables
integer i,j,k,l,margin
real tflt,sflt,thlo,thhi,thflt,
& uflt1,uflt2,uflt3,uflt4,uflt,
& vflt1,vflt2,vflt3,vflt4,vflt,
& wflt1,wflt2,wflt3,wflt4,wflt,
& vkflt,tkflt,skflt,xpos0,xpos1,xpos2,xpos3,
& ypos0,ypos1,ypos2,ypos3,depflt,plo,phi,q,
& depiso,wvhi,wvlo,uvfctr,qisop,ufltm,vfltm,
& alomin,alomax,alamin,alamax
real uturb,vturb,uturb1,vturb1,dlodx,dlody,dladx,dlady
real*4 rtab(200,2)
integer ifltll(3),jfltll(3),ngood(3)
integer kflt,k1,k2,k3,ier,ngrid,kold1w,kold2w,kold1,kold2,
& ntermn,i1,j1,ngoodi
integer ied,iede,ifladv
c
integer*4 seed,numran,inisee,iflag
integer*4 iseed1,iseed2
equivalence ( iseed1, iseed2, seed )
c
c --- local velocity component storage for synthetic moorings
real fltloc(nflt,3)
c
c --- 2-d local arrays used for the 16-point box horizontal interpolation
real varb2d(4,4),ptlon(4,4,3),ptlat(4,4,3)
logical maskpt(4,4,3),maskpi(4,4)
c
c --- arrays required for parallelization
real flt1(12*nflt),flt3(17*nflt)
integer iproc
c
real timer1, timer2, totime
real*8 wtime
c
include 'stmt_fns.h'
c
c -------------------------------------------------------------------------
c --- floats.f (hycom version 2.2)
c -------------------------------------------------------------------------
c
c --- synthetic floats, drifters, and mooring instruments
c
c --- optionally samples time series of dynamical/thermodynamical variables
c --- at the location of each float
c
c --- four float types are presently supported:
c --- 3-d lagrangian (vertically advected by diagnosed vertical velocity)
c --- isopycnic (remains on specified density surface)
c --- isobaric (remains at the pressure depth where it was released)
c --- stationary (synthetic instrument/mooring)
c
c --- horizontal advection, along with vertical advection of lagrangian floats,
c --- is performed using the MICOM runga-kutta-4 time interpolation algorithm
c --- developed by Zulema Garraffo
c
c --- works on any curvilinear grid
c
c --- float position is stored as longitude and latitude
c
c --- to horizontal advect the float, u and v are first converted to
c --- d(longitude)/dt and d(latitude)/dt as follows:
c
c --- u_lon = u*dlondx + v*dlondy
c --- v_lat = u*dlatdx + v*dlatdy
c
c --- to horizontally advect the float, terms u*dlondx and u*dlatdx are
c --- calculated on u grid points while terms v*dlondy and v*dlatdy are
c --- calculated on v grid ponts. these terms are then interpolated to
c --- the float locations from the surrounding u and v grid points,
c --- respectively
c
c --- since the terms v*dlondy and u*dlatdx are always zero when the model
c --- x,y axes are everywhere lines of constant latitude and longitude,
c --- respectively, logical variable 'nonlatlon' prevents horizontal
c --- interpolation of these terms in this case
c
c --- horizontal interpolation to the float location follows the MICOM
c --- procedure of Zulema Garraffo - second-order interpolation from the 16
c --- surrounding grid points is performed unless fewer than nptmin (presently
c --- set to 10 in subroutine intrph) water points are available, in which
c --- case nearest-neighbor interpolation is used.
c
c --- variables are horizontally interpolated from their native grid (p, u,
c --- or v)
c
c --- adapted for hycom by george halliwell
c --- code parallelized by remy baraille
c
c --- the second-order interpolation subroutine included here is the one
c --- included in the mariano and brown parameter matrix objective analysis
c --- algorithm to estimate the large scale trend surface - it is different
c --- from the algorithm used by Garraffo in MICOM
c
c --- float initialization is performed in floats_init - information
c --- about the input file containing initial float information is
c --- presented in that subroutine.
c
c --- variables in array flt(nfl,n) are:
c
c --- n = 1 longitude
c --- n = 2 latitude
c --- n = 3 float depth
c --- n = 4 water depth
c --- n = 5 temperature
c --- n = 6 salinity
c --- n = 7 water density
c --- n = 8 float start time (in days from start of model run)
c --- n = 9 float end time (in days from start of model run)
c
c --- time series of several variables are interpolated to the location of each
c --- float and saved every nflsam time steps when ioflag is set to 1
c
c --- variables output onto file float.out for each float are:
c
c --- 1. initial sequential float number
c --- 2. time (in days from start of model run)
c --- 3. model layer number that contains the float
c --- 4. longitude (u for synthetic moorings)
c --- 5. latitude (v for synthetic moorings)
c --- 6. float depth (w for synthetic moorings)
c --- 7. water depth
c --- 8. temperature
c --- 9. salinity
c --- 10. water density (remains constant for isopycnic floats)
c
call xctilr(dp( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_ps)
call xctilr(dpu( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_vs)
call xctilr(dpv( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_vs)
call xctilr(u( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_uv)
call xctilr(v( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_vv)
call xctilr(ubavg( 1-nbdy,1-nbdy, n),1, 1, 6,6, halo_uv)
call xctilr(vbavg( 1-nbdy,1-nbdy, n),1, 1, 6,6, halo_vv)
call xctilr(temp( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_ps)
call xctilr(saln( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_ps)
call xctilr(th3d( 1-nbdy,1-nbdy,1,n),1,kk, 6,6, halo_ps)
c
c --- calculate vertical velocity at interfaces as the sum of the
c --- vertical interface velocity estimated in cnuity.f and the
c --- advective component due to flow parallel to interfaces.
c --- wvelup and wveldn represent velocity at the top and bottom of
c --- layer k
c
margin = 6
c
c --- set pressure array at p points
do j=1-margin,jj+margin
do i=1-margin,ii+margin
if (SEA_P) then
do k=1,kk
k1=k+1
p(i,j,k+1)=p(i,j,k)+dp(i,j,k,n)
enddo !k
endif !ip
enddo !i
enddo !j
c
c --- set pressure array at u and v points
c --- calculate dpdx, dpdy at interfaces above and below layer k for
c --- estimating wvel within layer k
c
do k=1,kk
c
margin = 5
c
do j=1-margin,jj+margin
do i=1-margin,ii+margin
if (SEA_U) then
pu(i,j,k+1)=pu(i,j,k)+dpu(i,j,k,n)
if (wvelfl) then
dpdxup(i,j)=
& (p(i,j,k )*scpy(i,j)-p(i-1,j,k )*scpy(i-1,j))
& /scu2(i,j)
dpdxdn(i,j)=
& (p(i,j,k+1)*scpy(i,j)-p(i-1,j,k+1)*scpy(i-1,j))
& /scu2(i,j)
endif
endif !iu
if (SEA_V) then
pv(i,j,k+1)=pv(i,j,k)+dpv(i,j,k,n)
if (wvelfl) then
dpdyup(i,j)=
& (p(i,j,k )*scpx(i,j)-p(i,j-1,k )*scpx(i,j-1))
& /scv2(i,j)
dpdydn(i,j)=
& (p(i,j,k+1)*scpx(i,j)-p(i,j-1,k+1)*scpx(i,j-1))
& /scv2(i,j)
endif
endif !iv
enddo !i
enddo !j
c
c --- calculate the vertical velocity arrays
c
margin = 4
c
do j=1-margin,jj+margin
do i=1-margin,ii+margin
if (SEA_P) then
if (wvelfl) then
wveli(i,j,k+1)=deltfl*wveli(i,j,k+1)/delt1
if (dp(i,j,k,n).gt.tencm) then
wvelup(i,j,k)=-0.5*deltfl*
#if defined(STOKES)
& ((u(i ,j ,k,n)+usd(i ,j ,k)+
& ubavg(i ,j ,n))*dpdxup(i ,j )+
& (u(i+1,j ,k,n)+usd(i+1,j ,k)+
& ubavg(i+1,j ,n))*dpdxup(i+1,j )+
& (v(i ,j ,k,n)+vsd(i ,j ,k)+
& vbavg(i ,j ,n))*dpdyup(i ,j )+
& (v(i ,j+1,k,n)+vsd(i ,j+1,k)+
& vbavg(i ,j+1,n))*dpdyup(i ,j+1))
#else
& ((u(i ,j ,k,n)+ubavg(i ,j ,n))*dpdxup(i ,j )+
& (u(i+1,j ,k,n)+ ubavg(i+1,j ,n))*dpdxup(i+1,j )+
& (v(i ,j ,k,n)+vbavg(i ,j ,n) )*dpdyup(i ,j )+
& (v(i ,j+1,k,n)+vbavg(i ,j+1,n))*dpdyup(i ,j+1))
#endif
& /onem+wveli(i,j,k )
wveldn(i,j,k)=-0.5*deltfl*
#if defined(STOKES)
& ((u(i ,j ,k,n)+usd(i ,j ,k)+
& ubavg(i ,j ,n))*dpdxdn(i ,j )+
& (u(i+1,j ,k,n)+usd(i+1,j ,k)+
& ubavg(i+1,j,n))*dpdxdn(i+1,j )+
& (v(i ,j ,k,n)+vsd(i ,j ,k)+
& vbavg(i ,j ,n))*dpdydn(i ,j )+
& (v(i ,j+1,k,n)+vsd(i ,j+1,k)+
& vbavg(i ,j+1,n))*dpdydn(i ,j+1))
#else
& ((u(i ,j ,k,n)+ubavg(i ,j ,n))*dpdxdn(i ,j )+
& (u(i+1,j ,k,n)+ubavg(i+1,j ,n))*dpdxdn(i+1,j )+
& (v(i ,j ,k,n)+vbavg(i ,j ,n))*dpdydn(i ,j )+
& (v(i ,j+1,k,n)+vbavg(i ,j+1,n))*dpdydn(i ,j+1))
#endif
& /onem+wveli(i,j,k+1)
else
wvelup(i,j,k)=0.0
wveldn(i,j,k)=0.0
endif
endif
endif !ip
enddo !i
enddo !j
c
enddo !k
c
* call xctilr(wvelup(1-nbdy,1-nbdy,1),1, kk, 4,4, halo_ps)
* call xctilr(wveldn(1-nbdy,1-nbdy,1),1, kk, 4,4, halo_ps)
c
c --- set old velocity indices used for interpolation of float position
c --- perform full float update every two times this subroutine is called;
c --- at the intermediate times, just store old velocity fields for the next
c --- call
c
if (mod(nstepfl,nfladv).eq.nfldta) then
kold1 =kk
kold1w=kk+1
kold2 =0
kold2w=0
go to 10
else
kold1 =0
kold1w=0
kold2 =kk
kold2w=kk+1
endif
c
c --- turbulent horizontal velocity option
c
c --- to customize the parameter settings in blkdat.input that control
c --- the calculated turbulent velocity (tbvar, tdecri), refer to
c --- Griffa (1996), "Aplications of stochastic particles models to
c --- oceanographic problems" in "Stochatic Modelling in Physical
c --- Oceanography", pp. 114-140, Adler, Muller, and Rozovoskii, editors
c
if (turbvel) then
c
c --- initialize random seeds and create random number table
c --- if iflag=1 (iflag=0), then time() is (is not) used to generate seeds
c
iflag=1
call system_clock(inisee)
seed = 414957000-inisee
iede = iseed1
ied = iseed2
if (iede .lt. 0) iede = abs(iede)
if (ied .lt. 0) ied = abs(ied)
c
call rantab_ini(rtab,iseed1,iseed2,iflag)
endif
c
c --- get particle locations from processor 1
if (mnproc.eq.1) then
do nfl=1,nflt
do i=1,9
flt1( i+12*(nfl-1)) = flt(nfl,i)
enddo
flt1(10+12*(nfl-1)) = kfloat(nfl)
flt1(11+12*(nfl-1)) = iflnum(nfl)
flt1(12+12*(nfl-1)) = ifltyp(nfl)
enddo
endif !1st tile
call xcastr(flt1(1:12*nflt), 1)
if (mnproc.ne.1) then
do nfl=1,nflt
do i=1,9
flt(nfl,i) = flt1( i+12*(nfl-1))
enddo
kfloat(nfl) = flt1(10+12*(nfl-1))
iflnum(nfl) = flt1(11+12*(nfl-1))
ifltyp(nfl) = flt1(12+12*(nfl-1))
enddo
endif !1st tile
c
c ----------------------
c --- begin float loop
c ----------------------
c
do nfl=1,nflt
c
c --- skip if float has previously run aground or exceeded termination time
if(flt(nfl,1).le.-999.) then
go to 22
endif
c
c --- ier = error flag for horizontal interpolation
c --- ntermn = float termination flag
c --- -10 => reached termination time
c --- -5 => exited domain
c --- >0 => ran aground
c --- ifladv = float horizontal advection flag
ier=0
ntermn=0
ifladv=1
c
c --- suppress float advection if this is the first time step so that
c --- the initial position on the output file is exactly the position
c --- specified on the input file
if (nstepfl.eq.1) then
ifladv=0
endif
c
c --- check if model time has reached float deployment time
c --- when it does, again suppress float advection during first pass
if (timefl-flt(nfl,8).lt.-0.001) then
go to 22
endif
if (kfloat(nfl).lt.0) then
ifladv=0
endif
c
c --- check if model time has reached float termination time
if ( flt(nfl,9).gt.0.0 .and.
& timefl-flt(nfl,9).gt.0.001 ) then
ntermn=-10
endif
c
c ---------------------------------------------------------------------------
c --- search for the processor controlling the tile containing the grid point
c --- search for the surrounding 16 grid points on the p, u, and v grids
c ---------------------------------------------------------------------------
c
c --- processeur sur lequel va se derouler le calcul
c
iproc=0
c
margin = 0
c
do j=1-margin,jj+margin
do i=1-margin,ii+margin
if (SEA_P) then
c
c --- search for the surrounding p-grid points
if (i+i0.lt.itdm .and. j+j0.lt.jtdm) then
alomin=min(plon(i ,j),plon(i ,j+1),
& plon(i+1,j),plon(i+1,j+1))
alomax=max(plon(i ,j),plon(i ,j+1),
& plon(i+1,j),plon(i+1,j+1))
alamin=min(plat(i ,j),plat(i ,j+1),
& plat(i+1,j),plat(i+1,j+1))
alamax=max(plat(i ,j),plat(i ,j+1),
& plat(i+1,j),plat(i+1,j+1))
if (flt(nfl,1).ge.alomin.and.flt(nfl,1).lt.alomax.and.
& flt(nfl,2).ge.alamin.and.flt(nfl,2).lt.alamax) then
ifltll(1)=i
jfltll(1)=j
c
c --- determine if float has exited domain
if (i+i0.lt.2 .or. i+i0.gt.itdm-1 .or.
& j+j0.lt.2 .or. j+j0.gt.jtdm-1) then
ntermn=-5
endif
c
c --- search for the surrounding u-grid points
do i1=i-1,i+1
do j1=j-1,j+1
if (i1+i0.lt.itdm .and. j1+j0.lt.jtdm) then
alomin=min(ulon(i1 ,j1),ulon(i1 ,j1+1),
& ulon(i1+1,j1),ulon(i1+1,j1+1))
alomax=max(ulon(i1 ,j1),ulon(i1 ,j1+1),
& ulon(i1+1,j1),ulon(i1+1,j1+1))
alamin=min(ulat(i1 ,j1),ulat(i1 ,j1+1),
& ulat(i1+1,j1),ulat(i1+1,j1+1))
alamax=max(ulat(i1 ,j1),ulat(i1 ,j1+1),
& ulat(i1+1,j1),ulat(i1+1,j1+1))
if (flt(nfl,1).ge.alomin .and.
& flt(nfl,1).lt.alomax .and.
& flt(nfl,2).ge.alamin .and.
& flt(nfl,2).lt.alamax) then
ifltll(2)=i1
jfltll(2)=j1
c
c --- determine if float has exited domain
if (i1+i0.lt.2 .or. i1+i0.gt.itdm-1 .or.
& j1+j0.lt.2 .or. j1+j0.gt.jtdm-1) then
ntermn=-5
endif
go to 11
endif
endif
enddo !j1
enddo !i1
11 continue
c
c --- search for the surrounding v-grid points
do i1=i-1,i+1
do j1=j-1,j+1
if (i1+i0.lt.itdm .and. j1+j0.lt.jtdm) then
alomin=min(vlon(i1 ,j1),vlon(i1 ,j1+1),
& vlon(i1+1,j1),vlon(i1+1,j1+1))
alomax=max(vlon(i1 ,j1),vlon(i1 ,j1+1),
& vlon(i1+1,j1),vlon(i1+1,j1+1))
alamin=min(vlat(i1 ,j1),vlat(i1 ,j1+1),
& vlat(i1+1,j1),vlat(i1+1,j1+1))
alamax=max(vlat(i1 ,j1),vlat(i1 ,j1+1),
& vlat(i1+1,j1),vlat(i1+1,j1+1))
if (flt(nfl,1).ge.alomin .and.
& flt(nfl,1).lt.alomax .and.
& flt(nfl,2).ge.alamin .and.
& flt(nfl,2).lt.alamax) then
ifltll(3)=i1
jfltll(3)=j1
c
c --- determine if float has exited domain
if (i1+i0.lt.2 .or. i1+i0.gt.itdm-1 .or.
& j1+j0.lt.2 .or. j1+j0.gt.jtdm-1) then
ntermn=-5
endif
go to 12
endif
endif
enddo !j1
enddo !i1
12 continue
c
c --- set processor number for the tile containing the float and exit grid
c --- point search
iproc=mnproc
go to 13
endif
endif
endif !ip
enddo !i
enddo !j
c
13 continue
c
c --- float nfl is now updated by the processor running the tile containing
c --- the float
c
if (iproc.gt.0) then
c
if (nfl.eq.nfl_debug) then
write(lp,100) nfl,nflt,nstep,nstepfl
write(lp,101) nfl,ntermn,flt(nfl,1),flt(nfl,2),
& (ifltll(i)+i0,jfltll(i)+j0,i=1,3),
& mnproc,
& plon(ifltll(1),jfltll(1)),
& plat(ifltll(1),jfltll(1))
100 format(/'diagnostics for float',i6,' of',i6,', time step',
& i9/'float time step',i9)
101 format('float',i6,' ntermn:',i6,' position:',2(1pe12.4)/
& 'lower left points, p,u,v:',3(2i5,2x)/
& 'mnproc =',i4,' plon,plat =',1pe12.4,1pe12.4)
call flush(lp)
endif !nfl_debug
c
c --- if float has exited domain or run aground as determined previously,
c --- jump ahead to the float termination code bloci
if (ntermn.eq.-5 .or. ntermn.eq.1) then
go to 30
endif
c
c --- set ptlon, ptlat for all horizontal interpolations from each grid
ngrid=1
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
ptlon(i1,j1,ngrid)=plon(i,j)
ptlat(i1,j1,ngrid)=plat(i,j)
enddo
enddo
c
ngrid=2
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
ptlon(i1,j1,ngrid)=ulon(i,j)
ptlat(i1,j1,ngrid)=ulat(i,j)
enddo
enddo
c
ngrid=3
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
ptlon(i1,j1,ngrid)=vlon(i,j)
ptlat(i1,j1,ngrid)=vlat(i,j)
enddo
enddo
c
c --- the float is assumed to remain within the same layer that it was in
c --- during the previous update unless it is the first advection time step
c --- for the float, in which case it is initially assumed to be in layer 1
c
k=max(1,kfloat(nfl))
c
c --- mask p grid points that are on land or where the layer containing the
c --- float is a zero or nearly-zero thickness layer at the bottom
ngrid=1
ngood(ngrid)=0
do i1=1,4
i=i1+ifltll(ngrid)-2
do j1=1,4
j=j1+jfltll(ngrid)-2
if (depths(i,j).lt.0.01 .or. depths(i,j).eq.hugel .or.
& depths(i,j)*onem-p(i,j,k).lt.tencm) then
maskpt(i1,j1,ngrid)=.true.
else
maskpt(i1,j1,ngrid)=.false.
ngood(ngrid)=ngood(ngrid)+1
endif
enddo
enddo
c
if (nfl.eq.nfl_debug) then
write(lp,102) ngood(1)
102 format('initial masking: no. of good p points',i4)
endif !nfl_debug
c