-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathnswing.c
5726 lines (5118 loc) · 230 KB
/
nswing.c
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
/*--------------------------------------------------------------------
* $Id: nswing.c 9936 2016-11-17 17:57:28Z j $
*
* Copyright (c) 2012-2015 by J. Luis and J. M. Miranda
*
* This program is part of Mirone and is free software; you can redistribute
* it and/or modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* Contact info: w3.ualg.pt/~jluis/mirone
*--------------------------------------------------------------------*/
static char prog_id[] = "$Id: nswing.c 9936 2016-11-17 17:57:28Z j $";
/*
* Original Fortran version of core hydrodynamic code by J.M. Miranda and COMCOT
*
*
* The nesting grids algorithm
LOOP_OVER_ALL_CYCLES {
mass_L0
openbd
LOOP_OVER_N_NESTED {
LOOP_TIME_REFINE_1 {
interp_edges_L1
mass_L1
LOOP_TIME_REFINE_2 {
interp_edges_L2
mass_L2
LOOP_TIME_REFINE_3 {
...
}
moment_L2
upscale_L2
update_L2
}
moment_L1
upscale_L1
update_L1
}
}
moment_L0
update_L0
}
*
*
* To compile, do for example
* cl nswing.c -IC:\programs\compa_libs\netcdf_GIT\compileds\VC12_64\include
* C:\programs\compa_libs\netcdf_GIT\compileds\VC12_64\lib\netcdf.lib /DI_AM_C
* /DHAVE_NETCDF /nologo /D_CRT_SECURE_NO_WARNINGS /fp:precise /Ox
*
* Rewritten in C, mexified, added number options, etc... By
* Joaquim Luis - 2013
*
*/
#if defined(WIN32) || defined(_WIN32) || defined(_WIN64)
# define DO_MULTI_THREAD /* ISTO TEM DE SER AUTOMATIZADO, OU VIA COMPILA */
#else
# define strtok_s strtok_r
#endif
#define I_AM_MEX /* Build as a MEX */
#ifdef I_AM_C /* Build as a stand-alone exe */
# undef I_AM_MEX
#endif
#include <float.h>
#include <math.h>
#include <string.h>
#include <stdint.h>
#include <time.h>
#ifdef I_AM_MEX
# include "mex.h"
# include "mwsize.h"
#else
# include <stdio.h>
# include <stdlib.h>
# include <stdint.h>
# define mxCalloc calloc
# define mxMalloc malloc
# define mxRealloc realloc
# define mxFree free
# define mexPrintf(...) fprintf(stderr, __VA_ARGS__);
# ifndef NAN
# ifdef _MSC_VER
# include <ymath.h>
# define NAN _Nan._Double
# else
static const double _NAN = 20;//(HUGE_VAL-HUGE_VAL);
# define NAN _NAN
# endif
# endif
# define mxGetNaN() (NAN)
#endif
/* At least the Intel compiler screws anf the NAN macro returns 0.0 So we need this hack */
/* http://stackoverflow.com/questions/5714131/nan-literal-in-c */
union {uint64_t i; double d;} loc_nan = {0x7ff8000000000000};
#ifdef HAVE_NETCDF
# include <netcdf.h>
# define err_trap(status) if (status) {mexPrintf ("NSWING: error at line: %d\t and errorcode = %s\n", __LINE__, nc_strerror(status));}
#endif
#if defined(WIN32) || defined(_WIN32) || defined(_WIN64)
# include <windows.h>
# include <process.h>
#endif
#if HAVE_OPENMP
#include <omp.h>
#endif
#define FALSE 0
#define TRUE 1
#ifndef M_PI
# define M_PI 3.14159265358979323846
#endif
#ifndef M_PI_2
# define M_PI_2 1.57079632679489661923
#endif
#define D2R M_PI / 180.
#define R2D 180. / M_PI
#define GMT_CONV_LIMIT 1.0e-8 /* Fairly tight convergence limit or "close to zero" limit */
#define NORMAL_GRAV 9.806199203 /* Moritz's 1980 IGF value for gravity at 45 degrees latitude */
#define EQ_RAD 6378137.0 /* WGS-84 */
#define flattening 1.0/298.2572235630
#define LIMIT_DISCHARGE /* If defined the moment functions will limit the discharge */
#define ECC2 2 * flattening - flattening * flattening
#define ECC4 ECC2 * ECC2
#define ECC6 ECC2 * ECC4
#define EPS12 1e-12
#define EPS10 1e-10
#define EPS6 1e-6
#define EPS5 1e-5
#define EPS4_ 1e-4
#define EPS3 1e-3
#define EPS2 1e-2
#define EPS1 1e-1
static double EPS4 = EPS4_; /* Kinda trick to be able to change EPS4 via a command line option */
#define MAXRUNUP -50 /* Do not waste time computing flood above this altitude */
#define V_LIMIT 20 /* Upper limit of maximum velocity */
#define CNULL ((char *)NULL)
#define Loc_copysign(x,y) ((y) < 0.0 ? -fabs(x) : fabs(x))
#ifndef MIN
# define MIN(x, y) (((x) < (y)) ? (x) : (y)) /* min and max value macros */
#endif
#ifndef MAX
# define MAX(x, y) (((x) > (y)) ? (x) : (y))
#endif
#ifndef rint
# define rint(x) (floor((x)+0.5))
#endif
#ifndef irint
# define irint(x) ((int)rint(x))
#endif
#define ijs(i,j,n) ((i) + (j)*n)
#define ijc(i,j) ((i) + (j)*n_ptmar)
#define ij_grd(col,row,hdr) ((col) + (row)*hdr.nx)
typedef void (*PFV) (); /* PFV declares a pointer to a function returning void */
struct tracers { /* For tracers (oranges) */
double *x; /* x coordinate */
double *y; /* y coordinate */
};
struct srf_header { /* Surfer file hdr structure */
char id[4]; /* ASCII Binary identifier (DSAA/DSBB) */
short int nx; /* Number of columns */
short int ny; /* Number of rows */
double x_min; /* Minimum x coordinate */
double x_max; /* Maximum x coordinate */
double y_min; /* Minimum y coordinate */
double y_max; /* Maximum y coordinate */
double z_min; /* Minimum z value */
double z_max; /* Maximum z value */
};
struct grd_header { /* Generic grid hdr structure */
int nx; /* Number of columns */
int ny; /* Number of rows */
unsigned int nm; /* nx * ny */
double x_inc;
double y_inc;
double x_min; /* Minimum x coordinate */
double x_max; /* Maximum x coordinate */
double y_min; /* Minimum y coordinate */
double y_max; /* Maximum y coordinate */
double z_min; /* Minimum z value */
double z_max; /* Maximum z value */
int doCoriolis; /* Apply Coriolis if this != 0 */
double lat_min4Coriolis; /* PRECISA SOLUCAO. POR AGORA SERA Cte = 0 */
};
struct nestContainer { /* Container for the nestings */
int do_upscale; /* If false, do not upscale the parent grid */
int do_long_beach; /* If true, compute a mask with ones over the "dryed beach" */
int do_short_beach; /* If true, compute a mask with ones over the "innundated beach" */
int do_linear; /* If true, use linear approximation */
int do_max_level; /* If true, inform nestify() on the need to update max level at every inner iteration */
int do_max_velocity; /* If true, inform nestify() on the need to update max velocity at each inner iteration */
int do_Coriolis; /* If true, compute the Coriolis effect */
int out_velocity_x; /* To know if we must compute the vex,vey velocity arrays */
int out_velocity_y;
int out_momentum; /* To know if save the momentum in the 3D netCDF grid. Mutually exclusive with out_velocity_x|y */
int isGeog; /* 0 == Cartesian, otherwise Geographic coordinates */
int writeLevel; /* Store info about which level is (if) to be writen [0] */
int bnc_pos_nPts; /* Number of points in a external boundary condition file */
int bnc_var_nTimes; /* Number of time steps in the external boundary condition file */
int bnc_border[4]; /* Each will be set to TRUE if boundary condition on that border W->0, S->1, E->2, N->3 */
int level[10]; /* 0 Will mean base level, others the nesting level */
int LLrow[10], LLcol[10], ULrow[10], ULcol[10], URrow[10], URcol[10], LRrow[10], LRcol[10];
int incRatio[10];
short *long_beach[10]; /* Mask arrays for storing the "dry beaches" */
short *short_beach[10]; /* Mask arrays for storing the "dry beaches" */
float *work, *wmax; /* Auxiliary pointers (not direcly allocated) to compute max level of nested grids */
float *vmax; /* Pointer to array storing the max velocity */
double run_jump_time; /* Time to hold before letting the nested grids start to iterate */
double lat_min4Coriolis; /* South latitute when computing the Coriolis effect on a cartesian grid */
double manning_depth; /* Do not use manning if depth is deeper than this value */
double manning[10]; /* Manning coefficient. Set to zero if no friction */
double LLx[10], LLy[10], ULx[10], ULy[10], URx[10], URy[10], LRx[10], LRy[10];
double dt[10]; /* Time step at current level */
double *bat[10]; /* Bathymetry of current level */
double *fluxm_a[10], *fluxm_d[10]; /* t-1/2 & t+1/2 fluxes arrays along X */
double *fluxn_a[10], *fluxn_d[10]; /* t-1/2 & t+1/2 fluxes arrays along Y */
double *htotal_a[10], *htotal_d[10]; /* t-1/2 & t+1/2 total water depth */
double *vex[10], *vey[10]; /* X,Y velocity components */
double *etaa[10], *etad[10]; /* t-1/2 & t+1/2 water height (eta) arrays */
double *edge_col[10], *edge_colTmp[10];
double *edge_row[10], *edge_rowTmp[10];
double *edge_col_P[10], *edge_col_Ptmp[10];
double *edge_row_P[10], *edge_row_Ptmp[10];
double *r0[10], *r1m[10], *r1n[10], *r2m[10], *r2n[10], *r3m[10], *r3n[10], *r4m[10], *r4n[10];
double time_h;
double *bnc_pos_x;
double *bnc_pos_y;
double *bnc_var_t;
double **bnc_var_z;
double *bnc_var_zTmp;
double *bnc_var_z_interp;
struct grd_header hdr[10];
};
/* Argument struct for threading */
typedef struct {
struct nestContainer *nest; /* Pointer to a nestContainer struct */
int iThread; /* Thread index */
int lev; /* Level of nested grid */
} ThreadArg;
void no_sys_mem(char *where, unsigned int n);
int count_col(char *line);
int read_grd_info_ascii(char *file, struct srf_header *hdr);
int read_header_bin (FILE *fp, struct srf_header *hdr);
int write_grd_bin(char *name, double x_min, double y_min, double x_inc, double y_inc, unsigned int i_start,
unsigned int j_start, unsigned int i_end, unsigned int j_end, unsigned int nX, float *work);
int read_grd_ascii (char *file, struct srf_header *hdr, double *work, int sign);
int read_grd_bin(char *file, struct srf_header *hdr, double *work, int sign);
int read_maregs(struct grd_header hdr, char *file, unsigned int *lcum_p, char *names[]);
int read_tracers(struct grd_header hdr, char *file, struct tracers *oranges);
int count_n_maregs(char *file);
int decode_R(char *item, double *w, double *e, double *s, double *n);
int check_region(double w, double e, double s, double n);
double ddmmss_to_degree (char *text);
void openb(struct grd_header hdr, double *bat, double *fluxm_d, double *fluxn_d, double *etad, struct nestContainer *nest);
void wave_maker(struct nestContainer *nest);
void wall_it(struct nestContainer *nest);
void wall_two(struct nestContainer *nest, int ot1, int ot2, int in1, int in2);
int initialize_nestum(struct nestContainer *nest, int isGeog, int lev);
int intp_lin (double *x, double *y, int n, int m, double *u, double *v);
void inisp(struct nestContainer *nest);
void inicart(struct nestContainer *nest);
void interp_edges(struct nestContainer *nest, double *flux_L1, double *flux_L2, char *what, int lev, int i_time);
void sanitize_nestContainer(struct nestContainer *nest);
void nestify(struct nestContainer *nest, int nNg, int recursionLevel, int isGeog);
void resamplegrid(struct nestContainer *nest, int nNg);
void edge_communication(struct nestContainer *nest, int lev, int i_time);
void mass(struct nestContainer *nest, int lev);
void mass_sp(struct nestContainer *nest, int lev);
void mass_conservation(struct nestContainer *nest, int isGeog, int m);
void moment_conservation(struct nestContainer *nest, int isGeog, int m);
void update(struct nestContainer *nest, int lev);
void upscale(struct nestContainer *nest, double *out, int lev, int i_tsr);
void upscale_(struct nestContainer *nest, double *out, int lev, int i_tsr);
void replicate(struct nestContainer *nest, int lev);
void moment_M(struct nestContainer *nest, int lev);
void moment_N(struct nestContainer *nest, int lev);
void moment_sp_M(struct nestContainer *nest, int lev);
void moment_sp_N(struct nestContainer *nest, int lev);
void free_arrays(struct nestContainer *nest, int isGeog, int lev);
int check_paternity(struct nestContainer *nest);
int check_binning(double x0P, double x0D, double dxP, double dxD, double tol, double *suggest);
int read_bnc_file(struct nestContainer *nest, char *file);
int interp_bnc (struct nestContainer *nest, double t);
void total_energy(struct nestContainer *nest, float *work, int lev);
void power(struct nestContainer *nest, float *work, int lev);
void vtm (double lat0, double *t_c1, double *t_c2, double *t_c3, double *t_c4, double *t_e2, double *t_M0);
void deform (struct srf_header hdr, double x_inc, double y_inc, int isGeog, double fault_length,
double fault_width, double th, double dip, double rake, double d, double top_depth,
double xl, double yl, double *z);
void kaba_source(struct srf_header hdr, double x_inc, double y_inc, double x_min, double x_max,
double y_min, double y_max, int type, double *z);
void tm (double lon, double lat, double *x, double *y, double central_meridian, double t_c1,
double t_c2, double t_c3, double t_c4, double t_e2, double t_M0);
double uscal(double x1, double x2, double x3, double c, double cc, double dp);
double udcal(double x1, double x2, double x3, double c, double cc, double dp);
unsigned int gmt_bcr_prep (struct grd_header hdr, double xx, double yy, double wx[], double wy[]);
double GMT_get_bcr_z(double *grd, struct grd_header hdr, double xx, double yy);
void update_max(struct nestContainer *nest);
void update_max_velocity(struct nestContainer *nest);
#ifdef HAVE_NETCDF
void write_most_slice(struct nestContainer *nest, int *ncid_most, int *ids_most, unsigned int i_start,
unsigned int j_start, unsigned int i_end, unsigned int j_end, float *work, size_t *start,
size_t *count, double *slice_range, int isMost, int lev);
int open_most_nc(struct nestContainer *nest, float *work, char *basename, char *name_var, char hist[], int *ids,
unsigned int nx, unsigned int ny, double xMinOut, double yMinOut, int isMost, int lev);
int open_anuga_sww(struct nestContainer *nest, char *fname_sww, char history[], int *ids, unsigned int i_start,
unsigned int j_start, unsigned int i_end, unsigned int j_end, double xMinOut, double yMinOut, int lev);
void write_anuga_slice(struct nestContainer *nest, int ncid, int z_id, unsigned int i_start, unsigned int j_start,
unsigned int i_end, unsigned int j_end, float *work, size_t *start, size_t *count,
float *slice_range, int idx, int with_land, int lev);
int write_maregs_nc(struct nestContainer *nest, char *fname, float *work, double *t, unsigned int *lcum_p,
char *names[], char hist[], int n_maregs, unsigned int count_time_maregs_timeout, int lev);
int write_greens_nc(struct nestContainer *nest, char *fname, float *work, size_t *start, size_t *count,
double *t, unsigned int *lcum_p, char *names[], char hist[], int *ids, int n_maregs,
unsigned int n_times, int lev);
void err_trap_(int status);
#endif
#if defined(WIN32) || defined(_WIN32) || defined(_WIN64)
/* Prototypes for threading related functions */
unsigned __stdcall MT_cart(void *Arg_p);
unsigned __stdcall MT_sp(void *Arg_p);
int GetLocalNThread(void);
#endif
/* Function pointers to M & N moment components */
PFV call_moment[2];
PFV call_moment_sp[2];
int Return(int code) { /* To handle return codes between MEX and standalone code */
#ifdef I_AM_MEX
mexErrMsgTxt("\n"); /* Most of the cases (no_sys_mem) this instruction was executed already */
#endif
return(code);
}
/* --------------------------------------------------------------------------- */
/* Matlab Gateway routine */
#ifdef I_AM_MEX
# define Return(code) {mexErrMsgTxt("\n");}
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
#else
# define Return(code) {return(code);}
int main(int argc, char **argv) {
#endif
int writeLevel = 0; /* If save grids, will hold the saving level (when nesting) */
int i, j, k, n;
int start_i; /* Where the loop over argc starts: MEX -> 0; STANDALONE -> 1 */
int grn = 0, cumint = 0, decimate_max = 1, iprc, r_bin_b, r_bin_f, r_bin_mM, r_bin_mN;
int w_bin = TRUE, cumpt = FALSE, error = FALSE, do_2Dgrids = FALSE, do_maxs = FALSE;
int out_energy = FALSE, max_energy = FALSE, out_power = FALSE, max_power = FALSE;
int first_anuga_time = TRUE, out_sww = FALSE, out_most = FALSE, out_3D = FALSE;
int surf_level = TRUE, max_level = FALSE, max_velocity = FALSE, water_depth = FALSE;
int do_Okada = FALSE; /* For when one will compute the Okada deformation here */
int do_Kaba = FALSE; /* For when one will use prismatic sources */
int do_tracers = FALSE; /* For when doing Lagrangian tracers */
int out_maregs_nc = FALSE; /* For when maregs in output are written in netCDF */
int out_oranges_nc = FALSE; /* For when tracers in output are written in netCDF */
int do_HotStart = FALSE; /* For when doing a Hot Start */
int n_arg_no_char = 0;
int ncid, ncid_most[3], z_id = -1, ids[13], ids_ha[6], ids_ua[6], ids_va[6], ids_most[3];
int ncid_3D[3], ids_z[10], ids_3D[3], ncid_Mar, ids_Mar[8];
int n_of_cycles = 1010; /* Default number of cycles to compute */
int num_of_nestGrids = 0; /* Number of nesting grids */
int bat_in_input = FALSE, source_in_input = FALSE, write_grids = FALSE, isGeog = FALSE;
int maregs_in_input = FALSE, out_momentum = FALSE, got_R = FALSE;
int with_land = FALSE, IamCompiled = FALSE, do_nestum = FALSE, saveNested = FALSE, verbose = FALSE;
int out_velocity = FALSE, out_velocity_x = FALSE, out_velocity_y = FALSE, out_velocity_r = FALSE;
int out_maregs_velocity = FALSE;
int KbGridCols = 1, KbGridRows = 1; /* Number of rows & columns IF computing a grid of 'Kabas' */
int cntKabas = 0; /* Counter of the number of Kabas (prisms) already processed */
int n_mareg, n_ptmar, n_oranges, pos_prhs;
unsigned int *lcum_p = NULL, lcum = 0, ij, nx, ny;
unsigned int i_start, j_start, i_end, j_end, count_maregs_timeout = 0, count_time_maregs_timeout = 0;
size_t start0 = 0, count0 = 1, len, start1_A[2] = {0,0}, count1_A[2];
size_t start1_M[3] = {0,0,0}, count1_M[3], start_Mar[3] = {0,0,0}, count_Mar[2];
char *bathy = NULL; /* Name pointer for bathymetry file */
char hcum[256] = ""; /* Name of the cumulative hight file */
char maregs[256] = ""; /* Name of the maregraph positions file */
char **mareg_names = NULL; /* Array to hold optional maregraph names */
char *fname_sww = NULL; /* Name pointer for Anuga's .sww file */
char *basename_most = NULL; /* Name pointer for basename of MOST .nc files */
char *fname3D = NULL; /* Name pointer for the 3D netCDF file */
char *fonte = NULL; /* Name pointer for tsunami source file */
char *bnc_file = NULL; /* Name pointer for a boundary condition file */
char fname_mask_lbeach[256] = ""; /* Name pointer for the "long_beach" mask grid */
char fname_mask_sbeach[256] = ""; /* Name pointer for the "short_beach" mask grid */
char tracers_infile[256] = "", tracers_outfile[256] = ""; /* Names for in and out tracers files */
char stem[256] = "", prenome[128] = "", str_tmp[128] = "", fname_momentM[256] = "", fname_momentN[256] = "";
char history[512] = {""}; /* To hold the full command call to be saved in nc files as History */
char *pch;
char *nesteds[10] = {NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL};
char txt[128]; /* Auxiliary variable */
float *work = NULL, *workMax = NULL, *vmax = NULL, *wmax = NULL, *time_p = NULL;
float work_min = FLT_MAX, work_max = -FLT_MAX, *maregs_array = NULL, *maregs_array_t = NULL;
double *maregs_timeout = NULL, m_per_deg = 111317.1;
double *bat = NULL, *dep1 = NULL, *dep2 = NULL, *cum_p = NULL, *h = NULL;
double dfXmin = 0, dfYmin = 0, dfXmax = 0, dfYmax = 0, xMinOut, yMinOut;
double kaba_xmin = 0, kaba_xmax = 0, kaba_ymin = 0, kaba_ymax = 0;
double time_jump = 0, time0, time_for_anuga, prc;
double dt = 0; /* Time step for Base level grid */
double dx, dy, ds, dtCFL, etam, one_100, t;
double *eta_for_maregs, *vx_for_maregs, *vy_for_maregs, *htotal_for_maregs, *fluxm_for_maregs, *fluxn_for_maregs;
double *vx_for_oranges, *vy_for_oranges, *fluxm_for_oranges, *fluxn_for_oranges, *htotal_for_oranges; /* For tracers */
double f_dip, f_azim, f_rake, f_slip, f_length, f_width, f_topDepth, x_epic, y_epic; /* For Okada initial condition */
double add_const = 0, time_h = 0;
double dxKb = 0, dyKb = 0; /* Grid steps for when computing a grid of 'Kabas' */
double z_offset = 0; /* To apply to bathymetry to simulate a tide */
double manning[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; /* Manning coefficients */
double actual_range[6] = {1e30, -1e30, 1e30, -1e30, 1e30, -1e30};
float stage_range[2], xmom_range[2], ymom_range[2], *tmp_slice;
struct srf_header hdr_b, hdr_f, hdr_mM, hdr_mN;
struct grd_header hdr;
struct nestContainer nest;
struct tracers *oranges;
FILE *fp, *fp_oranges;
#ifdef I_AM_MEX
int argc;
unsigned nm;
char **argv, cmd[16] = "";
double *head, *tmp, x_inc, y_inc, x, y;
double *ptr_wb; /* Pointer to be used in the aguentabar */
mxArray *rhs[3];
#endif
clock_t tic;
call_moment[0] = (PFV)moment_M;
call_moment[1] = (PFV)moment_N;
call_moment_sp[0] = (PFV)moment_sp_M;
call_moment_sp[1] = (PFV)moment_sp_N;
#ifdef DO_MULTI_THREAD
if ((k = GetLocalNThread()) == 1) {
mexPrintf("NSWING: This version of the program is build for multi-threading but "
"this machine has only one core. Don't know what will happen here.\n");
}
#endif
sanitize_nestContainer(&nest);
#ifdef I_AM_MEX
argc = nrhs;
for (i = 0; i < nrhs; i++) { /* Check input to find how many arguments are of type char */
if (!mxIsChar(prhs[i])) {
argc--;
n_arg_no_char++; /* Number of arguments that have a type other than char */
}
}
if (n_arg_no_char >= 4) { /* Not all combinations are tested */
/* Bathymetry */
dep1 = mxGetPr(prhs[0]);
nx = mxGetN(prhs[0]);
ny = mxGetM(prhs[0]);
head = mxGetPr(prhs[1]); /* Get bathymetry header info */
hdr_b.x_min = head[0]; hdr_b.x_max = head[1];
hdr_b.y_min = head[2]; hdr_b.y_max = head[3];
hdr_b.z_min = head[4]; hdr_b.z_max = head[5];
hdr_b.nx = nx; hdr_b.ny = ny;
x_inc = head[7]; y_inc = head[8];
bat_in_input = TRUE;
/* Tsunami source */
h = mxGetPr(prhs[2]);
nx = mxGetN(prhs[2]);
ny = mxGetM(prhs[2]);
head = mxGetPr(prhs[3]); /* Get bathymetry header info */
hdr_f.x_min = head[0]; hdr_f.x_max = head[1];
hdr_f.y_min = head[2]; hdr_f.y_max = head[3];
hdr_f.z_min = head[4]; hdr_f.z_max = head[5];
hdr_f.nx = nx; hdr_f.ny = ny;
source_in_input = TRUE;
/* Now test if bat & source are compatible. If they are not, we go out right here. */
if (hdr_f.nx != hdr_b.nx || hdr_f.ny != hdr_b.ny) {
mexPrintf("Bathymetry & Source grids have different nx and/or ny\n");
mexPrintf("%d %d %d %d\n", hdr_f.nx, hdr_b.nx, hdr_f.ny, hdr_b.ny);
return;
}
}
if (n_arg_no_char >= 5 && mxIsCell(prhs[4])) {
const mxArray *mx_ptr;
num_of_nestGrids = mxGetNumberOfElements(prhs[4]);
if (num_of_nestGrids % 2 != 0)
mexErrMsgTxt("NSWING: Cell array with nested grids is not a even number of elements");
num_of_nestGrids /= 2;
for (k = 0; k < num_of_nestGrids; k++) {
mx_ptr = mxGetCell(prhs[4], k);
if (mx_ptr == NULL)
mexErrMsgTxt("NSWING: bloody cell element is empty!!!!!!!");
dep2 = mxGetPr(mx_ptr); /* Get the grid for this level */
nest.hdr[k+1].nx = mxGetN(mx_ptr);
nest.hdr[k+1].ny = mxGetM(mx_ptr);
nest.hdr[k+1].nm = (unsigned int)nest.hdr[k+1].nx * (unsigned int)nest.hdr[k+1].ny;
mx_ptr = mxGetCell(prhs[4], k + num_of_nestGrids);
head = mxGetData(mx_ptr); /* Get header info */
nest.hdr[k+1].x_min = head[0]; nest.hdr[k+1].x_max = head[1];
nest.hdr[k+1].y_min = head[2]; nest.hdr[k+1].y_max = head[3];
nest.hdr[k+1].z_min = head[4]; nest.hdr[k+1].z_max = head[5];
nest.hdr[k+1].x_inc = head[7]; nest.hdr[k+1].y_inc = head[8];
nm = nest.hdr[k+1].nx * nest.hdr[k+1].ny;
if ((nest.bat[k+1] = (double *)mxCalloc((size_t)nm, sizeof(double)) ) == NULL)
{no_sys_mem("(bat)", nm); Return(-1);}
for (i = 0; i < nest.hdr[k+1].ny; i++) {
for (j = 0; j < nest.hdr[k+1].nx; j++)
nest.bat[k+1][i*nest.hdr[k+1].nx+j] = -dep2[j*nest.hdr[k+1].ny+i];
}
}
do_nestum = TRUE;
}
/* See if we have a vector of maregraphs */
pos_prhs = 0;
if (n_arg_no_char == 5 && !mxIsEmpty(prhs[4]))
pos_prhs = 4;
else if (n_arg_no_char == 6 && !mxIsEmpty(prhs[5]))
pos_prhs = 5;
if (pos_prhs) {
double x_min, x_max, y_min, y_max;
tmp = mxGetPr(prhs[pos_prhs]);
n_mareg = mxGetM(prhs[pos_prhs]);
if (do_nestum) {
dx = nest.hdr[num_of_nestGrids].x_inc; dy = nest.hdr[num_of_nestGrids].y_inc;
x_min = nest.hdr[num_of_nestGrids].x_min; x_max = nest.hdr[num_of_nestGrids].x_max;
y_min = nest.hdr[num_of_nestGrids].y_min; y_max = nest.hdr[num_of_nestGrids].y_max;
nx = nest.hdr[num_of_nestGrids].nx;
}
else {
dx = x_inc; dy = x_inc;
x_min = hdr_b.x_min; x_max = hdr_b.x_max;
y_min = hdr_b.y_min; y_max = hdr_b.y_max;
nx = hdr_b.nx;
}
lcum_p = (unsigned int *)mxCalloc((size_t)(n_mareg), sizeof(unsigned int));
mareg_names = mxCalloc((size_t)(n_mareg), sizeof(char *));
for (ij = j = 0; ij < n_mareg; ij++) {
x = tmp[ij]; y = tmp[ij+n_mareg]; /* Matlab vectors are stored by columns */
if (x < x_min || x > x_max || y < y_min || y > y_max)
continue;
lcum_p[j] = (unsigned int)(irint((y - y_min) / dy) ) * nx + (unsigned int)irint((x - x_min) / dx);
j++;
}
if (j > 0) {
maregs_in_input = TRUE;
cumpt = TRUE;
}
else {
mexPrintf("NSWING: Warning, none of the provided maregraph locations lies inside (inner?) grid. Ignoring them\n");
mxFree(lcum_p);
lcum_p = NULL;
}
}
/* get the length of the input string */
argv = (char **)mxCalloc(argc, sizeof(char *));
for (i = 0; i < argc; i++)
argv[i] = (char *)mxArrayToString(prhs[i+n_arg_no_char]);
start_i = 0;
#else
start_i = 1; /* For stand-alone first index is prog name, so we want to start at 1 */
#endif /* end #ifdef I_AM_MEX */
for (i = start_i; i < argc; i++) {
if (argv[i][0] == '-') {
switch (argv[i][1]) {
case 'c':
add_const = atof(&argv[i][2]);
break;
case 'e':
IamCompiled = TRUE;
break;
case 'f': /* */
isGeog = TRUE;
break;
case 'n': /* Write MOST files (*.nc) */
basename_most = &argv[i][2];
out_most = TRUE;
break;
case 'A': /* Name for the Anuga .sww netCDF file */
fname_sww = &argv[i][2];
out_sww = TRUE;
if (argv[i][2] == 'l') /* Output land nodes in SWW file */
with_land = TRUE;
break;
case 'B': /* File with the boundary condition (experimental) */
bnc_file = &argv[i][2];
break;
case 'C': /* Output momentum grids */
nest.do_Coriolis = TRUE;
if (argv[i][2])
nest.lat_min4Coriolis = atof(&argv[i][2]);
break;
case 'D':
water_depth = TRUE;
surf_level = FALSE;
break;
case 'E': /* Compute total Energy or Power*/
if (strstr(argv[i], "EPS4=")) { /* Secreet option to change the EPS4 value */
EPS4 = atof(&argv[i][6]);
break;
}
if (argv[i][2] == 'p') {
if (argv[i][3] == 'm')
max_power = TRUE;
else
out_power = TRUE;
}
else {
if (argv[i][2] == 'm')
max_energy = TRUE;
else
out_energy = TRUE;
}
if ((pch = strstr(&argv[i][2],",")) != NULL) {
decimate_max = atoi(++pch);
if (decimate_max < 1) decimate_max = 1;
}
break;
case 'F': /* Okada parameters to compute Initial condition */
if (argv[i][2] == 'k') {
char *lost_str1 = NULL, *lost_str2 = NULL;
int have_RC = FALSE;
do_Kaba = TRUE;
k = 3;
if (argv[i][3] == 'c') {
k = 4;
do_Kaba++; /* Signal kaba_source function that the above is actually x/y/nx/ny */
}
n = sscanf(&argv[i][k], "%lf/%lf/%lf/%lf/%s/%lf", &kaba_xmin, &kaba_xmax, &kaba_ymin, &kaba_ymax, txt, &dyKb);
if (n > 4) {
lost_str1 = strdup(txt); /* Save it to restore argv[i] for use in the nc History attrib */
pch = strchr(txt, 'x');
if (pch != NULL) {
KbGridCols = atoi(&pch[1]);
pch[0] = '\0'; /* Strip the 'x...' part */
KbGridRows = atoi(txt);
have_RC = TRUE;
}
else {
dxKb = atof(txt);
if (n == 5) dyKb = dxKb;
}
/* Now strip the 5th (&6th) fields from argv[i] before calling decode_R */
if (n == 6) {
pch = strrchr(argv[i], '/');
lost_str2 = strdup(pch); /* Iden lost_str1 */
pch[0] = '\0';
}
pch = strrchr(argv[i], '/'); pch[0] = '\0';
}
k -= 2; /* Decrease k because the decode_R() function expects 2 chars at begining of string (-R...) */
error += decode_R(&argv[i][k], &kaba_xmin, &kaba_xmax, &kaba_ymin, &kaba_ymax);
if (have_RC) {
dxKb = (kaba_xmax - kaba_xmin); dyKb = (kaba_ymax - kaba_ymin);
}
if (dxKb != 0 && !have_RC) { /* Here dx/dy were given instead of RxC, so we have to compute them */
KbGridCols = irint((kaba_xmax - kaba_xmin) / dxKb);
KbGridRows = irint((kaba_ymax - kaba_ymin) / dyKb);
}
if (KbGridRows * KbGridCols > 1) out_maregs_nc = TRUE; /* Otherwise we can still save in ascii */
if (n > 4) {
strcat(argv[i], "/"); strcat(argv[i], lost_str1);
mxFree(lost_str1);
}
if (lost_str2) {
strcat(argv[i], lost_str2); /* Here the leading slash was still in lost_str2 */
mxFree(lost_str2);
}
}
else {
do_Okada = TRUE;
n = sscanf(&argv[i][2], "%lf/%lf/%lf/%lf/%lf/%lf/%lf/%lf/%lf",
&f_dip, &f_azim, &f_rake, &f_slip, &f_length, &f_width, &f_topDepth, &x_epic, &y_epic);
if (n != 9) {
mexPrintf("NSWING: Error, -F option, must provide all 9 parameters.\n");
error++;
}
else { /* Convert fault dimensions to meters (that's what is used by deform) */
f_length *= 1000;
f_width *= 1000;
f_topDepth *= 1000;
}
}
break;
case 'G': /* Write grids at grn intervals */
case 'Z': /* Write one single 3D netCDF at grn intervals */
sscanf(&argv[i][2], "%s,%d", stem, &grn);
if ((pch = strstr(stem,",")) != NULL) {
grn = atoi(&pch[1]);
pch[0] = '\0'; /* Strip the ",num" part */
}
if ((pch = strstr(stem,"+")) != NULL) {
writeLevel = atoi(pch++);
if (writeLevel < 0) writeLevel = 0;
pch--;
pch[0] = '\0'; /* Hide the +lev from the stem string */
saveNested = TRUE;
}
if (argv[i][1] == 'G')
write_grids = TRUE;
else {
out_3D = TRUE;
fname3D = stem;
if (fname3D[strlen(fname3D)-3] != '.' && fname3D[strlen(fname3D)-4] != '.')
strcat(fname3D, ".nc"); /* If no 2 or 3 letters extension, add .nc */
}
break;
case 'H': /* Hot start stuff */
if (!argv[i][2]) /* Output momentum grids */
out_momentum = TRUE;
else if (argv[i][2] == 's' && argv[i][3] == ',') /* NOT YET. Maybe it will be -Hs,time_to_stop */
;
else {
sscanf(&argv[i][2], "%s", str_tmp);
if ((pch = strstr(str_tmp,",")) != NULL) {
pch[0] = '\0';
strcpy(fname_momentM, str_tmp); /* File names of moment M & N files */
strcpy(fname_momentN, ++pch);
if ((pch = strstr(fname_momentN,",")) != NULL) { /* We have the hot start time as well */
pch[0] = '\0'; /* The end of the moment N file name */
time_h += atof(++pch); /* Increment the starting time (was zero) */
}
do_HotStart = TRUE;
}
else {
mexPrintf("NSWING: Error, -H option (Hot start), must provide names of moment_X, moment_Y files.\n");
error++;
}
}
break;
case 'J': /* Jumping options. Accept either -Jn, -J+m, -Jn+m or -Jn -J+m */
sscanf(&argv[i][2], "%s", str_tmp);
if ((pch = strstr(str_tmp,"+")) != NULL) {
sscanf((++pch), "%lf", &nest.run_jump_time);
pch--;
pch[0] = '\0'; /* Put the string end where before was the '+' char */
}
if (argv[i][2])
sscanf(&argv[i][2], "%lf", &time_jump);
break;
case 'L': /* Use linear approximation or Lagrangean tracers*/
if (!argv[i][2])
nest.do_linear = TRUE;
else {
sscanf(&argv[i][2], "%s", str_tmp);
if (str_tmp[strlen(str_tmp)-2] == '+') { /* Output tracers file will be in netCDF */
out_oranges_nc = TRUE;
str_tmp[strlen(str_tmp)-2] = '\0';
}
if ((pch = strstr(str_tmp,",")) != NULL) {
char *pch2;
pch[0] = '\0';
strcpy(tracers_infile, str_tmp);
strcpy(tracers_outfile, ++pch); /* NEED TO TEST IF WE GOT A FNAME */
}
else {
mexPrintf("NSWING: Error, -L option, must provide at least the tracers file name\n");
error++;
}
do_tracers = TRUE;
out_velocity_x = out_velocity_y = TRUE;
}
break;
case 'M': /* Max/min options */
if (argv[i][2] == '-') { /* Compute a mask with ones over the "dried beach" */
nest.do_long_beach = TRUE;
if (argv[i][3])
sscanf(&argv[i][3], "%s", fname_mask_lbeach);
else
strcpy(fname_mask_lbeach, "long_beach.grd");
}
else if (argv[i][2] == '+') { /* Compute a mask with ones over the "innundated beach" */
nest.do_short_beach = TRUE;
if (argv[i][3])
sscanf(&argv[i][3], "%s", fname_mask_sbeach);
else
strcpy(fname_mask_sbeach, "short_beach.grd");
}
else
max_level = TRUE;
break;
case 'N': /* Number of cycles to compute */
n_of_cycles = atoi(&argv[i][2]);
break;
case 'O': /* Time interval and fname for maregraph data. Use only when mareg locations were sent in as arg (MEX) */
sscanf(&argv[i][2], "%s", str_tmp);
if ((pch = strstr(str_tmp,",")) != NULL) {
pch[0] = '\0';
cumint = atoi(str_tmp);
strcpy(hcum, ++pch);
}
else {
mexPrintf("NSWING: Error, -O option, must provide interval and output maregs file name\n");
error++;
}
break;
case 'Q': /* Vertical offset (simulate tide) */
if (argv[i][2])
sscanf(&argv[i][2], "%lf", &z_offset);
break;
case 'R':
error += decode_R(argv[i], &dfXmin, &dfXmax, &dfYmin, &dfYmax);
got_R = TRUE;
break;
case 'S': /* Output velocity grids */
strcpy(str_tmp, &argv[i][2]);
if ((pch = strstr(str_tmp,"+m")) != NULL) { /* Velocity at maregraphs */
out_maregs_velocity = TRUE;
out_velocity_x = out_velocity_y = TRUE;
for (k = 0; k < 2; k++) { /* Strip the '+m' from the str_tmp string */
len = strlen(pch);
for (j = 0; j < (int)len; j++)
pch[j] = pch[j+1];
}
}
if ((pch = strstr(str_tmp,"+s")) != NULL) { /* Comput max speed */
max_velocity = TRUE;
out_velocity_x = out_velocity_y = TRUE;
for (k = 0; k < 2; k++) { /* Strip the '+s' from the str_tmp string */
len = strlen(pch);
for (j = 0; j < (int)len; j++)
pch[j] = pch[j+1];
}
}
if (str_tmp[0] == 'x') { /* Only X component */
out_velocity = out_velocity_x = TRUE;
if (!(out_maregs_velocity || max_velocity)) out_velocity_y = FALSE;
}
else if (str_tmp[0] == 'y') { /* Only Y component */
out_velocity = out_velocity_y = TRUE;
if (!(out_maregs_velocity || max_velocity)) out_velocity_x = FALSE;
}
else if (str_tmp[0] == 'r') { /* Speed (velocity module) -- NOT YET -- */
out_velocity = out_velocity_r = TRUE;
if (!(out_maregs_velocity || max_velocity)) out_velocity_x = out_velocity_y = FALSE;
}
else if (str_tmp[0] == 'n') { /* No grids, only vx,vy in maregs */
out_velocity_x = out_velocity_y = TRUE;
out_velocity = FALSE;
}
else { /* Both X & Y*/
out_velocity = out_velocity_x = out_velocity_y = TRUE;
}
break;
case 't': /* Time step of simulation */
dt = atof(&argv[i][2]);
nest.dt[0] = dt;
break;
case 'T': /* File with time interval (n steps), maregraph positions and optional output fname */
if (cumpt) {
mexPrintf("NSWING: Error, this option is not to be used when maregraphs were transmitted in input\n");
mexPrintf(" Ignoring it.\n");
break;
}
sscanf(&argv[i][2], "%s", str_tmp);
if (str_tmp[strlen(str_tmp)-2] == '+') { /* Output maregs file will be in netCDF */
out_maregs_nc = TRUE;
str_tmp[strlen(str_tmp)-2] = '\0';
}
if ((pch = strstr(str_tmp,",")) != NULL) {
char *pch2;
pch[0] = '\0';
if ((pch2 = strstr(str_tmp,".")) != NULL)
mexPrintf("NSWING: WARNING, 'int' in option -T<int> must be an integer number. Expect surprises.\n");
cumint = atoi(str_tmp);
if ((pch2 = strstr(++pch,",")) != NULL) {
pch2[0] = '\0';
strcpy(maregs, pch);
strcpy(hcum, ++pch2);
}
else
strcpy(maregs, pch);
}
else {
mexPrintf("NSWING: Error, -T option, must provide at least a interval and maregs file name\n");
error++;
}
cumpt = TRUE;
maregs_in_input = FALSE;
#ifndef HAVE_NETCDF
if (out_maregs_nc) {
mexPrintf("NSWING: Error, -T cannot choose an netCDF output because this exe was not linked to netCDF.\n");
out_maregs_nc = FALSE;
}
#endif
break;
case 'X': /* Manning coeffs */
k = 0;
sscanf(&argv[i][2], "%s", str_tmp);
if ((pch = strstr(str_tmp,"+")) != NULL) {
nest.manning_depth = -atof(++pch); /* Reverse sense right away because bat will be pos down */
pch--; pch[0] = '\0'; /* Remove traces of this option in string */
}
if ((pch = strstr(str_tmp,",")) != NULL) {
char t[16] = "";
pch[0] = '\0';
nest.manning[k++] = atof(str_tmp);
strcpy(t, ++pch);
while ((pch = strstr(t,",")) != NULL) {
pch[0] = '\0';
nest.manning[k++] = atof(t);
strcpy(t, ++pch);
}
nest.manning[k] = atof(t); /* Last one doesn't end with a comma so was skipped above */
}
else
nest.manning[k] = atof(&argv[i][2]);
if (k == 0) /* Only one set. Replicate it to the others (being used or not) */
for (n = 1; n < 10; n++)
nest.manning[n] = nest.manning[0];
break;
case 'U':
nest.do_upscale = TRUE;
break;
case 'V':
verbose = TRUE;
break;
case '1':
nesteds[0] = &argv[i][2];
break;
case '2':
case '3':
case '4':
case '5':
case '6':
case '7':
case '8':
case '9':
case '10':
nesteds[atoi(&argv[i][1])-1] = &argv[i][2];
break;
default:
mexPrintf("NSWING: Unknown option %s\n", argv[i]);
error = TRUE;
break;
}
}
else {
if (argv[i][0] == ' ' && argv[i][1] == '\0') /* An empty option (not uncommon from the Matlab side) */
continue;
if (bathy == NULL)
bathy = argv[i];
else if (fonte == NULL)
fonte = argv[i];
else if (bathy != NULL && fonte != NULL) {
mexPrintf("NSWING: Wrong option %s (misses the minus sign)\n", argv[i]);
error = TRUE;
}
}
}
if (argc <= 1 || error) {