forked from chiaraholgate/QIBT_shared
-
Notifications
You must be signed in to change notification settings - Fork 2
/
BTrIMS.f90
4097 lines (3128 loc) · 150 KB
/
BTrIMS.f90
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
!%%%%%%% AUTHORSHIP %%%%%%%
! This program was written by Jason P. Evans (UNSW) and modified by Chiara M. Holgate (ANU).
! Please seek permission before publishing all or part of this code, from jason.evans@unsw.edu.au and/or chiara.holgate@anu.edu.au.
!%%%%%%% PURPOSE %%%%%%%
! This program calculates the quasi-isentropic back trajectories of water vapor using a method based on Dirmeyer & Brubaker, 1999, "Contrasting evaporative moisture sources during the drought of 1988 and the flood of 1993", Journal of Geophysical Research, 104 D16 pg 19,383-19,397.
!
!%%%%%%% INPUT %%%%%%%
! Input data are taken from NARCliM output.
! Data here:
! /srv/ccrc/data33/z3481416/CCRC-WRF3.6.0.5-SEB/ERA-Interim/R2_nudging/out/
! And copied to here: /g/data/hh5/tmp/w28/jpe561/back_traj/
!%%Model expects:%%
! % Data to range between 0deg and 360deg.
! % Rainfall: [mm], 3d.
! % Latent heat [W/m2], 3d, which is converted to evaporation by program.
! % Temperature: actual temperature [K], 4d.
! % Pressure: total pressure [Pa], 4d, where first vertical level is at the top of the model.
! % Surface pressure [Pa], 3d.
! % U,V wind speed [m/s] components, 4d. If you're using vertical wind speeds, you need W [m/s] too.
! % Water-equivalent variables: in wrf, this is QCLD, QRAIN, QSNOW, QICE, QVAPOR [kg/kg], 4d.
!%%%%%%% OUTPUT %%%%%%%
! The program outputs a daily 2d grid of water vapour contribution of each grid cell to each precipitation cell.
! It also outputs the total rainfall depth in each cell where it rained, and the coordinates of each rain grid cell.
!
! The output file is dimensioned (rec,lat,lon) where rec is the number of grid cells where it rained that day,
! lat & lon are the sizes of the domain in the lat and lon directions. Attributes include location and time and amount
! of precip in the pixel of interest (y,x,day,month,year,precip) each of which is a 1D array of length rec.
!
!%%%%%%% EXPLANATION OF TIME-RELATED VARIABLES %%%%%%%
! "set" = user defined values
! "calcd" = calculated within program
! totdays = number of days to run simulation forward, based on set start/end dates(calcd)
! totbtadays = number of days to back-track for (set)
! tstep = number of minutes for back-track time step (set)
! daytsteps = 1440/tstep = 48 = number of simulation time steps in a day (calcd)
! totsteps = daytsteps*(totbtadays+1) = total number of simulation time steps over period (calcd)
! datatstep = 180 = input file time step in minutes (set)
! datadaysteps = 1440/datatstep = 8 = number of input file time steps in a day (calcd)
! indatatsteps = datatstep/tstep = number of simulation time steps per input time step (calcd)
! datatotsteps = (datadaysteps*(totbtadays+1)) + 1 = total number of input time steps over the back-track period (calcd)
!.............we want the event day + totbtadays before it + the time step after ! (ie 0 hour time step which is the last in each file)
! ttdataday = = ((tt-1)/indatatsteps) + 1 = position of parcel time step in input file (calcd)
! ttdata = datatotsteps - datadaysteps - 1 + ttdataday = input file time step from the beginning of the loaded files (calcd)
!%%%%%%% ASSUMPTIONS %%%%%%%
! All forms of water are included (vapour, rain, ice, snow, cloud water) in the parcel's mixing ratio.
! The PBL is NOT split from the upper atmosphere, i.e. the whole column is assumed to be well-mixed at the time step scale.
! Height of parcel release is determined randomly from a humidity-weighted vertical profile, i.e the vertical distribution of water vapour indicates where the rain forms.
! Time of parcel release is determined randomly through a precipitation-weighted time profile.
! The only source of parcel moisture is surface evaporation.
! The only sink for parcel moisture is precipitation.
! Program assumes data ranges between 0deg and 360deg.
! If input data of different structure is to be used for this program, subroutines (e.g. get_data, get_data_mixtot) will need to be changed.
! Note that I have not checked the isentropic routines "implicit_back_traj" or potential temp routines. If you want to use that version of the program (instead of moving parcels vertically using w), then those routines will need to be thoroughly checked.
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MODULE global_data
IMPLICIT NONE
SAVE
!
!*******user modified variables**********************
!
INTEGER :: sday,smon,syear !start day for calculations
INTEGER :: edday,edmon,edyear !end day for calculations (Exclusive. Must be at least one day after start day)
INTEGER :: totdays
INTEGER, PARAMETER :: totbtadays = 15 !number of days of data to keep for bta; i.e. how far back in time to calc.
!must be less than days you have input data for
INTEGER, PARAMETER :: tstep = 15 !number of minutes for back trajectory time step (simultion time step)
!must divide evenly into number of minutes in day 1440 and number of minutes in MM5 time step (here 180)
INTEGER, PARAMETER :: nparcels = 10 !set the number of parcels to release if it rains
REAL, PARAMETER :: minpre = 1 !min daily precip to deal with (mm)
INTEGER, PARAMETER :: bdy = 6 !boundary layers to ignore; trajectories will be tracked to this boundary
CHARACTER(LEN=50), PARAMETER :: diri = "/g/data/hh5/tmp/w28/jpe561/back_traj/"
CHARACTER(LEN=50), PARAMETER :: diri_era5 = "/g/data/w28/jpe561/BTrIMS/"
! CHARACTER(LEN=50), PARAMETER :: diri = "/srv/ccrc/data03/z3131380/PartB/Masks/"
! CHARACTER(LEN=100), PARAMETER :: diro = "/g/data/xc0/user/Holgate/QIBT/exp02/"
CHARACTER(LEN=100) :: diro
CHARACTER(LEN=100), PARAMETER :: dirdata_atm = "/g/data/hh5/tmp/w28/jpe561/back_traj/wrfout/"
CHARACTER(LEN=100), PARAMETER :: dirdata_era5 = "/g/data/rt52/era5/"
CHARACTER(LEN=100), PARAMETER :: dirdata_land = "/g/data/hh5/tmp/w28/jpe561/back_traj/wrfhrly/"
! CHARACTER(LEN=100), PARAMETER :: dirdata_atm = "/srv/ccrc/data33/z3481416/CCRC-WRF3.6.0.5-SEB/ERA-Interim/R2_nudging/out/"
! CHARACTER(LEN=100), PARAMETER :: dirdata_land = "/srv/ccrc/data03/z3131380/PartB/NARCliM_postprocess/"
INTEGER, PARAMETER :: numthreads = 48 !set the number of parallel openmp threads
LOGICAL, PARAMETER :: peak = .FALSE. !does the daylist indicate storm peaks (TRUE) or whole days (FALSE)
LOGICAL, PARAMETER :: wshed = .TRUE. !only calculate trajectories for watershed
CHARACTER(LEN=50), PARAMETER :: fwshed = "NARCliM_AUS_land_sea_mask.nc"
CHARACTER(LEN=50), PARAMETER :: fwshed_era5 = "Pakistan_mask_int_to180.nc"
!set to "" if no watershed
!0 outside watershed, >0 inside
REAL, PARAMETER :: min_del_q = 0.0001 !the minimum change in parcel mixing ratio (kg/kg) to be considered "real"
REAL, PARAMETER :: delta_coord = 0.0001 ! 1/10000th degree - for floating point calculations
LOGICAL, PARAMETER :: eachParcel = .FALSE. !output the data along the trajectory of each parcel
!****************************************************
!
INTEGER :: daytsteps,totsteps,indatatsteps,datadaysteps,datatotsteps
INTEGER :: dim_i,dim_j,dim_k,fdim_i,fdim_j,ssdim
INTEGER :: dim_i_start, dim_j_start, dim_k_start
INTEGER :: dim_i_end, dim_j_end, dim_k_end
INTEGER :: mon,year,dd,totpts
INTEGER :: day
! Additional variable to define the number of time intervals in the input data, now that the input data is monthly instead of daily
INTEGER :: datansteps
REAL, PARAMETER :: Lv = 2.25E6 !latent heat of vaporization of water (Jkg-1)
REAL, PARAMETER :: g = 9.8 !gravity (m.s-2)
REAL, PARAMETER :: P0 = 100000 !reference surface pressure (Pa)
REAL, PARAMETER :: Rd = 287.053 !ideal gas constant for dry air (J/kgK)
REAL, PARAMETER :: Cp = 1004.67 !heat capacity of air at constant pressure (J/kgK)
REAL, PARAMETER :: Rv = 461.5 !gas constant of water vapor (J/kgK)
REAL, PARAMETER :: Cl = 4400 !heat capacity of liquid water at ~-20C (J/kgK)
REAL, PARAMETER :: pi = 3.14159265
REAL, PARAMETER :: deg_dist = 111. !average distance of 1 degree lat is assumed to be 111km
REAL, PARAMETER :: water_density = 1000 ! density of water (kg/m3)
END MODULE global_data
MODULE util
IMPLICIT NONE
CONTAINS
SUBROUTINE handle_err(status)
!---------------------------------
! handle any errors from the fortran 90 netCDF interface
!---------------------------------
USE netcdf
IMPLICIT NONE
integer, intent (in) :: status
print *,'status = ',status
if(status /= nf90_noerr) then
print *, trim(nf90_strerror(status))
stop "Stopped"
end if
END SUBROUTINE handle_err
!***********************************************************************
INTEGER FUNCTION string_to_int(string)
!------------------------------------------
! converts a character string to an integer
!--------------------------------------------
IMPLICIT NONE
character (len=*), intent(in) :: string
! local constant
integer, parameter :: zero = iachar("0")
integer :: i, sign, integ
character (len=50) :: str
str = trim(string)
integ = 0
select case (str(1:1))
case ("-")
sign = -1
str = str(2:)
case ("+")
sign = 1
str = str(2:)
case ("0":"9")
sign = 1
end select
do i=len(trim(str)),1,-1
select case (str(i:i))
case ("0":"9")
integ = integ + (iachar(string(i:i))-zero)*10**(len(trim(str))-i)
case default
print *, "cannot convert a non-integer character to an integer!!"
return
end select
end do
string_to_int = integ
end FUNCTION string_to_int
!***********************************************************************
CHARACTER(LEN=50) FUNCTION int_to_string(integ)
!----------------------------------------------
! converts an integer to a character string
!----------------------------------------------
IMPLICIT NONE
integer, intent(in) :: integ
! local constant
integer, parameter :: zero = iachar("0")
character (len=50) :: str
integer :: i, inte
str=" "
inte = integ
do i=1,50
str(50-i:50-i) = achar(mod(abs(inte),10)+zero)
inte = int(inte/10)
if (abs(inte)<1) exit
end do
if (integ<0) str(50-i-1:50-i) = "-"
int_to_string = adjustl(str)
end FUNCTION int_to_string
!***********************************************************************
CHARACTER(LEN=50) FUNCTION real_to_string(num)
!----------------------------------------------
! converts an real to a character string
! here I allow a maximum of 10 decimal places
!----------------------------------------------
IMPLICIT NONE
!real, intent(in) :: num
INTEGER, intent(in) :: num
! local
integer :: i, whole, frac
real :: fracnum
whole = num !AINT(num)
fracnum = num - whole
do i=1,10
if (MOD(fracnum,1.)==0) exit
fracnum = fracnum * 10.
end do
frac = AINT(fracnum)
real_to_string = TRIM(int_to_string(whole))//"."//TRIM(int_to_string(frac))
end FUNCTION real_to_string
!***********************************************************************
INTEGER FUNCTION julian(year,mon,day)
! From http://aa.usno.navy.mil/faq/docs/JD_Formula.php
IMPLICIT NONE
INTEGER, INTENT(IN) :: year,mon,day
julian= day-32075+1461*(year+4800+(mon-14)/12)/4+367*(mon-2-(mon-14)/12*12)/12-3*((year+4900+(mon-14)/12)/100)/4
END FUNCTION julian
!***********************************************************************
SUBROUTINE GREGORIAN(JD,YEAR,MONTH,DAY)
! From http://aa.usno.navy.mil/faq/docs/JD_Formula.php
IMPLICIT NONE
!REAL, INTENT(IN) :: JD
INTEGER :: JD
INTEGER, INTENT(OUT) :: YEAR, MONTH, DAY!, HOUR, MINUTE, SECOND
REAL :: JT
INTEGER :: I,J,K,L,N
L = INT(JD)+68569!JD= K-32075+1461*(I+4800+(J-14)/12)/4+367*(J-2-(J-14)/12*12)/12-3*((I+4900+(J-14)/12)/100)/4
N = 4*L/146097
L = L-(146097*N+3)/4
I = 4000*(L+1)/1461001
L = L-1461*I/4+31
J = 80*L/2447
K = L-2447*J/80
L = J/11
J = J+2-12*L
I = 100*(N-49)+I+L
YEAR = I
MONTH = J
DAY = K
! JT = DMOD(JD,1.D0)*24.D0
! HOUR = INT(JT)
! JT = DMOD(JT,1.D0)*60.D0
! MINUTE = INT(JT)
! JT = DMOD(JT,1.D0)*60.D0
! SECOND = NINT(JT)
!
! IF (SECOND == 60) THEN
! SECOND = SECOND-60
! MINUTE = MINUTE+1
! END IF
END SUBROUTINE GREGORIAN
!***********************************************************************
INTEGER FUNCTION simlength(startday,startmon,startyear,endday,endmon,endyear)
!-----------------------------------------------------
! given the start and end dates of the desired simulation
! period (top of global data), calculate the number of days
! in the period
! Functions: don't need to specify what comes out, e.g. fn_name(in,in,in)
! Subroutines: do specify in and out, e.g. sbrtn_name(in,in,in,out,out)
!-----------------------------------------------------
USE global_data
IMPLICIT NONE
INTEGER, intent(in) :: startday,startmon,startyear,endday,endmon,endyear
INTEGER :: start_jd, end_jd
start_jd = julian(startyear,startmon,startday)
end_jd = julian(endyear,endmon,endday)
simlength = end_jd - start_jd
END FUNCTION simlength
!***********************************************************************
SUBROUTINE day_month_year(simday)
!----------------------------------------------
! given the simulation day, calculate the corresponding month and year
! simulation day one is the first day
!-----------------------------------------------------
USE global_data
IMPLICIT NONE
INTEGER,INTENT(IN) :: simday !simulation day
INTEGER :: jday,simdayyear,simdaymonth,simdayday
if (simday==1) then
year=syear
mon=smon
day=sday
else
! Find julian day of simday, being the start day + an increment of days
jday=julian(syear,smon,sday)+(simday-1)
! Convert the julian day to a gregorian day
call gregorian(jday,simdayyear,simdaymonth,simdayday)
year=simdayyear
mon=simdaymonth
day=simdayday
end if
END SUBROUTINE day_month_year
!***********************************************************************
SUBROUTINE all_positive_longitude(lon2d,lon2d_corrected)
!------------------------------------------------
! NARCliM longitude 2d grid (XLONG) is negative east of dateline.
! This subroutine converts the negative values to positive.
!------------------------------------------------
REAL, DIMENSION(:,:) :: lon2d,lon2d_corrected
lon2d_corrected=lon2d
WHERE (lon2d < 0)
lon2d_corrected = lon2d+360
END WHERE
END SUBROUTINE all_positive_longitude
!***********************************************************************
FUNCTION to_iso_date(y,m,d)
!----------------------------------------------------
! Convert year month day integers to iso date string
!----------------------------------------------------
INTEGER, INTENT(IN) :: y,m,d
CHARACTER(len=8) :: to_iso_date
write(to_iso_date,'(i4.4,i2.2,i2.2)') y,m,d
END FUNCTION to_iso_date
FUNCTION month_end(y,m)
!----------------------------------------------------
! Lookup table for ends of months
!----------------------------------------------------
INTEGER, INTENT(IN) :: y,m
INTEGER :: month_end
SELECT CASE (m)
CASE(1,3,5,7,8,10,12)
month_end = 31
RETURN
CASE(4,6,9,11)
month_end = 30
RETURN
END SELECT
IF ( MODULO(y,4) == 0 ) then
IF ( MODULO(y,100) == 0 .and. MODULO(y,400) /= 0 ) then
month_end = 28
RETURN
ELSE
month_end = 29
RETURN
END IF
ELSE
month_end = 28
RETURN
END IF
END FUNCTION month_end
SUBROUTINE array_extents(arr,in_start,in_end,i_start,i_end,reverse,periodic)
USE global_data, ONLY: delta_coord
IMPLICIT NONE
REAL, DIMENSION(:), INTENT(IN) :: arr
REAL, INTENT(IN) :: in_start, in_end
INTEGER, INTENT(OUT) :: i_start, i_end
LOGICAL,OPTIONAL,INTENT(IN) :: reverse
LOGICAL,OPTIONAL,INTENT(IN) :: periodic
!!! Locals
INTEGER :: i
REAL :: start, end
i_start = -1
i_end = -1
if( .not.present(periodic) .or. .not.periodic ) then
if ( in_start > in_end ) then
!!! Swap bounds if necessary
end = in_start
start = in_end
else
start = in_start
end = in_end
end if
else
start = in_start
end = in_end
end if
if( present(reverse) .and. reverse ) then
do i=1,SIZE(arr)
if ( i_start == -1 ) then
if ( abs(arr(i) - end) < delta_coord ) i_start = i
else if ( i_end == -1 ) then
if ( abs(arr(i) - start) < delta_coord ) i_end = i
else
exit
end if
end do
else
do i=1,SIZE(arr)
if ( i_start == -1 ) then
if ( abs(arr(i) - start) < delta_coord ) i_start = i
else if ( i_end == -1 ) then
if ( abs(arr(i) - end) < delta_coord ) i_end = i
else
exit
end if
end do
end if
if ( i_start == -1 ) i_start = 1
if ( i_end == -1 ) then
if ( present(periodic) .and. periodic ) then
!!! If we haven't found and 'end', check if we need to wrap around
if ( end < start ) then
!!! Redo the loop from the start
do i=1,SIZE(arr)
if ( i_end == -1 ) then
if ( abs(arr(i) - end) < delta_coord ) i_end = i
else
exit
end if
end do
else
i_end = size(arr)
end if
else
i_end = size(arr)
end if
end if
END SUBROUTINE array_extents
END MODULE util
!**************************************************************
!****************************************************************
MODULE bt_subs
IMPLICIT NONE
CONTAINS
!***********************************************************************
SUBROUTINE new_out_file(outncid,wvcid,wvc2id,xlocid,ylocid,dayid,opreid,daynum,lat2d,lon2d)
!----------------------------------------------
!create the output netcdf file and prepare it to accept data
!--------------------------------------------------------
USE global_data
USE util
USE netcdf
IMPLICIT NONE
INTEGER, INTENT(INOUT) :: outncid,wvcid,wvc2id,xlocid,ylocid,dayid,opreid
!REAL, INTENT(IN) :: daynum
INTEGER, INTENT(IN) :: daynum
REAL,INTENT(IN),DIMENSION(:,:) :: lat2d,lon2d
INTEGER :: status,jdimid,idimid,gwvcdimid,latid,lonid
!
!create the file
!
!differentiate whether we are doing whole days or around storm peaks
if (peak) then
print *,'we are doing peaks here!'
if (mon<10) then
status = nf90_create(TRIM(diro)//"bt."//TRIM(int_to_string(year))//"0" &
//TRIM(int_to_string(mon))//"_"//TRIM(real_to_string(daynum))// &
".nc",nf90_clobber,outncid)
if (status /= NF90_NOERR) call handle_err(status)
else
status = nf90_create(TRIM(diro)//"bt."//TRIM(int_to_string(year)) &
//TRIM(int_to_string(mon))//"_"//TRIM(real_to_string(daynum))// &
".nc",nf90_clobber,outncid)
if (status /= NF90_NOERR) call handle_err(status)
end if
else
if (mon<10) then
status = nf90_create(TRIM(diro)//"bt."//TRIM(int_to_string(year))//"0" &
//TRIM(int_to_string(mon))//"_"//TRIM(int_to_string(INT(daynum)))// &
".nc",nf90_clobber,outncid)
if (status /= NF90_NOERR) call handle_err(status)
else
status = nf90_create(TRIM(diro)//"bt."//TRIM(int_to_string(year)) &
//TRIM(int_to_string(mon))//"_"//TRIM(int_to_string(INT(daynum)))// &
".nc",nf90_clobber,outncid)
if (status /= NF90_NOERR) call handle_err(status)
end if
end if
print *,'outfile=',TRIM(diro)//"bt."//TRIM(int_to_string(year))//"0" &
//TRIM(int_to_string(mon))//"_"//TRIM(int_to_string(INT(daynum)))// &
".nc"
!
!define dimensions
!
status = nf90_def_dim(outncid,"j_cross",dim_j,jdimid)
if (status /= NF90_NOERR) call handle_err(status)
status = nf90_def_dim(outncid,"i_cross",dim_i,idimid)
if (status /= NF90_NOERR) call handle_err(status)
status = nf90_def_dim(outncid,"gridcell_wvc",nf90_unlimited,gwvcdimid)
if (status /= NF90_NOERR) call handle_err(status)
!
!define the variable
!
status = nf90_def_var(outncid,"wv_cont",nf90_float,(/jdimid,idimid,gwvcdimid/),wvcid)
if (status /= NF90_NOERR) call handle_err(status)
!turning off apbl output
!status = nf90_def_var(outncid,"wv_cont_apbl",nf90_float,(/jdimid,idimid,gwvcdimid/),wvc2id)
!if (status /= NF90_NOERR) call handle_err(status)
status = nf90_def_var(outncid,"x_loc",nf90_int,(/gwvcdimid/),xlocid)
if (status /= NF90_NOERR) call handle_err(status)
status = nf90_def_var(outncid,"y_loc",nf90_int,(/gwvcdimid/),ylocid)
if (status /= NF90_NOERR) call handle_err(status)
status = nf90_def_var(outncid,"day",nf90_float,(/gwvcdimid/),dayid)
if (status /= NF90_NOERR) call handle_err(status)
status = nf90_def_var(outncid,"pre",nf90_float,(/gwvcdimid/),opreid)
if (status /= NF90_NOERR) call handle_err(status)
status = nf90_def_var(outncid,"latitcrs",nf90_float,(/jdimid,idimid/),latid)
if (status /= NF90_NOERR) call handle_err(status)
status = nf90_def_var(outncid,"longicrs",nf90_float,(/jdimid,idimid/),lonid)
if (status /= NF90_NOERR) call handle_err(status)
!
!define attributes
!
status = nf90_put_att(outncid,wvcid,"long_name","Water Vapor Contribution")
if (status /= NF90_NOERR) call handle_err(status)
status = nf90_put_att(outncid,wvcid,"units","proportion of precipitation")
if (status /= NF90_NOERR) call handle_err(status)
status = nf90_put_att(outncid,wvcid,"num_boundary_layers",bdy)
if (status /= NF90_NOERR) call handle_err(status)
status = nf90_put_att(outncid,wvcid,"parcels_per_grid_point",nparcels)
if (status /= NF90_NOERR) call handle_err(status)
!turn off apbl output
!status = nf90_put_att(outncid,wvc2id,"long_name","Water Vapor Contribution above PBL")
!if (status /= NF90_NOERR) call handle_err(status)
!status = nf90_put_att(outncid,wvc2id,"units","proportion of precipitation")
!if (status /= NF90_NOERR) call handle_err(status)
!status = nf90_put_att(outncid,wvc2id,"num_boundary_layers",bdy)
!if (status /= NF90_NOERR) call handle_err(status)
!status = nf90_put_att(outncid,wvc2id,"parcels_per_grid_point",nparcels)
!if (status /= NF90_NOERR) call handle_err(status)
status = nf90_put_att(outncid,xlocid,"long_name","x index location of precipitation (from 0)")
if (status /= NF90_NOERR) call handle_err(status)
status = nf90_put_att(outncid,ylocid,"long_name","y index location of precipitation (from 0)")
if (status /= NF90_NOERR) call handle_err(status)
status = nf90_put_att(outncid,dayid,"long_name","days since "// &
TRIM(int_to_string(sday))//"/"//TRIM(int_to_string(smon))//"/"// &
TRIM(int_to_string(syear)))
if (status /= NF90_NOERR) call handle_err(status)
status = nf90_put_att(outncid,opreid,"long_name","precipitation")
if (status /= NF90_NOERR) call handle_err(status)
status = nf90_put_att(outncid,opreid,"units","mm")
if (status /= NF90_NOERR) call handle_err(status)
status = nf90_put_att(outncid,latid,"long_name","LATITUDE (SOUTH NEGATIVE)")
if (status /= NF90_NOERR) call handle_err(status)
status = nf90_put_att(outncid,latid,"units","degrees")
if (status /= NF90_NOERR) call handle_err(status)
status = nf90_put_att(outncid,lonid,"long_name","LONGITUDE (WEST NEGATIVE)")
if (status /= NF90_NOERR) call handle_err(status)
status = nf90_put_att(outncid,lonid,"units","degrees")
if (status /= NF90_NOERR) call handle_err(status)
!
!leave define mode
!
status = nf90_enddef(outncid)
status = nf90_put_var(outncid,latid,lat2d,start=(/1,1/),count=(/dim_j,dim_i/))
if(status /= nf90_NoErr) call handle_err(status)
status = nf90_put_var(outncid,lonid,lon2d,start=(/1,1/),count=(/dim_j,dim_i/))
if(status /= nf90_NoErr) call handle_err(status)
END SUBROUTINE new_out_file
!***********************************************************************
!***********************************************************************
REAL FUNCTION lin_interp(var,fac)
!---------------------------------------
!linearly interpolate between the values of var
!fac is the proporational distance from the first value
!------------------------------------------
IMPLICIT NONE
REAL, INTENT(IN), DIMENSION(2) :: var
REAL, INTENT(IN) :: fac
lin_interp = var(1)*(1-fac) + var(2)*fac
END FUNCTION lin_interp
!***********************************************************************
FUNCTION lin_interp2D(var,fac)
!---------------------------------------
!linearly interpolate between the values of var (last dimension must have size 2)
!fac is the proporational distance from the first value
!------------------------------------------
IMPLICIT NONE
REAL, INTENT(IN), DIMENSION(:,:,:) :: var
REAL, INTENT(IN) :: fac
REAL, DIMENSION(SIZE(var,1),SIZE(var,2)):: lin_interp2D
lin_interp2D = var(:,:,1)*(1-fac) + var(:,:,2)*fac
END FUNCTION lin_interp2D
!***********************************************************************
FUNCTION lin_interp3D(var,fac)
!---------------------------------------
!linearly interpolate between the values of var (last dimension must have size 2)
!fac is the proporational distance from the first value
!------------------------------------------
IMPLICIT NONE
REAL, INTENT(IN), DIMENSION(:,:,:,:) :: var
REAL, INTENT(IN) :: fac
REAL, DIMENSION(SIZE(var,1),SIZE(var,2),SIZE(var,3)) :: lin_interp3D
lin_interp3D = var(:,:,:,1)*(1-fac) + var(:,:,:,2)*fac
END FUNCTION lin_interp3D
!***********************************************************************
SUBROUTINE parcel_release_time(precip,npar,par_release)
!---------------------------------------------------
! here we calculate the times of day to release our parcels
! based on a random precipitation weighted sampling
!-----------------------------------------------------
USE global_data
IMPLICIT NONE
REAL, DIMENSION(:), INTENT(IN) :: precip
INTEGER, DIMENSION(:), INTENT(OUT) :: par_release
INTEGER, INTENT(IN) :: npar
REAL, DIMENSION(SIZE(precip)*indatatsteps) :: cumm_precip
REAL, DIMENSION(npar) :: rand_nums
INTEGER :: tt,rr,ss,rec
par_release = 0
call RANDOM_NUMBER(rand_nums)
cumm_precip = 0.
rec = 0
do tt = 1,SIZE(precip)
do ss = 1,indatatsteps
rec = rec + 1
if (rec==1) then
cumm_precip(1) = precip(1)/indatatsteps
else
cumm_precip(rec) = cumm_precip(rec-1) + precip(tt)/indatatsteps
end if
end do
end do
cumm_precip = cumm_precip/cumm_precip(SIZE(cumm_precip))
do rr = 1,npar
do tt = 1,SIZE(cumm_precip)
if (cumm_precip(tt)>rand_nums(rr)) then
par_release(tt) = par_release(tt) + 1
EXIT
end if
end do
end do
END SUBROUTINE parcel_release_time
!***********************************************************************
SUBROUTINE parcel_release_height(pw,par_lev)
!----------------------------------------------
! calculate the height to release the parcel from
! based on precipitable water weighted random sampling
!-----------------------------------------------
IMPLICIT NONE
REAL, INTENT(IN), DIMENSION(:) :: pw ! This is the precipitable water, accumulated from the ground up, in the column at that point in time.
INTEGER, INTENT(OUT) :: par_lev
REAL :: rand_num
INTEGER :: kk
call RANDOM_NUMBER(rand_num)
! Take random number as a random proportion of the total pw in the column at that time; pw(1) is the pw at the top of the atm column, which represents the accumulated pw over the column below it (same as TPW).
rand_num = rand_num*pw(1)
do kk = SIZE(pw),1,-1 ! kk is then between 29 and 1
if (pw(kk)>rand_num) then
par_lev = kk
EXIT
end if
end do
!print *,"release height,pw, ",par_lev,pw
! For testing purposes only: take random number as a purely random model level, not weighted by pw.
!rand_num = 1 + FLOOR(size(pw)*rand_num)
!par_lev = rand_num
END SUBROUTINE parcel_release_height
!***********************************************************************
SUBROUTINE lin_interp_inMM5tsteps(var)
!------------------------------------------
!linearly interpolate through time inside MM5 time steps
!----------------------------------------------
USE global_data
IMPLICIT NONE
REAL, INTENT(INOUT), DIMENSION(:,:,:,:) :: var
INTEGER :: i
do i = 1,indatatsteps-1
var(:,:,:,i+1::indatatsteps) = (1-(i*1./indatatsteps))*var(:,:,:,1:datadaysteps+1-indatatsteps:indatatsteps) &
& + (i*1./indatatsteps)*var(:,:,:,1+indatatsteps::indatatsteps)
end do
END SUBROUTINE lin_interp_inMM5tsteps
!***********************************************************************
SUBROUTINE calc_pw(mix,pres,surf_pres,ptop,pw)
!------------------------------------------
! calculate the precipitable water from input data fields
! only for the day of current interest
! save as accumulated field from the ground up (4D field)
! This is used to calc parcel initial height.
!-------------------------------------------------
USE global_data
IMPLICIT NONE
REAL, INTENT(IN), DIMENSION(:,:,:,:) :: mix,pres
REAL, INTENT(IN), DIMENSION(:,:,:) :: surf_pres
REAL, INTENT(OUT), DIMENSION(:,:,:,:) :: pw
REAL, DIMENSION(dim_j,dim_i,dim_k,SIZE(mix,4)) :: dp
INTEGER :: k
REAL, INTENT(IN) :: ptop
!
! calculate the change in pressure (Pa) represented by each point
!
!for highest level
dp(:,:,1,:) = SUM(pres(:,:,:2,:),3)/2. - ptop
!need to account for posibility that pressure levels go below the ground
!for the middle levels
do k = 2,dim_k-1
where (pres(:,:,k+1,:) <= surf_pres(:,:,:))
dp(:,:,k,:) = (pres(:,:,k+1,:) - pres(:,:,k-1,:)) /2. !dp(:,:,k,:) = SUM(pres(:,:,k-1:k+1:2,:),3)/2.
elsewhere (pres(:,:,k,:) <= surf_pres(:,:,:))
!for the lowest level above surface
dp(:,:,k,:) = surf_pres(:,:,:) - (pres(:,:,k,:) + pres(:,:,k-1,:))/2.
elsewhere
dp(:,:,k,:) = 0.
end where
end do
!for the lowest level
where (pres(:,:,dim_k,:) <= surf_pres(:,:,:))
dp(:,:,dim_k,:) = surf_pres(:,:,:) - SUM(pres(:,:,dim_k-1:,:),3)/2.
elsewhere
dp(:,:,dim_k,:) = 0.
end where
!mass in mm
pw(:,:,:,::indatatsteps) = dp*mix/g
!interpolate inside input data time steps
call lin_interp_inMM5tsteps(pw)
!accumulate from the bottom up. The precipitable water is then the total moisture in the column below it.
do k = dim_k-1,1,-1 ! i.e. from level 28 to 1
!pw(:,:,k,2:) = pw(:,:,k+1,2:) + pw(:,:,k,2:)! THE PW IS ACCUMULATED FROM THE SECOND TS ON, SO THE FIRST 10MIN IS NOT ACCUMULATED. WHY?? This only matters if tt=1. Here I change it to do all timesteps.
pw(:,:,k,:) = pw(:,:,k+1,:) + pw(:,:,k,:)
end do
! print *,shape(pw)
! print *,"pw ",pw(1,1,dim_k,:),surf_pres(1,1,1),ptop
END SUBROUTINE calc_pw
!***********************************************************************
SUBROUTINE calc_tpw(mix,pres,surf_pres,ptop,tpw)
!------------------------------------------
! calculate the total precipitable water from input data fields
! at every level and time (3D field)
!-------------------------------------------------
USE global_data
IMPLICIT NONE
REAL, INTENT(IN), DIMENSION(:,:,:,:) :: mix,pres
REAL, INTENT(IN), DIMENSION(:,:,:) :: surf_pres
REAL, INTENT(OUT), DIMENSION(:,:,:) :: tpw
REAL, DIMENSION(dim_j,dim_i,dim_k,datatotsteps) :: dp
INTEGER :: j,i,k,t
REAL, INTENT(IN) :: ptop
!
! calculate the change in pressure (Pa) represented by each point
!
!for highest level
dp(:,:,1,:) = SUM(pres(:,:,:2,:),3)/2. - ptop
!for the middle levels
!do k = 2,dim_k-1
! dp(:,:,k,:) = (pres(:,:,k+1,:) - pres(:,:,k-1,:)) /2.
! ! equiv to (p2+p3)/2 - (p1+p2)/2 = (p3-p1)/2
!end do
!for the lowest level
!dp(:,:,dim_k,:) = surf_pres(:,:,:) - SUM(pres(:,:,dim_k-1:,:),3)/2.
!need to account for posibility that pressure levels go below the ground
!for the middle levels
do k = 2,dim_k-1
where (pres(:,:,k+1,:) <= surf_pres(:,:,:))
dp(:,:,k,:) = (pres(:,:,k+1,:) - pres(:,:,k-1,:)) /2. !dp(:,:,k,:) = SUM(pres(:,:,k-1:k+1:2,:),3)/2.
elsewhere (pres(:,:,k,:) <= surf_pres(:,:,:))
!for the lowest level above surface
dp(:,:,k,:) = surf_pres(:,:,:) - (pres(:,:,k,:) + pres(:,:,k-1,:))/2.
elsewhere
dp(:,:,k,:) = 0.
end where
end do
!for the lowest level also accounting for possibility of being below ground
where (pres(:,:,dim_k,:) <= surf_pres(:,:,:))
dp(:,:,dim_k,:) = surf_pres(:,:,:) - SUM(pres(:,:,dim_k-1:,:),3)/2.
elsewhere
dp(:,:,dim_k,:) = 0.
end where
!mass in mmtpw
do j = 1,dim_j
do i = 1,dim_i
do t = 1,datatotsteps
tpw(j,i,t) = SUM(dp(j,i,:,t)*mix(j,i,:,t)/g)
end do
end do
end do
END SUBROUTINE calc_tpw
!***********************************************************************
SUBROUTINE calc_tpw_pbl(mix,pres,surf_pres,tpw,pbl_lev)
!------------------------------------------
! calculate the total precipitable water in the pbl from MM5 fields
! at every level and time