#!/usr/bin/perl -w
# vim:ts=4:sw=4:sts=4:si

#Script for aa kjore flextra trajektoriemodell for $numdays dager, bakover

use Getopt::Long;
use File::Copy qw(mv);
$numdays=7;
$date=0;
$plot_only=0;
$run_ftp = 0;
$work_dir = '/storage/home/flextra/flextra5.0/';
#$metfile_size = 62958720;
$metfile_size = 64298808;


my $ret = GetOptions(
		     "--days=f",\$numdays,
		     "--date=f",\$date,
		     "--plot!",\$plot_only,
		     "--run_ftp!",\$run_ftp,
		     "--work_dir=s",\$work_dir,
		     "--help",\&usage,
		     "<>",\&usage);

@daysinmonth = (31,28,31,30,31,30,31,31,30,31,30,31);

#Get current day
my ($sec,$min,$hour,$day,$mon,$year,$weekday,$yday,$isdst) = localtime(time);
if ($date > 0) {
#Run flextra for $date
    $day = $date % 100;
    $date = ($date - $day)/100;
    $mon = $date % 100;
    $year = ($date-$mon)/100;
    $year_2d = $year % 100;
} else {
    $mon++;
    $year_2d = $year;
    if ($year_2d >= 100) {$year_2d=$year_2d-100;}
    $year=$year+1900;
    &subtract_days($year,$mon,$day,3);
}


#Lese fil med stasjoner og koordinater
#open(FP,'/storage/nilu/Atmos/Projects/Active_projects/o7726_SOGE/traj/emma/list.txt');   
open(FP,'station_list_part1.txt');
#open(FP,'/nilu/reg/soge/traj/emma/station_list.txt');
@stations = <FP>;
close(FP);
foreach $l (@stations) {
    chop $l;
    $l =~ s/\cM//g;
}




################################################################################
#   Kjøre modellen
################################################################################


if ($plot_only == 0) {
    printf STDERR "Running flextra for %d/%d/%d:\n",$day,$mon,$year;


#Generere AVAILABLE-fil for de siste $numdays+1 dager
#    print STDERR "...Generating AVAILABLE\n";
#    $met_errors = &print_available($year,$mon,$day,$numdays);


#Sjekke om metdata er paa plass
    $met_errors = check_metdata($year,$mon,$day,$numdays);
    if ($met_errors > 0) {
	print STDERR "Error in met data. Aborting....\n";
	exit(1);
    }

#Generere COMMAND-fil
    print STDERR "...Generating COMMAND\n";
    &print_command($year,$mon,$day,$numdays);

#Generere STARTPOINTS-fil
    printf STDERR "...Generating STARTPOINTS\n";
    &print_startpoints($year,$mon,$day,$numdays,$height);
    
    printf STDERR "...Running flextra for %d days\n",$numdays;
    $command = $work_dir . "/FLEXTRA_ECMWF";
    my $ret_flextra = `$command`;
    print $ret_flextra;
    print "\n";
    printf STDERR "...Done flextra\n";

#flytte filer fra output til modelldata/<aar>/<stasjon>/:

    foreach $i (2 .. $#stations) {
	($code,$station_name,$country,$region,$progam_acr,$long,$lat,$area,$unit,@heights) = split(/ +/,$stations[$i]); # Split on blanks
    
	foreach $height (@heights) {

	    $filename = sprintf "TI_%s_%d_%04d%02d%02d",$code,$height,$year,$mon,$day;
	    $dest_dir = sprintf "/viper/nadir/traj/modelldata/%04d/%s/",$year,$code;

		############################ <fiksing> ############################
	    my $src_file = "output/$filename";
	    my $dst_file = "$dest_dir/$filename";

	    if ( -f "$src_file" ) {
			chdir("/storage/nilu/Atmos/Projects/Active_projects/o7726_SOGE/traj/emma/version5/");
			my @stat = stat($src_file);
			my $perm = $stat[2] & 07777;
			$perm |= 055;
			chmod $perm,$src_file;
			mv($src_file,$dst_file);
	    } else {
#			print STDERR "Source file $src_file not found\n";
		}
		############################ </fiksing> ###########################

	    #`cd /nilu/reg/soge/traj/emma; chmod a+rx output/$filename; mv output/$filename $dest_dir`;
	}
    }
}




