Index: /issm/trunk-jpl/test/MITgcm/code_4003/CPP_OPTIONS.h =================================================================== --- /issm/trunk-jpl/test/MITgcm/code_4003/CPP_OPTIONS.h (revision 26663) +++ /issm/trunk-jpl/test/MITgcm/code_4003/CPP_OPTIONS.h (revision 26664) @@ -54,8 +54,4 @@ #define INCLUDE_PHIHYD_CALCULATION_CODE -C o Include/exclude sound speed calculation code -C o (Note that this is a diagnostic from Del Grasso algorithm, not derived from EOS) -#undef INCLUDE_SOUNDSPEED_CALC_CODE - C-- Vertical mixing code options: @@ -78,8 +74,4 @@ #undef EXCLUDE_PCELL_MIX_CODE -C o Exclude/allow to use isotropic 3-D Smagorinsky viscosity as diffusivity -C for tracers (after scaling by constant Prandtl number) -#undef ALLOW_SMAG_3D_DIFFUSIVITY - C-- Time-stepping code options: @@ -93,7 +85,4 @@ #undef ALLOW_ADAMSBASHFORTH_3 -C o Include/exclude Quasi-Hydrostatic Stagger Time-step AdamsBashforth code -#undef ALLOW_QHYD_STAGGER_TS - C-- Model formulation options: @@ -105,7 +94,4 @@ C this implies that grid-cell thickness (hFactors) varies with time #define NONLIN_FRSURF -C o Disable code for rStar coordinate and/or code for Sigma coordinate -c#define DISABLE_RSTAR_CODE -c#define DISABLE_SIGMA_CODE C o Include/exclude nonHydrostatic code @@ -117,5 +103,5 @@ C-- Algorithm options: -C o Include/exclude code for Non Self-Adjoint (NSA) conjugate-gradient solver +C o Use Non Self-Adjoint (NSA) conjugate-gradient solver #undef ALLOW_CG2D_NSA @@ -131,39 +117,8 @@ C-- Retired code options: -C o ALLOW isotropic scaling of harmonic and bi-harmonic terms when -C using an locally isotropic spherical grid with (dlambda) x (dphi*cos(phi)) -C *only for use on a lat-lon grid* -C Setting this flag here affects both momentum and tracer equation unless -C it is set/unset again in other header fields (e.g., GAD_OPTIONS.h). -C The definition of the flag is commented to avoid interference with -C such other header files. -C The preferred method is specifying a value for viscAhGrid or viscA4Grid -C in data which is then automatically scaled by the grid size; -C the old method of specifying viscAh/viscA4 and this flag is provided -C for completeness only (and for use with the adjoint). -c#define ISOTROPIC_COS_SCALING - -C o This flag selects the form of COSINE(lat) scaling of bi-harmonic term. -C *only for use on a lat-lon grid* -C Has no effect if ISOTROPIC_COS_SCALING is undefined. -C Has no effect on vector invariant momentum equations. -C Setting this flag here affects both momentum and tracer equation unless -C it is set/unset again in other header fields (e.g., GAD_OPTIONS.h). -C The definition of the flag is commented to avoid interference with -C such other header files. -c#define COSINEMETH_III - -C o Use "OLD" UV discretisation near boundaries (*not* recommended) -C Note - only works with pkg/mom_fluxform and "no_slip_sides=.FALSE." -C because the old code did not have no-slip BCs -#undef OLD_ADV_BCS - C o Use LONG.bin, LATG.bin, etc., initialization for ini_curviliear_grid.F C Default is to use "new" grid files (OLD_GRID_IO undef) but OLD_GRID_IO C is still useful with, e.g., single-domain curvilinear configurations. #undef OLD_GRID_IO - -C o Use old EXTERNAL_FORCING_U,V,T,S subroutines (for backward compatibility) -#undef USE_OLD_EXTERNAL_FORCING C-- Other option files: Index: /issm/trunk-jpl/test/MITgcm/code_4003/DIAGNOSTICS_SIZE.h =================================================================== --- /issm/trunk-jpl/test/MITgcm/code_4003/DIAGNOSTICS_SIZE.h (revision 26663) +++ /issm/trunk-jpl/test/MITgcm/code_4003/DIAGNOSTICS_SIZE.h (revision 26664) @@ -23,5 +23,4 @@ PARAMETER( diagSt_size = 10*Nr ) - CEH3 ;;; Local Variables: *** CEH3 ;;; mode:fortran *** Index: /issm/trunk-jpl/test/MITgcm/code_4003/SIZE.h =================================================================== --- /issm/trunk-jpl/test/MITgcm/code_4003/SIZE.h (revision 26663) +++ /issm/trunk-jpl/test/MITgcm/code_4003/SIZE.h (revision 26664) @@ -43,6 +43,6 @@ INTEGER Nr PARAMETER ( - & sNx = 1, - & sNy = 50, + & sNx = 3, + & sNy = 25, & OLx = 3, & OLy = 3, @@ -50,5 +50,5 @@ & nSy = 1, & nPx = 1, - & nPy = 4, + & nPy = 8, & Nx = sNx*nSx*nPx, & Ny = sNy*nSy*nPy, Index: /issm/trunk-jpl/test/MITgcm/input_4003/data =================================================================== --- /issm/trunk-jpl/test/MITgcm/input_4003/data (revision 26663) +++ /issm/trunk-jpl/test/MITgcm/input_4003/data (revision 26664) @@ -62,9 +62,8 @@ #Time stepping parameters &PARM03 -#nIter0=2880, -#nTimeSteps=288, -#startTime=0., + startTime=0., + nTimeSteps=192, #endTime=2592000., - deltaT=300.0, + deltaT=450.0, forcing_In_AB = .FALSE., abEps=0.1, @@ -72,12 +71,4 @@ chkptFreq = 2592000., dumpFreq = 86400., - monitorFreq= 7200., - monitorSelect=1, - monitorFreq=1800., - dumpFreq = 10800., -#- for testing purpose: - nIter0=2898, - nTimeSteps=20, -#dumpFreq = 1200., monitorFreq=1., & @@ -88,5 +79,5 @@ usingSphericalPolarGrid=.TRUE., delR=90*10., - delX=1*.125, + delX=3*.125, delY=200*.0078125, xgOrigin = 0., Index: /issm/trunk-jpl/test/MITgcm/input_4003/data.diagnostics =================================================================== --- /issm/trunk-jpl/test/MITgcm/input_4003/data.diagnostics (revision 26664) +++ /issm/trunk-jpl/test/MITgcm/input_4003/data.diagnostics (revision 26664) @@ -0,0 +1,69 @@ +# Diagnostic Package Choices +#----------------- +# for each output-stream: +# filename(n) : prefix of the output file name (only 8.c long) for outp.stream n +# frequency(n):< 0 : write snap-shot output every |frequency| seconds +# > 0 : write time-average output every frequency seconds +# timePhase(n) : write at time = timePhase + multiple of |frequency| +# averagingFreq(n) : frequency (in s) for periodic averaging interval +# averagingPhase(n): phase (in s) for periodic averaging interval +# repeatCycle(n) : number of averaging intervals in 1 cycle +# levels(:,n) : list of levels to write to file (Notes: declared as REAL) +# when this entry is missing, select all common levels of this list +# fields(:,n) : list of diagnostics fields (8.c) (see "available_diagnostics.log" +# file for the list of all available diag. in this particular config) +#----------------- + &DIAGNOSTICS_LIST +# diag_mnc = .FALSE., + dumpAtLast = .TRUE., + fields(1:13,1) = 'ETAN ','RSURF ','oceTAUX ','oceTAUY ', + 'oceQnet ','oceFWflx','MXLDEPTH', + 'SHIfwFlx','SHIhtFlx','SHIgammT','SHIgammS', + 'SHI_mass','SHIuStar', +# 'surForcT','surForcS','TFLUX ','SFLUX ','oceFreez', +# 'TRELAX ','SRELAX ', +# fields(1,1)='ETAN' + filename(1) = 'surfDiag', + frequency(1) = 86400., +#- + fields(1:10,2) = 'UVEL ','VVEL ','WVEL ', + 'THETA ','SALT ','RHOAnoma', 'CONVADJ' + 'hFactorC','hFactorW','hFactorS', + filename(2) = 'dynDiag', + fileFlags(2) = ' h ', + frequency(2) = 86400., +#- + fields(1:7,3) = 'ETAN ','RSURF ', + 'oceQnet ','oceFWflx', + 'SHIfwFlx','SHIhtFlx','SHIuStar', +# filename(3) = 'surfInst', + frequency(3) = -86400., + timePhase(3) = 0., +#- + fields(1:8,4) = 'UVEL ','VVEL ','WVEL ', + 'THETA ','SALT ', + 'hFactorC','hFactorW','hFactorS', +# filename(4) = 'dynInst', + frequency(4) = -86400., + timePhase(4) = 0., + & + +#-------------------- +# Parameter for Diagnostics of per level statistics: +#-------------------- +# diagSt_mnc (logical): write stat-diags to NetCDF files (default=diag_mnc) +# diagSt_regMaskFile : file containing the region-mask to read-in +# nSetRegMskFile : number of region-mask sets within the region-mask file +# set_regMask(i) : region-mask set-index that identifies the region "i" +# val_regMask(i) : region "i" identifier value in the region mask +#--for each output-stream: +# stat_fName(n) : prefix of the output file name (max 80c long) for outp.stream n +# stat_freq(n):< 0 : write snap-shot output every |stat_freq| seconds +# > 0 : write time-average output every stat_freq seconds +# stat_phase(n) : write at time = stat_phase + multiple of |stat_freq| +# stat_region(:,n) : list of "regions" (default: 1 region only=global) +# stat_fields(:,n) : list of selected diagnostics fields (8.c) in outp.stream n +# (see "available_diagnostics.log" file for the full list of diags) +#-------------------- + &DIAG_STATIS_PARMS + & Index: /issm/trunk-jpl/test/MITgcm/input_4003/data.obcs =================================================================== --- /issm/trunk-jpl/test/MITgcm/input_4003/data.obcs (revision 26664) +++ /issm/trunk-jpl/test/MITgcm/input_4003/data.obcs (revision 26664) @@ -0,0 +1,32 @@ +# Open-boundaries + &OBCS_PARM01 + OB_Jnorth=3*200, + useOBCSprescribe=.TRUE., +#OBWuFile='uVel_obc.bin', + OBNvFile='vVel_obc.bin', + OBNtFile='temp_obc.bin', + OBNsFile='salt_obc.bin', + OBCS_u1_adv_T=1, + OBCS_u1_adv_S=1, +#-- + useOBCSsponge=.TRUE., + useLinearSponge = .TRUE., +#- default is =T for these specific switches: +#OBCSsponge_VatNS=.TRUE., + OBCSsponge_UatNS=.FALSE., +#-- + OBCSprintDiags = .FALSE., + & + +# Orlanski parameters (skipped with #undef ALLOW_ORLANSKI) +#&OBCS_PARM02 +#Cmax=0.45, +#cVelTimeScale=1000., +#& + +# Sponge layer parameters + &OBCS_PARM03 + spongeThickness = 10, + Vrelaxobcsbound = 864000., + Urelaxobcsbound = 864000., + & Index: /issm/trunk-jpl/test/MITgcm/input_4003/data.pkg =================================================================== --- /issm/trunk-jpl/test/MITgcm/input_4003/data.pkg (revision 26664) +++ /issm/trunk-jpl/test/MITgcm/input_4003/data.pkg (revision 26664) @@ -0,0 +1,6 @@ +# Packages (lines beginning "#" are comments) + &PACKAGES + useOBCS=.TRUE., + useShelfIce=.TRUE., + useDiagnostics = .TRUE., + & Index: /issm/trunk-jpl/test/NightlyRun/test4003.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test4003.m (revision 26663) +++ /issm/trunk-jpl/test/NightlyRun/test4003.m (revision 26664) @@ -5,11 +5,10 @@ steps=1:12; steps=1:6; -steps=7; final_time=1/365; %To download and recompile MITgcm from scratch: -%!rm -rf ${ISSM_DIR}/test/MITgcm/install -%!rm -rf ${ISSM_DIR}/test/MITgcm/build/* -%!rm -rf Models +!rm -rf ${ISSM_DIR}/test/MITgcm/install +!rm -rf ${ISSM_DIR}/test/MITgcm/build/* +!rm -rf Models %Organizer @@ -21,13 +20,13 @@ % {{{ Parameters: if perform(org,'Parameters'), - Nx=1; %number of longitude cells + Nx=3; %number of longitude cells Ny=200; %number of latitude cells Nz=90; %number of MITgcm vertical cells nPx=1; %number of MITgcm processes to use in x direction - nPy=4; %number of MITgcm processes to use in y direction - xgOrigin=1; %origin of longitude - ygOrigin=1; %origin of latitude - dLong=1; %longitude grid spacing - dLat=1; %latitude grid spacing + nPy=8; %number of MITgcm processes to use in y direction + xgOrigin=0; %origin of longitude + ygOrigin=-75.5; %origin of latitude + dLong=.125; %longitude grid spacing + dLat=dLong/16; %latitude grid spacing delZ=10; %thickness of vertical levels icefront_position_ratio=.75; @@ -47,5 +46,5 @@ obcUvel = -0.1; % open boundary velocity (m/s) - MITgcmDeltaT=600; % MITgcm time step in seconds + MITgcmDeltaT=450; % MITgcm time step in seconds y2s=31536000; % year to seconds conversion, i.e., seconds per year @@ -55,5 +54,5 @@ async_step_MITgcm_multiplier=1; % used to reduce run time for MItgcm - % bedrock/bathymetry + % bathymetry hmax=1000; trough_depth=200; @@ -76,22 +75,114 @@ loaddata(org,'Parameters'); - %create lat,long + + v0 = 2e3; + h0 = 800; + + hfacMin = 0.2; + + eos = 'jmd95z'; + prec = 'real*8'; + + %create lat,lon + latg = ygOrigin+[0:Ny-1]*dLat; + latc = latg+dLat/2; + long = xgOrigin+[0:Nx-1]*dLong; + lonc = long+dLong/2; + zC=-delZ*([1:Nz]-0.5); + zF=-delZ*[0:Nz]; + lat=(ygOrigin+dLat/2):dLat:(ygOrigin+Ny*dLat); - long=(xgOrigin+dLong/2):dLong:(xgOrigin+Nx*dLong); - [lat long]=meshgrid(lat,long); - - longmin=min(long(:)); - longmax=max(long(:)); + lon=(xgOrigin+dLong/2):dLong:(xgOrigin+Nx*dLong); + [lat lon]=meshgrid(lat,lon); + + longmin=min(lon(:)); + longmax=max(lon(:)); latmin=min(lat(:)); latmax=max(lat(:)); - %create bedrock/bathymetry: - bedrock=zeros(Nx,Ny); - bedrock=hmax-deltah*tanh(pi*(2*(lat-latmin)./(latmax-latmin)-1))+ ... - trough_depth*cos(2*pi*long./(longmax-longmin)); +% Gravity +gravity= 9.81; +rhoConst= 1030; + +% Nominal depth of model (meters) +H = -900; %water depth in the ice shelf cavity +Hmin = -600; % deepest point of cavern +Hmax = -300; % shallowest point of cavern +jEnd = Ny*3/4; % where ice-shelf ends +dHdy = (Hmax-Hmin)/dLat/(jEnd-2); %Slope of ice shelf + + %create bathymetry: + bathymetry = ones(Nx,Ny)*H; + bathymetry(:,1) = 0; %save bathymetry file for MITgcm - bathymetry=bedrock-sea_level; - savedata(org,lat,long,bathymetry); + savedata(org,lat,lon,bathymetry); + +j2=jEnd+1; +hIce=bathymetry; +for i=1:Nx + hIce(i,:)=Hmin+dHdy*[-1:Ny-2]*dLat; + hIce(i,1)=0; hIce(i,j2:Ny)=0; +end + +var=([1:Ny]-2)/(jEnd-2); +dMdt_fy=bathymetry; +for i=1:Nx + dMdt_fy(i,:)=-cos(pi*var); + dMdt_fy(i,1)=0; dMdt_fy(i,j2:Ny)=0; +end + +regMsk=ones(Nx,Ny); +regMsk(:,1)=0; regMsk(:,j2:Ny)=2; + +%- rate of change due to ice-stream dynamics +rateDyn=rhoConst*0.1/3600; sfx='r02'; + +dMdt=rateDyn*dMdt_fy; + +dz = delZ*ones(1,Nz); +zgp1 = [0,cumsum(dz)]; +zc = .5*(zgp1(1:end-1)+zgp1(2:end)); +zg = zgp1(1:end-1); +dz = diff(zgp1); + +T_sfc = -1.9; +T_bot = 2; +del_T = (T_bot - T_sfc)/(59*delZ); +tref=0*dz; +for k = 1:Nz; + tref(k) = T_sfc + del_T*((k-20)*delZ); + tref(k)= max(T_sfc,min(tref(k),T_bot)); +end + +S_sfc = 34.2; +S_bot = 34.7; +del_S = (S_bot - S_sfc)/(59*delZ); +sref=0*dz; +for k = 1:Nz; + sref(k) = S_sfc + del_S*((k-20)*delZ); + sref(k)= max(S_sfc,min(sref(k),S_bot)); +end +pEOS=-rhoConst*gravity*zC; % in Pa +pEOS=pEOS*1.e-4; % in dBar +rhoAn=densjmd95(sref,tref,pEOS); +rhoAn=rhoAn-rhoConst; + +pF=-rhoConst*gravity*zF*1.e-4; % in dBar +rhoUp=densjmd95(sref,tref,pF(2:end)); +rhoDw=densjmd95(sref,tref,pF(1:Nz)); +dRho=rhoUp(1:Nz-1)-rhoDw(2:Nz); +NSq=-gravity*dRho/delZ/rhoConst; + +mnV=min(NSq); MxV=max(NSq); Avr=mean(NSq); + +zax=[1:Nz]; + +v1=2.5e-2; +var=1+Nz-2*zax; var=var/(Nz-1); +vobc=v1*var; + +rhoAvr=rhoConst-1.345; +mIce=-rhoAvr*hIce; end @@ -153,11 +244,50 @@ !\mkdir run !\cp ../MITgcm/build/mitgcmuv run - !\cp ../MITgcm/install/verification/shelfice_2d_remesh/input/* run !\cp ../MITgcm/input_4003/* run - !\cp ../MITgcm/input_4003/data_uncoupled run/data !\cp ../MITgcm/input_4003/eedata_uncoupled run/eedata cd run + + namF='bathy_flat.bin'; + fid=fopen(namF,'w','b'); fwrite(fid,bathymetry,prec);fclose(fid); + + namF='shelficeTopo.Lin.bin'; + fid=fopen(namF,'w','b'); fwrite(fid,hIce,prec);fclose(fid); + + namF='under_Ice_mask.bin'; + fid=fopen(namF,'w','b'); fwrite(fid,regMsk,prec);fclose(fid); + + namF=sprintf('%s.%s.%s','shelfice_dMdt',sfx,'bin'); + fid=fopen(namF,'w','b'); fwrite(fid,dMdt,prec);fclose(fid); + + namF='temp_obc.bin'; + tref=[tref; tref; tref]; + fid=fopen(namF,'w','b'); fwrite(fid,tref,prec);fclose(fid); + + namF='salt_obc.bin'; + sref=[sref; sref; sref]; + fid=fopen(namF,'w','b'); fwrite(fid,sref,prec);fclose(fid); + + namF='vVel_obc.bin'; + vobc=[vobc; vobc; vobc]; + fid=fopen(namF,'w','b'); fwrite(fid,vobc,prec);fclose(fid); + + var=zeros(Nx,Ny,Nz); + for i=1:Nx, for j=1:Ny + var(i,j,:)=tref(1,:); + end, end + namF='temp_ini.bin'; + fid=fopen(namF,'w','b'); fwrite(fid,var,prec);fclose(fid); + + for i=1:Nx, for j=1:Ny + var(i,j,:)=sref(1,:); + end, end + namF='salt_ini.bin'; + fid=fopen(namF,'w','b'); fwrite(fid,var,prec);fclose(fid); + + namF='shelficeMass.Lin.bin'; + fid=fopen(namF,'w','b'); fwrite(fid,mIce,prec);fclose(fid); + eval(['!mpirun -np ' int2str(nPx*nPy) ' ./mitgcmuv']); - cd .. + cd .. end % }}} @@ -671,43 +801,2 @@ end % }}} - -%Fields and tolerances to track changes -%fnm=['run/SHICE_fwFlux.0000011400.data']; -%melting_rate_1=readbin(fnm,[Nx Ny]); -%fnm=['run/SHICE_fwFlux.0000011450.data']; -%melting_rate_2=readbin(fnm,[Nx Ny]); -%fnm=['run/SHICE_fwFlux.0000011500.data']; -%melting_rate_3=readbin(fnm,[Nx Ny]); -%fnm=['run/SHICE_fwFlux.0000011550.data']; -%melting_rate_4=readbin(fnm,[Nx Ny]); -%field_names ={'Base1','Melting1','Vx2','Vy2','Thickness2','Base2','MaskOceanLevelset2','FloatingiceMeltingRate2',... -% 'Melting2','Vx3','Vy3','Thickness3','Base3','MaskOceanLevelset3','FloatingiceMeltingRate3',... -% 'Melting3','Vx4','Vy4','Thickness4','Base4','MaskOceanLevelset4','FloatingiceMeltingRate4','Melting4'}; -%field_tolerances={2e-13,1e-13,7e-13,7e-13,1e-13,1e-13,1e-13,1e-12,... -% 1e-13, 1e-12, 1e-12, 1e-13, 1e-13, 1e-13, 1e-12,... -% 1e-13, 1e-12, 1e-12, 1e-13, 1e-13, 1e-13, 1e-12, 1e-13 }; -%field_values={... -% (md.results.TransientSolution(1).Base),... -% (melting_rate_1(:)),... -% (md.results.TransientSolution(2).Vx),... -% (md.results.TransientSolution(2).Vy),... -% (md.results.TransientSolution(2).Thickness),... -% (md.results.TransientSolution(2).Base),... -% (md.results.TransientSolution(2).MaskOceanLevelset),... -% (md.results.TransientSolution(2).BasalforcingsFloatingiceMeltingRate),... -% (melting_rate_2(:)),... -% (md.results.TransientSolution(3).Vx),... -% (md.results.TransientSolution(3).Vy),... -% (md.results.TransientSolution(3).Thickness),... -% (md.results.TransientSolution(3).Base),... -% (md.results.TransientSolution(3).MaskOceanLevelset),... -% (md.results.TransientSolution(3).BasalforcingsFloatingiceMeltingRate),... -% (melting_rate_3(:)),... -% (md.results.TransientSolution(4).Vx),... -% (md.results.TransientSolution(4).Vy),... -% (md.results.TransientSolution(4).Thickness),... -% (md.results.TransientSolution(4).Base),... -% (md.results.TransientSolution(4).MaskOceanLevelset),... -% (md.results.TransientSolution(4).BasalforcingsFloatingiceMeltingRate),... -% (melting_rate_4(:)),... -% };