%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Build a CROCO bulk file % % Extrapole and interpole surface data to get surface boundary % conditions for CROCO (forcing oct_netcdf file) % % Data input format (oct_netcdf): % taux(T, Y, X) % T : time [Months] % Y : Latitude [degree north] % X : Longitude [degree east] % % Data source : IRI/LDEO Climate Data Library % (Atlas of Surface Marine Data 1994) % % http://ingrid.ldgo.columbia.edu/ % http://iridl.ldeo.columbia.edu/SOURCES/.DASILVA/ % % Further Information: % http://www.croco-ocean.org % % This file is part of CROCOTOOLS % % CROCOTOOLS is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published % by the Free Software Foundation; either version 2 of the License, % or (at your option) any later version. % % CROCOTOOLS 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 General Public License for more details. % % You should have received a copy of the GNU General Public License % along with this program; if not, write to the Free Software % Foundation, Inc., 59 Temple Place, Suite 330, Boston, % MA 02111-1307 USA % % Copyright (c) 2005 by Patrick Marchesiello and Pierrick Penven % e-mail:Patrick.Marchesiello@ird.fr % % Updated 2006/09/29 by Pierrick Penven (add a test for the plots) % Updated 2006/10/02 by Pierrick Penven (add the 'tmp file' for % oct_ext_data) % Updated 2006/10/05 by Pierrick Penven (add coads_dir) % Updated 25-Oct-2006 by Pierrick Penven (uwnd and vwnd) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all close all %%%%%%%%%%%%%%%%%%%%% USERS DEFINED VARIABLES %%%%%%%%%%%%%%%%%%%%%%%% % crocotools_param % % Load air-sea parameters % as_consts % % sat : Surface atmospheric temperature % airdens : Surface atmospheric density % w3 : Wind speed at 10 meters % qsea : Sea level specific humidity % rh : relative humidity % precip : precipitation rate % shortrad : Short wave radiation % longrade : Outgoing long wave radiation % sat_file =[coads_dir,'sat.cdf']; sat_name ='sat'; sst_file =[coads_dir,'sst.cdf']; sst_name ='sst'; airdens_file =[coads_dir,'airdens.cdf']; airdens_name ='airdens'; u3_file =[coads_dir,'u3.cdf']; u3_name ='u3'; v3_file =[coads_dir,'v3.cdf']; v3_name ='v3'; w3_file =[coads_dir,'w3.cdf']; w3_name ='w3'; qsea_file =[coads_dir,'qsea.cdf']; qsea_name ='qsea'; rh_file =[coads_dir,'rh.cdf']; rh_name ='rh'; precip_file =[coads_dir,'precip.cdf']; precip_name ='precip'; srf_file =[coads_dir,'shortrad.cdf']; srf_name ='shortrad'; lrf_file =[coads_dir,'longrad.cdf']; lrf_name ='longrad'; taux_file =[coads_dir,'taux.cdf']; taux_name ='taux'; tauy_file =[coads_dir,'tauy.cdf']; tauy_name ='tauy'; % % %%%%%%%%%%%%%%%%%%% END USERS DEFINED VARIABLES %%%%%%%%%%%%%%%%%%%%%%% % % Title % disp(' ') disp(CROCO_title) % % Read in the grid % disp(' ') disp(' Read in the grid...') ncid = netcdf.open(grdname, 'NC_NOWRITE'); Lp=netcdf.inqDimLen(ncid, netcdf.inqDimID(ncid, 'xi_rho')); Mp=netcdf.inqDimLen(ncid, netcdf.inqDimID(ncid, 'eta_rho')); lon=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'lon_rho')); lat=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'lat_rho')); lonu=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'lon_u')); latu=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'lat_u')); lonv=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'lon_v')); latv=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'lat_v')); angle=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'angle')); netcdf.close(ncid); cosa=cos(angle); sina=sin(angle); % % Create the forcing file % disp(' ') disp(' Create the bulk forcing file...') oct_create_bulk(blkname,grdname,CROCO_title,coads_time,coads_cycle); % % Loop on time % ncid = netcdf.open(blkname, 'NC_WRITE'); for tindex=1:length(coads_time) time=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'bulk_time')); netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'tair'), tindex,:,:-1, 1, oct_ext_data(sat_file,sat_name,tindex,...); % [conv] 0-based lon,lat,time,Roa,1); end for tindex=1:length(coads_time) time=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'bulk_time')); % percent -> fraction netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'rhum'), tindex,:,:-1, 1, 0.01*oct_ext_data(rh_file,rh_name,tindex,...); % [conv] 0-based lon,lat,time,Roa,1); end for tindex=1:length(coads_time) time=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'bulk_time')); % mm/(3hour) -> centimeter day-1 netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'prate'), tindex,:,:-1, 1, 0.8*oct_ext_data(precip_file,precip_name,tindex,...); % [conv] 0-based lon,lat,time,Roa,1); end for tindex=1:length(coads_time) time=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'bulk_time')); radlw=oct_ext_data(lrf_file,lrf_name,tindex,... lon,lat,time,Roa,1); netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'radlw'), tindex,:,:-1, 1, radlw); % [conv] 0-based % radlw_in: substract upward gray-body longwave flux % and make it positive downward sst= oct_ext_data(sst_file,sst_name,tindex,... lon,lat,time,Roa,1); lwup=emiss_lw.*sigmaSB.*((sst+CtoK).^4); netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'radlw_in'), tindex,:,:-1, 1, -(radlw-lwup)); % [conv] 0-based end for tindex=1:length(coads_time) time=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'bulk_time')); netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'radsw'), tindex,:,:-1, 1, oct_ext_data(srf_file,srf_name,tindex,...); % [conv] 0-based lon,lat,time,Roa,1); end % % Compute wind rotated and at u,v points % for tindex=1:length(coads_time) time=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'bulk_time')); netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'wspd'), tindex,:,:-1, 1, oct_ext_data(w3_file,w3_name,tindex,...); % [conv] 0-based lon,lat,time,Roa,1); end for tindex=1:length(coads_time) time=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'bulk_time')); uwnd = oct_ext_data(u3_file,u3_name,tindex,... lon,lat,time,Roa,1); vwnd = oct_ext_data(v3_file,v3_name,tindex,... lon,lat,time,Roa,1); u10=oct_rho2u_2d(uwnd.*cosa+vwnd.*sina); v10=oct_rho2v_2d(vwnd.*cosa-uwnd.*sina); netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'uwnd'), tindex,:,:-1, 1, u10); % [conv] 0-based netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'vwnd'), tindex,:,:-1, 1, v10); % [conv] 0-based end % % Compute wind stress rotated and at u,v points % for tindex=1:length(coads_time) time=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'bulk_time')); tx=oct_ext_data(taux_file,taux_name,tindex,... lon,lat,time,Roa,2); ty=oct_ext_data(tauy_file,tauy_name,tindex,... lon,lat,time,Roa,2); netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'sustr'), tindex,:,:-1, 1, oct_rho2u_2d(tx.*cosa + ty.*sina)); % [conv] 0-based netcdf.putVar(ncid, netcdf.inqVarID(ncid, 'svstr'), tindex,:,:-1, 1, oct_rho2v_2d(ty.*cosa - tx.*sina)); % [conv] 0-based end netcdf.close(ncid); % % Make a few plots % if makeplot==1 disp(' ') disp(' Make a few plots...') oct_test_forcing(blkname,grdname,'tair',[1 4 7 10],3,coastfileplot) figure oct_test_forcing(blkname,grdname,'rhum',[1 4 7 10],3,coastfileplot) figure oct_test_forcing(blkname,grdname,'prate',[1 4 7 10],3,coastfileplot) figure oct_test_forcing(blkname,grdname,'uwnd',[1 4 7 10],3,coastfileplot) figure oct_test_forcing(blkname,grdname,'vwnd',[1 4 7 10],3,coastfileplot) figure oct_test_forcing(blkname,grdname,'wspd',[1 4 7 10],3,coastfileplot) figure oct_test_forcing(blkname,grdname,'radlw',[1 4 7 10],3,coastfileplot) figure oct_test_forcing(blkname,grdname,'radlw_in',[1 4 7 10],3,coastfileplot) figure oct_test_forcing(blkname,grdname,'radsw',[1 4 7 10],3,coastfileplot) end % % End % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%