Welcome, Guest
Username: Password: Remember me
Forum on HARMONIE Surface development

TOPIC: Incorrect values in soil moisture output

Incorrect values in soil moisture output 2 years 3 weeks ago #1980

  • Bert van Ulft
  • Bert van Ulft's Avatar
  • OFFLINE
  • Expert Boarder
  • Posts: 98
  • Thank you received: 22
Dear all,

in cy40h1.2.tg2 experiments the soil moisture contains points with a value of 1.0E+19 (the red spots in attached image). These values do not occur in a cold start cycle, but in the next cycle with data assimilation these values occur. They show up along the coast and lakes (glaciers?), only in the top soil layer, and only for liquid water (X00[12]WG1), not for ice. The attached image is for patch 2, it also exists in patch 1, but only for 1 point.
It looks a bit like an incorrect missing value is used somewhere, as the default missing value, used above sea, is 1.0E+20. I don't know if it has any effect on the forecast.
The problem can easily be shown with gl (e.g. gl -f ICMSHFULL+0000.sfx -m | grep WG[123]). Correct output after cold start:
X001WG1 >  20171002_12:00+000h00m         000     1.000E-003  241.996E-003  421.737E-003   42.843E-003
X002WG1 >  20171002_12:00+000h00m         000     1.000E-003  242.193E-003  421.737E-003   42.612E-003
X001WG2 >  20171002_12:00+000h00m         000     1.000E-003  226.766E-003  423.184E-003   43.052E-003
X002WG2 >  20171002_12:00+000h00m         000    13.977E-003  219.704E-003  416.031E-003   44.565E-003
X001WG3 >  20171002_12:00+000h00m         000     1.000E-003  223.815E-003  447.865E-003   42.642E-003
X002WG3 >  20171002_12:00+000h00m         000    51.360E-003  230.481E-003  434.635E-003   41.223E-003

Incorrect output in cycle with data assimilation:
X001WG1 >  20171002_15:00+000h00m         000     1.912E-003   83.917E+012   10.000E+018   28.968E+015
X002WG1 >  20171002_15:00+000h00m         000    11.602E-003   44.448E+015   10.000E+018  665.209E+015
X001WG2 >  20171002_15:00+000h00m         000     1.000E-003  225.869E-003  420.351E-003   42.857E-003
X002WG2 >  20171002_15:00+000h00m         000     1.000E-003  218.335E-003  415.227E-003   46.572E-003
X001WG3 >  20171002_15:00+000h00m         000     2.176E-003  223.893E-003  442.194E-003   42.614E-003
X002WG3 >  20171002_15:00+000h00m         000    51.369E-003  230.524E-003  432.472E-003   41.207E-003

It looks like it is an error in the surface data assimilation. The soil moisture at the end of the cold start cycle looks OK, while the ICMSHANAL+0000.sfx file of the first cycle with data assimilation already contains the incorrect values. As far as I have seen all points that get value 1.e19 had a missing value before the data assimilation step.
Has anybody else seen this problem and perhaps already fixed it?

best wishes,

Bert
Attachments:
Last Edit: 2 years 3 weeks ago by Bert van Ulft. Reason: More info

Incorrect values in soil moisture output 2 years 3 weeks ago #1993

  • Bert van Ulft
  • Bert van Ulft's Avatar
  • OFFLINE
  • Expert Boarder
  • Posts: 98
  • Thank you received: 22
Dear all,

Patrick let me know that Mariken and Trygve had already detected and fixed this error in the harmonie-40h1.2_EKF branch in hirlam.org/trac/changeset/16413/branches...m_nature_isba_oi.F90.
The modifications exclude points with missing values using where statements. The relevant lines are:
Old r15742:
370	 	  ZWSINC(:) = ZWS(:) - PWS(:)*(RD1*XRHOLW)     
371	 	  ZWPINC(:) = ZWP(:) - PWP(:)*(PD2(:)*XRHOLW)  
372	 	  ZTLINC(:) = ZTL(:) - PTL(:)*(PD2(:)*XRHOLW)  
373	 	 
374	 	  PWS(:)  = ZWS(:)/(RD1*XRHOLW) 
375	 	  PWP(:)  = ZWP(:)/(PD2(:)*XRHOLW) 
376	 	  PTL(:)  = ZTL(:)/(PD2(:)*XRHOLW) 

Corrected code r16413:
 	387	  WHERE(PWS(:)/=XUNDEF) 
 	388	     ZWSINC(:) = ZWS(:) - PWS(:)*(RD1*XRHOLW)     
 	389	     PWS(:)    = ZWS(:)/(RD1*XRHOLW) 
 	390	  ENDWHERE 
 	391	  WHERE(PWP(:)/=XUNDEF) 
 	392	     ZWPINC(:) = ZWP(:) - PWP(:)*(PD2(:)*XRHOLW)  
 	393	     PWP(:)    = ZWP(:)/(PD2(:)*XRHOLW) 
 	394	  ENDWHERE 
 	395	  WHERE(PTL(:)/=XUNDEF) 
 	396	     ZTLINC(:) = ZTL(:) - PTL(:)*(PD2(:)*XRHOLW)  
 	397	     PTL(:)    = ZTL(:)/(PD2(:)*XRHOLW) 
 	398	  ENDWHERE 

Thanks Patrick, Mariken and Trygve.

Bert
Time to create page: 0.061 seconds