| 1 | gmaze | 1.1 | % | 
| 2 | gmaze | 1.2 | % [EKL] = compute_EKLx(SNAPSHOT) | 
| 3 | gmaze | 1.1 | % | 
| 4 |  |  | % Here we compute the Ekman Layer Depth as: | 
| 5 |  |  | % EKL = 0.7 sqrt( TAUx/RHO )/f | 
| 6 |  |  | % | 
| 7 |  |  | % where: | 
| 8 |  |  | %  TAUx is the amplitude of the zonal surface wind-stress (N/m2) | 
| 9 |  |  | %  RHO is the density of seawater (kg/m3) | 
| 10 |  |  | %  f is the Coriolis parameter (kg/m3) | 
| 11 |  |  | %  EKL is the Ekman layer depth (m) | 
| 12 |  |  | % | 
| 13 |  |  | % Files names are: | 
| 14 |  |  | % INPUT: | 
| 15 |  |  | % ./netcdf-files/<SNAPSHOT>/<netcdf_RHO>.<netcdf_domain>.<netcdf_suff> | 
| 16 |  |  | % ./netcdf-files/<SNAPSHOT>/<netcdf_TAUX>.<netcdf_domain>.<netcdf_suff> | 
| 17 |  |  | % OUTPUT | 
| 18 |  |  | % ./netcdf-files/<SNAPSHOT>/<netcdf_EKLx>.<netcdf_domain>.<netcdf_suff> | 
| 19 |  |  | % | 
| 20 |  |  | % with netcdf_* as global variables | 
| 21 |  |  | % netcdf_EKLx = 'EKLx' by default | 
| 22 |  |  | % | 
| 23 |  |  | % 12/04/06 | 
| 24 |  |  | % gmaze@mit.edu | 
| 25 |  |  |  | 
| 26 | gmaze | 1.2 | function varargout = compute_EKLx(snapshot) | 
| 27 | gmaze | 1.1 |  | 
| 28 |  |  | global sla toshow | 
| 29 |  |  | global netcdf_suff netcdf_domain | 
| 30 |  |  | global netcdf_TAUX netcdf_RHO netcdf_EKLx | 
| 31 |  |  | pv_checkpath | 
| 32 |  |  | global EKL Tx Ty TAU RHO f | 
| 33 |  |  |  | 
| 34 |  |  |  | 
| 35 |  |  | % NETCDF file name: | 
| 36 |  |  | filTx  = netcdf_TAUX; | 
| 37 |  |  | filRHO = netcdf_RHO; | 
| 38 |  |  |  | 
| 39 |  |  | % Path and extension to find them: | 
| 40 |  |  | pathname = strcat('netcdf-files',sla); | 
| 41 |  |  | ext = netcdf_suff; | 
| 42 |  |  |  | 
| 43 |  |  | % Load files: | 
| 44 |  |  | ferfile = strcat(pathname,sla,snapshot,sla,filTx,'.',netcdf_domain,'.',ext); | 
| 45 |  |  | ncTx    = netcdf(ferfile,'nowrite'); | 
| 46 |  |  | Tx      = ncTx{4}(1,:,:); | 
| 47 |  |  |  | 
| 48 |  |  | ferfile = strcat(pathname,sla,snapshot,sla,filRHO,'.',netcdf_domain,'.',ext); | 
| 49 |  |  | ncRHO   = netcdf(ferfile,'nowrite'); | 
| 50 |  |  | RHO     = ncRHO{4}(1,:,:); | 
| 51 |  |  | [RHOlon RHOlat RHOdpt] = coordfromnc(ncRHO); | 
| 52 |  |  |  | 
| 53 |  |  |  | 
| 54 |  |  | %%%%%%%%%%%%%%%%%%%%%%%%%%%% | 
| 55 |  |  | % Get EKL | 
| 56 |  |  | %%%%%%%%%%%%%%%%%%%%%%%%%%%% | 
| 57 |  |  |  | 
| 58 |  |  | % Dim: | 
| 59 |  |  | if toshow, disp('dim'), end | 
| 60 |  |  | nx = length(RHOlon); | 
| 61 |  |  | ny = length(RHOlat); | 
| 62 |  |  | ynz = length(RHOdpt); | 
| 63 |  |  |  | 
| 64 |  |  | % Pre-allocate: | 
| 65 |  |  | if toshow, disp('pre-allocate'), end | 
| 66 |  |  | EKL = zeros(ny,nx); | 
| 67 |  |  |  | 
| 68 |  |  | % Planetary vorticity: | 
| 69 |  |  | f = 2*(2*pi/86400)*sin(RHOlat*pi/180); | 
| 70 |  |  | [a f c]=meshgrid(RHOlon,f,RHOdpt); clear a c | 
| 71 |  |  | f = permute(f,[3 1 2]); | 
| 72 |  |  | f = squeeze(f(1,:,:)); | 
| 73 |  |  |  | 
| 74 |  |  | % Windstress amplitude: | 
| 75 |  |  | TAU = sqrt( Tx.^2 ); | 
| 76 |  |  |  | 
| 77 |  |  | % Ekman Layer Depth: | 
| 78 |  |  | EKL = 0.7* sqrt(TAU ./ RHO) ./f; | 
| 79 |  |  | %EKL = 1.7975 * sqrt( TAU ./ RHO ./ f ); | 
| 80 |  |  |  | 
| 81 |  |  | %%%%%%%%%%%%%%%%%%%%%%%%%%%% | 
| 82 |  |  | % Record | 
| 83 |  |  | %%%%%%%%%%%%%%%%%%%%%%%%%%%% | 
| 84 |  |  | if toshow, disp('record'), end | 
| 85 |  |  |  | 
| 86 |  |  | % General informations: | 
| 87 |  |  | if ~isempty('netcdf_EKLx') | 
| 88 |  |  | netfil = netcdf_EKLx; | 
| 89 |  |  | else | 
| 90 |  |  | netfil = 'EKLx'; | 
| 91 |  |  | end | 
| 92 |  |  | units      = 'm'; | 
| 93 |  |  | ncid       = 'EKLx'; | 
| 94 |  |  | longname   = 'Ekman Layer Depth from TAUx'; | 
| 95 |  |  | uniquename = 'EKLx'; | 
| 96 |  |  |  | 
| 97 |  |  | % Open output file: | 
| 98 |  |  | nc = netcdf(strcat(pathname,sla,snapshot,sla,netfil,'.',netcdf_domain,'.',ext),'clobber'); | 
| 99 |  |  |  | 
| 100 |  |  | % Define axis: | 
| 101 |  |  | nx = length(RHOlon) ; | 
| 102 |  |  | ny = length(RHOlat) ; | 
| 103 |  |  | nz = 1 ; | 
| 104 |  |  |  | 
| 105 |  |  | nc('X') = nx; | 
| 106 |  |  | nc('Y') = ny; | 
| 107 |  |  | nc('Z') = nz; | 
| 108 |  |  |  | 
| 109 |  |  | nc{'X'}            = ncfloat('X'); | 
| 110 |  |  | nc{'X'}.uniquename = ncchar('X'); | 
| 111 |  |  | nc{'X'}.long_name  = ncchar('longitude'); | 
| 112 |  |  | nc{'X'}.gridtype   = nclong(0); | 
| 113 |  |  | nc{'X'}.units      = ncchar('degrees_east'); | 
| 114 |  |  | nc{'X'}(:)         = RHOlon; | 
| 115 |  |  |  | 
| 116 |  |  | nc{'Y'}            = ncfloat('Y'); | 
| 117 |  |  | nc{'Y'}.uniquename = ncchar('Y'); | 
| 118 |  |  | nc{'Y'}.long_name  = ncchar('latitude'); | 
| 119 |  |  | nc{'Y'}.gridtype   = nclong(0); | 
| 120 |  |  | nc{'Y'}.units      = ncchar('degrees_north'); | 
| 121 |  |  | nc{'Y'}(:)         = RHOlat; | 
| 122 |  |  |  | 
| 123 |  |  | nc{'Z'}            = ncfloat('Z'); | 
| 124 |  |  | nc{'Z'}.uniquename = ncchar('Z'); | 
| 125 |  |  | nc{'Z'}.long_name  = ncchar('depth'); | 
| 126 |  |  | nc{'Z'}.gridtype   = nclong(0); | 
| 127 |  |  | nc{'Z'}.units      = ncchar('m'); | 
| 128 |  |  | nc{'Z'}(:)         = RHOdpt(1); | 
| 129 |  |  |  | 
| 130 |  |  | % And main field: | 
| 131 |  |  | nc{ncid}               = ncfloat('Z', 'Y', 'X'); | 
| 132 |  |  | nc{ncid}.units         = ncchar(units); | 
| 133 |  |  | nc{ncid}.missing_value = ncfloat(NaN); | 
| 134 |  |  | nc{ncid}.FillValue_    = ncfloat(NaN); | 
| 135 |  |  | nc{ncid}.longname      = ncchar(longname); | 
| 136 |  |  | nc{ncid}.uniquename    = ncchar(uniquename); | 
| 137 |  |  | nc{ncid}(:,:,:)        = EKL; | 
| 138 |  |  |  | 
| 139 |  |  |  | 
| 140 |  |  |  | 
| 141 |  |  | % Close files: | 
| 142 |  |  | close(ncTx); | 
| 143 |  |  | close(ncRHO); | 
| 144 |  |  | close(nc); | 
| 145 |  |  |  | 
| 146 |  |  |  | 
| 147 | gmaze | 1.2 |  | 
| 148 |  |  | % Output: | 
| 149 |  |  | output = struct('EKL',EKL,'lat',RHOlat,'lon',RHOlon); | 
| 150 |  |  | switch nargout | 
| 151 |  |  | case 1 | 
| 152 |  |  | varargout(1) = {output}; | 
| 153 |  |  | end |