################################################################################
#   Plotte data
################################################################################


if ($plot_only) {

#Generere plott
    printf STDERR "...Plotting trajectories\n";

    open(FP,"|/data/rsi/envi_3.5/idl_5.5/bin/idl");
#open(FP,"|/data/idl-5.3/idl_5.3/bin/idl");
#open(FP,"|idl");

    foreach $i (2 .. $#stations) {
	($code,$station_name,$country,$region,$progam_acr,$long,$lat,$area,$unit,@heights) = split(/ +/,$stations[$i]); # Split on blanks
	
	$zlevs="[";
	foreach $i (0 .. $#heights-1) {
	    $zlevs = $zlevs . $heights[$i] . ",";
	}
	$zlevs = $zlevs . $heights[$#heights] . "]";

	printf FP "pl_soge,%d,%d,%d,code='%s',zlevs=%s,name='%s',ztype=%d,area=%d,/png\n",$year,$mon,$day,$code,$zlevs,$station_name,$unit,$area;

    }
    print FP "exit\n";
    close(FP);
    printf STDERR "...Done plotting  AMF!\n";
}    
 



################################################################################
#   Kjøre ftp tilwebben
################################################################################


if ($run_ftp) {
    
#Fikse skriveaksess på pngfiler og ftp til tarantula
    $y_tmp = $year;
    $m_tmp = $mon;
    $d_tmp = $day;
	
#ftp til tarantula
    open(FP,">ftp.sh");
    
    print FP "#!/bin/sh\n";
    print FP "ftp -i sky.nilu.no << /EOF\n";
    print FP "bin\n";
	
    $fil = sprintf "%4d%02d%02d", $year,$mon,$day;
    $fil = $fil . "*.png";
    `cd /nilu/reg/soge/traj/emma/png; chmod g+w $year/*/$fil`;
    
    foreach $i (2 .. $#stations) {
	($code,$station_name,$country,$region,$progam_acr,$long,$lat,$area,$unit,@heights) = split(/ +/,$stations[$i]); # Split on blanks
#	print FP "cd /WWWroot/Trajectories/files_emma/png/$year/$code\n";
	print FP "cd /wwwroot/Trajectories/files/png/$year/$code\n";
	print FP "lcd /nilu/reg/soge/traj/emma/png/$year/$code\n";
	printf FP "mput %s\n",$fil;
    }
    
    subtract_days($year,$mon,$day,1);
    $fil = sprintf "%4d%02d%02d", $year,$mon,$day;
    $fil = $fil . "*.png";
    `cd /nilu/reg/soge/traj/emma/png; chmod g+w $year/*/$fil`;
    
    foreach $i (2 .. $#stations) {
        ($code,$station_name,$country,$region,$progam_acr,$long,$lat,$area,$unit,@heights) = split(/ +/,$stations[$i]); # Split on blanks
#	print FP "cd /WWWroot/Trajectories/files_emma/png/$year/$code\n";
	print FP "cd /wwwroot/Trajectories/files/png/$year/$code\n";
	print FP "lcd /nilu/reg/soge/traj/emma/png/$year/$code\n";
	printf FP "mput %s\n",$fil;
    }
    
    print FP "bye\n";
    print FP "/EOF\n";
    close(FP);
	
    $year = $y_tmp;
    $mon = $m_tmp;
    $day = $d_tmp;
    
    printf STDERR "...Running ftp\n";
	
    `sh ftp.sh`;
    `zip_filer.pl`;  
}




################################################################################
#   Subrutiner
################################################################################

sub check_metdata {
    $year = shift @_;
    $mon = shift @_;
    $day = shift @_;
    $numdays = shift @_;

    $error_files = 0;

    subtract_days($year,$mon,$day,$numdays+1);

    for ($i=0;$i<$numdays+1;$i++) {
	add_days($year,$mon,$day,1);
	foreach $j (0,3,6,9,12,15,18,21) { 
	    $fil = sprintf "EN%02d%02d%02d%02d",$year%100,$mon,$day,$j;

#	    $metfil = "/xnilu_wrk/flex_wrk/WIND_FIELDS/all_fields_oper_ecmwf/" . $fil;
	   $metfil = "/xnilu_wrk/projects/FLEXPART/flex_wrk/WIND_FIELDS/all_fields_oper_ecmwf/" . $fil;
	    if ( -e $metfil) {

		$size = -s $metfil;
		if ($size != $metfile_size) {
#		    print "$metfil: wrong size: $size\n";
#		    $error_files++;
		}
	    } else {
		print "$metfil: does not exist\n";
		$error_files++;
	    }
	}
    }

    $error_files;
}


sub print_command{
    $year = shift @_;
    $mon = shift @_;
    $day = shift @_;
    $numdays = shift @_;

#    open(FP,">/nilu2/home/flextra/flextra3.3/options/COMMAND") || die "Unable to open COMMAND\n";
    open(FP,">" . $work_dir . "/options/COMMAND") || die "Unable to open COMMAND\n";

    print FP "********************************************************************************\n";
    print FP "*                                                                              *\n";
    print FP "*    Input file for the trajectory model FLEXTRA: Please select your options   *\n";
    print FP "*                                                                              *\n";
    print FP "********************************************************************************\n";
    print FP "\n";
    print FP "'Test0001'                          LABEL FOR THE MODEL RUN\n";
    print FP "-1                                  DIRECTION (1 FORWARD, -1 BACKWARD TRAJECT.)\n";
    printf FP "%3d0000           HHHMISS           LENGTH OF AN INDIVIDUAL TRAJECTORY\n",$numdays * 24;
    subtract_days($year,$mon,$day,1);
    printf FP "%4d%02d%02d 180000   YYYYMMDD HHMISS   BEGINNING DATE\n",$year,$mon,$day;
    add_days($year,$mon,$day,1);
    printf FP "%4d%02d%02d 120000   YYYYMMDD HHMISS   ENDING DATE\n",$year,$mon,$day;
    print FP " 030000           HHHMISS           TIME INTERVAL BETWEEN STARTING TIMES OF TRAJECTORIES\n";
    print FP "1 10800           i  SSSSS          i>0: INTERPOLATED OUTPUT OF TRAJECTORY EVERY SSSSS SECONDS\n";
    print FP "0  0.5 2.0 0.08 0.08 0.20           NUMBER, DISTANCE (GRID UNITS), TIME CONSTANT (WIND FIELD INTERVAL UNITS) AND INTERPOLATION ERRORS (IN U, V AND W) OF UNCERTAINTY TRAJECT.\n";
    print FP "1                 INTERPOLATION     1 = IDEAL INTERPOLATION   >1 = LINEAR INTERPOLATION\n";
    print FP "5.0               CFL               TIMESTEP CRITERION HORIZONTAL AND VERTICAL\n";
    print FP "5.0               CFLT              TIMESTEP CRITERION TIME GAP BETWEEN INPUT WIND FIELDS\n";
    print FP "1                 MODE              1 NORMAL MODE, 2 CET MODE, 3 FLIGHT MODE\n";
    print FP "==================================================================================\n";
    print FP "1. Comment to identify the current model run\n";
    print FP "\n";
    print FP "2. Direction of trajectories (1 means forward trajectories, -1 backward)\n";
    print FP "\n";
    print FP "3. Temporal lengths of the trajectories in the format HHHMISS, where HHH is \n";
    print FP "   HOUR, MI is MINUTE and SS is SECOND\n";
    print FP "\n";
    print FP "4. Beginning date and time of trajectory calculation. Must be given in format\n";
    print FP "   YYYYMMDD HHMISS, where YYYY is YEAR, MM is MONTH, DD is DAY, HH is HOUR,\n";
    print FP "   MI is MINUTE and SS is SECOND. All times are in UTC.\n";
    print FP "\n";
    print FP "5. Ending date and time of trajectory calculation. Same format as 4.\n";
    print FP "\n";
    print FP "6. Time interval between two trajectory calculations. Same format as 3.\n";
    print FP "\n";
    print FP "7. Options for the trajectory output:     \n";
    print FP "   0 = original data in irregular time intervals\n";
    print FP "   1 = constant time intervals, interpolated output every SSSSS seconds\n";
    print FP "   2 = 0 plus 1\n";
    print FP "\n";
    print FP "8. Six parameters have to be inputted. The first is the number of\n";
    print FP "   uncertainty trajectories. They are starting in a distance from\n";
    print FP "   the starting point of the reference trajectory as given by the\n";
    print FP "   second parameter (in grid units).\n";
    print FP "   Additionally, random errors may be added at each time step of the\n";
    print FP "   trajectory calculation. Using a Langevin equation, they are relaxed \n";
    print FP "   with a time constant (in units of the wind field interval, third\n";
    print FP "   parameter) specified by the user. These random errors are\n";
    print FP "   thought to reflect typical wind errors caused, for instance, by\n";
    print FP "   interpolation. The magnitude of these errors (in relative units,\n";
    print FP "   relative to the wind velocity) must be specified by the user for\n";
    print FP "   the three wind components u, v and w (last three parameters).\n";
    print FP "\n";
    print FP "9. Kind of interpolation\n";
    print FP "   1  - horizontal interpolation bicubic\n";
    print FP "        vertical interpolation polynomial\n";
    print FP "        temporal interpolation linear\n";
    print FP "   >1 - horizontal interpolation bilinear         \n";
    print FP "        vertical interpolation linear\n";
    print FP "        temporal interpolation linear\n";
    print FP " \n";
    print FP "10.cfl criterion horizontal/vertical\n";
    print FP "   factor by which time step must be shorter than that determined\n";
    print FP "   from the CFL criterion, i.e.\n";
    print FP "\n";
    print FP "                delta_t1=delta x/u/cfl\n";
    print FP "                delta_t2=delta y/v/cfl\n";
    print FP "                delta_t3=delta z/w/cfl\n";
    print FP "\n";
    print FP "                delta_t(space) = min(delta_t1,delta_t2,delta_t3)\n";
    print FP "\n";
    print FP "11. cfl criterion time\n";
    print FP "    factor by wich time is shorter than time interval of the wind\n";
    print FP "    fields\n";
    print FP "\n";
    print FP "                delta_t(time) = delta_T(input wind)/cflt\n";
    print FP "\n";
    print FP "     The time step used for trajectory calculation is the minimum of\n";
    print FP "     delta_t(space) and delta_t(time)\n";
    print FP "\n";
    print FP "     cfl and cflt must not be less than 1!\n";
    print FP "\n";
    print FP "12.  1 NORMAL mode -> read file STARTPOINTS and calculate a time\n";
    print FP "     series of trajectories starting all from the same starting\n";
    print FP "     points\n";
    print FP "     2 CET mode -> read file STARTCET and calculate trajectories\n";
    print FP "     starting uniformly spaced from a user-defined domain\n";
    print FP "     (for a single starting time)\n";
    print FP "     3 FLIGHT mode -> read file STARTFLIGHT and calculate\n";
    print FP "     trajectories starting neither uniformly spaced nor with\n";
    print FP "     constant time intervals (as needed, for instance, to start\n";
    print FP "     trajectories along an aircraft leg)\n";

    close(FP);
}

sub print_startpoints{
    $year = shift @_;
    $mon = shift @_;
    $day = shift @_;
    $numdays = shift @_;
    $height  = shift @_;

#    open(FP,">/nilu2/home/flextra/flextra3.3/options/STARTPOINTS") || die "Unable to open STARTPOINTS\n";
    open(FP,">" . $work_dir . "/options/STARTPOINTS") || die "Unable to open STARTPOINTS\n";

    print FP "**********************************************************************\n";
    print FP "*                                                                    *\n";
    print FP "*                 TRAJECTORY MODEL                                   *\n";
    print FP "*                 DEFINITION OF STARTING/ENDING POINTS               *\n";
    print FP "*                                                                    *\n";
    print FP "*  The first 7 characters of the comment are also used as filenames. *\n";
    print FP "*  Therefore, they cannot be blank and they must be different for    *\n";
    print FP "*  each starting point.                                              *\n";
    print FP "*                                                                    *\n";
    print FP "*  Kind of trajectory: 1 = 3 dimensional                             *\n";
    print FP "*                      2 = on model layers                           *\n";
    print FP "*                      3 = mixing layer                              *\n";
    print FP "*                      4 = isobaric                                  *\n";
    print FP "*                      5 = isentropic                                *\n";
    print FP "*                                                                    *\n";
    print FP "**********************************************************************\n";
    print FP "*                                                                    *\n";
    print FP "*  Unit of z coordinate: 1 = Meters above sea level                  *\n";
    print FP "*                        2 = Meters above ground                     *\n";
    print FP "*                        3 = Hectopascal                             *\n";
    print FP "*                                                                    *\n";
    print FP "*  For mixing layer trajectories (kind 3), the z coordinate must be  *\n";
    print FP "*  given in m.a.g.l. (option 2)                                      *\n";
    print FP "*                                                                    *\n";
    print FP "**********************************************************************\n";
    print FP "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n";

    foreach $i (2 .. $#stations) {
	($code,$station_name,$country,$region,$progam_acr,$long,$lat,$area,$unit,@heights) = split(/ +/,$stations[$i]); # Split on blanks

	foreach $height (@heights) {
	    printf FP " %6.2f                    F Longitude [DEG]\n",$long;
	    printf FP " %6.2f                    F Latitude [DEG]\n",$lat;
	    print  FP "1                          I Kind of trajectory (see file header)\n";
	    printf FP "%d                          I Unit of z coordinate\n",$unit;
	    printf FP " %6.1f                    F z-coordinate (see file header)\n",$height;
	    printf FP "'%s_%d_%04d%02d%02d'     A Name of starting point\n",$code,$height,$year,$mon,$day;
	    print  FP "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n";
	}
    }
    close(FP);
}

sub add_days {
    $year = shift @_;
    $mon = shift @_;
    $day = shift @_;
    $daystoadd = shift @_;

    if ($year % 4 == 0) {
	$daysinmonth[1] = 29;
    } else {
	$daysinmonth[1] = 28;
    }

    $day = $day + $daystoadd;

    while ($day > $daysinmonth[$mon-1]) {
	$day = $day - $daysinmonth[$mon-1];
	$mon++;
	
	if ($mon == 13) {
	    $mon = 1;
	    $year++;
	    if ($year % 4 == 0) {
		$daysinmonth[1] = 29;
	    } else {
		$daysinmonth[1] = 28;
	    }
	}
    }
}

sub subtract_days {
    $year = shift @_;
    $mon = shift @_;
    $day = shift @_;
    $daystosubtract = shift @_;

    if ($year % 4 == 0) {
	$daysinmonth[1] = 29;
    } else {
	$daysinmonth[1] = 28;
    }

    $day = $day - $daystosubtract;

    while ($day < 1) {
	$day = $day + $daysinmonth[$mon-2];
	$mon--;
	
	if ($mon == 0) {
	    $mon = 12;
	    $year--;
	    if ($year % 4 == 0) {
		$daysinmonth[1] = 29;
	    } else {
		$daysinmonth[1] = 28;
	    }
	}
    }
}


sub usage {
    print STDERR "$0 runs flextra trajectory model\n";
    print STDERR "Options:\n";
    printf STDERR "-days: Number of days for the trajectory. Default is %d days\n",$numdays;
    print STDERR "-date: Run model for a given date, eg. 20020130. Default is two days back\n";
    print STDERR "-plot: Run the IDL routine for plotting the trajectories\n";
    print STDERR "-run_ftp: Run the ftp command for transfering the png files to tarantula\n";
    print STDERR "-help: Prints this message\n";
#    print STDERR "\n";
    exit();
}