Commit 75b6f144674ce9aedbecf6ee72d75d8123d672c3

Authored by dumoda01
1 parent 5949b5b3
Exists in master and in 1 other branch snow

Script pour la génération d'un ensemble de forçages NARR. Work-in-progress.

Showing 1 changed file with 189 additions and 0 deletions   Show diff stats
scripts/matlab/narr2gotm.m 0 → 100644
... ... @@ -0,0 +1,189 @@
  1 +function narr2gotm(name,year)
  2 +% cf. ftp://ftp.cdc.noaa.gov/Datasets/NARR/monolevel/
  3 +% cf. GOTM
  4 +% Uses SNCTools
  5 +
  6 +% Dany DUMONT
  7 +% 14/05/2008: Creation
  8 +% 23/02/2011: Put year as an argument
  9 +% 23/11/2011: Generalize the method to create many members for different
  10 +% regions
  11 +% ______________________________________________________________________
  12 +
  13 +yrstr = num2str(year)
  14 +
  15 +if mod(year,400) == 0
  16 + nod = 366;
  17 +elseif mod(year,100) == 0
  18 + nod = 365;
  19 +elseif mod(year,4) == 0
  20 + nod = 366;
  21 +else
  22 + nod = 365;
  23 +end
  24 +
  25 +% Time scale given in Matlab serial date number
  26 +init_date = datenum([year 01 01 00 00 00]);
  27 +time = [init_date:1/24:init_date+nod-1/24]';
  28 +timevec = datevec(time);
  29 +tlength = length(time);
  30 +
  31 +%% Horizontal coordinates
  32 +
  33 +ncfile = ['latlon.narr.',name,'.nc'];
  34 +lon = nc_varget(ncfile,'LON');
  35 +lat = nc_varget(ncfile,'LAT');
  36 +
  37 +%% Choosing members
  38 +
  39 +% Indices of chosen locations
  40 +i_loc = [5 5 6 6 6 7 7 7 8 8 8 9];
  41 +j_loc = [4 5 3 4 5 3 4 5 3 4 5 3];
  42 +nloc = length(i_loc);
  43 +
  44 +figure
  45 +gray = [0.55 0.55 0.55];
  46 +ax = worldmap([67 75],[-135 -115]);
  47 +tightmap
  48 +gridm('-')
  49 +mlabel('on')
  50 +plabel('on')
  51 +geoshow('landareas.shp', 'FaceColor',gray)
  52 +setm(gca,'FontName','Times','FontSize',12,'FontWeight','bold', ...
  53 + 'MLabelLocation',[-135 -125], ...
  54 + 'PLabelLocation',[70 75], ...
  55 + 'PLineLocation',[70 75], ...
  56 + 'PLineLimit', [-135 -115], ...
  57 + 'MLineLocation',[-135 -130 -120], ...
  58 + 'MLineLimit', [67 75], ...
  59 + 'FLineWidth',1, ...
  60 + 'GColor', [0.4 0.4 0.4])
  61 +geoshow('worldlakes.shp', 'FaceColor', 'white');
  62 +plotm(lat,lon,'ok')
  63 +
  64 +return
  65 +
  66 +
  67 +%% Momentum flux
  68 +
  69 +ncfile = ['uwnd.narr.',name,'.',yrstr,'.nc'];
  70 +uwnd = nc_varget(ncfile,'UWND');
  71 +
  72 +ncfile = ['vwnd.narr.',name,'.',yrstr,'.nc'];
  73 +vwnd = nc_varget(ncfile,'VWND');
  74 +
  75 +% interpolation des donn???es ??? chaque heure
  76 +uwnd = interp(uwnd,3);
  77 +vwnd = interp(vwnd,3);
  78 +
  79 +% Wind stress calculation
  80 +% cf. Stewart, R.H. (2002) Introduction to Physical Oceanography
  81 +rho = 1.3; % density of air (kg m^-3)
  82 +U = sqrt(uwnd.^2 + vwnd.^2); % velocity amplitude (m s^-1)
  83 +C10 = nan .* ones(length(U),1);
  84 +for i = 1:length(U)
  85 + if (U(i) >= 3 && U(i) < 6)
  86 + C10(i) = (0.29 + 3.1./U(i) + 7.7./U(i).^2).*1e-3;
  87 + elseif U(i) >= 6
  88 + C10(i) = (0.60 + 0.070.*U(i)).*1e-3; % Yelland and Taylor (1996)
  89 + else
  90 + C10(i) = 2.18e-3;
  91 + end
  92 +end
  93 + % C10 = (0.75 + 0.065.*U).*1e-3; % Garratt (1977)
  94 +dir = atan2(vwnd,uwnd); % wind direction (rad)
  95 +tau = rho.*C10.*U.*U; % wind stress amplitude (N m^-2)
  96 +taux = tau.*cos(dir); % zonal wind stress (N m^-2)
  97 +tauy = tau.*sin(dir); % meridional wind stress (N m^-2)
  98 +% figure; plot(U,C10,'ob')
  99 +% figure; plot(U,tau,'ok')
  100 +clear U dir tau C10 rho
  101 +
  102 +% Creation of the file 'narr_momentumflux.dat'
  103 +% file with tx and ty given in N/m^2
  104 +momentumflux = nan .* ones(tlength,8);
  105 +momentumflux(:,1:6) = timevec;
  106 +momentumflux(:,7) = taux;
  107 +momentumflux(:,8) = tauy;
  108 +
  109 +fid = fopen(['narr_hourly_momentumflux_',yrstr,'.dat'], 'wt');
  110 +fprintf(fid, '%4.0f-%02.0f-%02.0f %02.0f:%02.0f:%02.0f\t%f\t%f\n', ...
  111 + momentumflux');
  112 +fclose(fid);
  113 +
  114 +
  115 +% ----------------------------------------------------------------------
  116 +% Incoming shortwave radiative flux and surface heat flux
  117 +% from NARR 3-hourly reanalysis
  118 +% ----------------------------------------------------------------------
  119 +ncfile = ['dswrf_',yrstr,'_narr_amdgulf.nc'];
  120 +dswrf = nc_varget(ncfile,'DSWRF'); % W/m^2
  121 +% interpolation des donn???es ??? chaque heure
  122 +dswrf = interp(dswrf,3);
  123 +a = find(dswrf < 0.0);
  124 +dswrf(a) = 0.0;
  125 +
  126 +ncfile = ['uswrf_',yrstr,'_narr_amdgulf.nc'];
  127 +uswrf = nc_varget(ncfile,'USWRF'); % W/m^2
  128 +% interpolation des donn???es ??? chaque heure
  129 +uswrf = interp(uswrf,3);
  130 +
  131 +ncfile = ['dlwrf_',yrstr,'_narr_amdgulf.nc'];
  132 +dlwrf = nc_varget(ncfile,'DLWRF'); % W/m^2
  133 +% interpolation des donn???es ??? chaque heure
  134 +dlwrf = interp(dlwrf,3);
  135 +
  136 +ncfile = ['ulwrf_',yrstr,'_narr_amdgulf.nc'];
  137 +ulwrf = nc_varget(ncfile,'ULWRF'); % W/m^2
  138 +% interpolation des donn???es ??? chaque heure
  139 +ulwrf = interp(ulwrf,3);
  140 +
  141 +ncfile = ['shtfl_',yrstr,'_narr_amdgulf.nc'];
  142 +shtfl = nc_varget(ncfile,'SHTFL'); % W/m^2
  143 +% interpolation des donn???es ??? chaque heure
  144 +shtfl = interp(shtfl,3);
  145 +
  146 +ncfile = ['lhtfl_',yrstr,'_narr_amdgulf.nc'];
  147 +lhtfl = nc_varget(ncfile,'LHTFL'); % W/m^2
  148 +% interpolation des donn???es ??? chaque heure
  149 +lhtfl = interp(lhtfl,3);
  150 +
  151 +% Transmitted shortwave radiative flux
  152 +% that include the extinction due to the ice cover
  153 +% (Data from Canadian Ice Service)
  154 +eval(['load cis_icec_',yrstr])
  155 +swr = (dswrf - uswrf).*(1.0-icec);
  156 +for ii = 1:length(swr)
  157 + if swr(ii) < 0.0;
  158 + swr(ii) = 0.0;
  159 + end
  160 +end
  161 +% No ice effect
  162 +% swr = (dswrf - uswrf);
  163 +
  164 +% Total heat flux (minus the shortwave)
  165 +% positive = towards the ocean surface
  166 +% qt = dlwrf - ulwrf + shtfl + lhtfl; (????)
  167 +% that include the insulation due to the ice cover
  168 +% (Data from D. Barber group)
  169 +qt = (dlwrf - ulwrf + shtfl + lhtfl) .* (1.0-icec);
  170 +% No ice effect
  171 +% qt = dlwrf - ulwrf + shtfl + lhtfl;
  172 +
  173 +figure; plot(time,qt)
  174 +figure; plot(time,swr)
  175 +
  176 +% Creation of the file 'heatflux.dat'
  177 +% file with swr and total heat flux in W/m^2
  178 +init_date = datenum([year 01 01 00 00 00]);
  179 +time = [init_date:1/24:init_date+nod-1/24]';
  180 +timevec = datevec(time);
  181 +tlength = length(time);heatflux = nan .* ones(tlength,8);
  182 +heatflux(:,1:6) = timevec;
  183 +heatflux(:,7) = swr;
  184 +heatflux(:,8) = qt;
  185 +
  186 +fid = fopen(['narr_hourly_heatflux_ice_',yrstr,'.dat'], 'wt');
  187 +fprintf(fid, '%4.0f-%02.0f-%02.0f %02.0f:%02.0f:%02.0f\t%f\t%f\n', ...
  188 + heatflux');
  189 +fclose(fid);
... ...