upgrade "global_ocean.cs32x15.viscA4" test

to test "addMassFile" (and selectAddFluid=1,) with implicDiv2DFlow= 0.5

verification/global_ocean.cs32x15/code/CPP_OPTIONS.h

@@ -28,7 +28,7 @@ C o Allow to account for heating due to friction (and momentum dissipation)
 
 C o Allow mass source or sink of Fluid in the interior
 C   (3-D generalisation of oceanic real-fresh water flux)
-#undef ALLOW_ADDFLUID
+#define ALLOW_ADDFLUID
 
 C o Include pressure loading code
 #define ATMOSPHERIC_LOADING

verification/global_ocean.cs32x15/input.viscA4/data

@@ -26,12 +26,14 @@
  vectorInvariantMomentum=.TRUE.,
  implicitFreeSurface=.TRUE.,
  exactConserv=.TRUE.,
+#-- just to test (not physically justified):
  selectP_inEOS_Zc=3,
  nonHydrostatic=.TRUE.,
  selectNHfreeSurf=1,
  useNHMTerms=.TRUE.,
- hFacInf=0.2,
- hFacSup=2.0,
+ implicDiv2DFlow=0.5,
+ selectAddFluid=1,
+#---
  useRealFreshWaterFlux=.TRUE.,
  allowFreezing=.TRUE.,
  hFacMin=.1,
@@ -64,9 +66,11 @@
  taveFreq    =31104000.,
  dumpFreq    =31104000.,
  monitorFreq =86400.,
- periodicExternalForcing=.TRUE.,
- externForcingPeriod=2592000.,
- externForcingCycle=31104000.,
+#-- just to test: switch off periodic forcing --> constant forcing
+# periodicExternalForcing=.TRUE.,
+# externForcingPeriod=2592000.,
+# externForcingCycle=31104000.,
+#--
 # 2 months restoring timescale for temperature
  tauThetaClimRelax = 5184000.,
 # 2yrs restoring timescale for salinity
@@ -95,5 +99,8 @@
  thetaClimFile  ='lev_surfT_cs_12m.bin',
  saltClimFile   ='lev_surfS_cs_12m.bin',
  surfQnetFile   ='shiQnet_cs32.bin',
- EmPmRFile      ='shiEmPR_cs32.bin',
+#EmPmRFile      ='shiEmPR_cs32.bin',
+#- just to test using addMassFile:
+ addMassFile    ='addMass_k1.bin',
+#EmPmRFile      ='large_Precip.bin',
  &

verification/global_ocean.cs32x15/input.viscA4/data.diagnostics

@@ -68,5 +68,10 @@
 #                (see "available_diagnostics.log" file for the full list of diags)
 #--------------------
  &DIAG_STATIS_PARMS
+# diagSt_mnc    = .TRUE.,
+ stat_fields(1:4,1)  = 'ETAN    ','PHI_SURF','RSURF   ','PHI_NH  ',
+  stat_fName(1) = 'dynStDiag',
+   stat_freq(1) = 3600.,
+  stat_phase(1) = 0.,
  &
 

verification/global_ocean.cs32x15/input.viscA4/mk_addMass.m

@@ -0,0 +1,74 @@
+
+kwr=0;
+rDir='../input/';
+
+inpSfx='.bin'; prec='real*8';
+outSfx=inpSfx; precOut=prec;
+%outSfx='.r4_bin'; precOut='real*4';
+
+gDir='grid/';
+G=load_grid(gDir,0);
+nx=G.dims(1); ny=G.dims(2); nr=G.dims(3);
+
+mskC=ceil(G.hFacC);
+rAc=G.rAc.*mskC(:,:,1); areaG=sum(rAc(:));
+
+%for j=1:length(listInpF)
+% namF=char(listInpF(j));
+  namF='shiEmPR_cs32'; n3d=1; j=1;
+  fprintf('-- processing file (#%i) %s ',j,namF);
+ %fprintf('(length= %i ) :',length(namF));
+  fprintf(': n3d= %i ',n3d);
+  var=rdda([rDir,namF,inpSfx],[nx ny n3d],1,prec,'b');
+  fprintf(' Make changes\n');
+  emp0=var;
+
+  var=emp0; titv=namF;
+  var=var.*mskC(:,:,1); mnV=min(var(:)); MxV=max(var(:));
+  var=var.*rAc; Avr=sum(var(:))/areaG;
+  fprintf(' field %15s : min,max= %e, %e , Mean= %e\n',titv,mnV,MxV,Avr);
+
+  %- Multiply by 10 and substract 1.e-7 (~ 1.cm/d , excess of Precip):
+  emp1=10*emp0-ones(nx,ny)*1.e-7;
+
+  var=emp1; titv='- New-';
+  var=var.*mskC(:,:,1); mnV=min(var(:)); MxV=max(var(:));
+  var=var.*rAc; Avr=sum(var(:))/areaG;
+  fprintf(' field %15s : min,max= %e, %e , Mean= %e\n',titv,mnV,MxV,Avr);
+
+  %-- convert to addMass:
+  rhoFresh=1000.;
+  addM=-rhoFresh*emp1; %- mass input in kg/m2/s
+  addM=addM.*rAc;      %- in kg/s
+
+  var=addM; titv='addMass';
+  var=var.*mskC(:,:,1); mnV=min(var(:)); MxV=max(var(:));
+  var=var.*rAc; Avr=sum(var(:))/areaG;
+  fprintf(' field %15s : min,max= %e, %e , Mean= %e\n',titv,mnV,MxV,Avr);
+
+  if kwr > 0,
+   fnam=[namF,outSfx];
+   fnam=['large_Precip',outSfx]; var=emp1;
+   fprintf('  writing file: %s ...',fnam);
+   fid=fopen(fnam,'w','b');
+   fwrite(fid,var,precOut); fclose(fid);
+   fprintf(' done\n');
+
+   fnam=['addMass_k1',outSfx]; var=zeros(nx,ny,nr); var(:,:,1)=addM;
+   fprintf('  writing file: %s ...',fnam);
+   fid=fopen(fnam,'w','b');
+   fwrite(fid,var,precOut); fclose(fid);
+   fprintf(' done\n');
+  end
+
+%end
+%return
+
+figure(1);clf; colormap('jet');
+ var=emp0; ccB=[-1 1]*5.e-8;
+%var=emp1; ccB=[-1 1]*5.e-7;
+ imagesc([1:nx],[1:ny],var');
+ set(gca,'YDir','normal');colorbar;
+ caxis(ccB);
+
+return