Several issues in modlsm.f90

[adoyenne]

Hello,

Here, I would like to describe several issues I revealed in the LSM module modlsm.f90 and my tentative suggestions to solve them. Note that each of them results in a segmentation fault:

Probably, some of these problems are happen in some grid points only, so they are domain-specific. Also, some of these issues may be interconnected, so resolving one of them could potentially address several others. In my case, I faced with them when simulating a fine domain:

# fine domain 156.25m resolution
     x0 = 874468 # 3.25 deg E
    y0 = 926741 # 51.5 deg N   t
    itot = 1536
    jtot = 768
    dx   = 156.25
    dy   = 156.25
    nprocx = 16 # 256 cores = 2 nodes
    nprocy = 16

1. Prognostic phiw(i,j,k) (soil water content) that is read from restart files when the model runs from warmstart for some grid points may turn into highly negative values causing overflow and segmentation fault in subroutine calc_thermal_properties

! Heat conductivity at saturation (from IFS code..)    
        
   lambda_T_sat =   gamma_T_matrix**(1.-theta_sat(si)) &
                               * gamma_T_water**phiw(i,j,k) &
                               * 2.2**(theta_sat(si) - phiw(i,j,k))

My temporal solution was just to set the soil moisture content to 0 in the first place where it is used, if it turns into negative value, i.e., in subroutine calc_theta_mean(tile):

do k=1, kmax_soil
        do j=2,j1
            do i=2,i1
                
                
                
                if((phiw(i,j,k) <= 1e-5)) phiw(i,j,k)=0.          !if warmstart, 
    											!prognostic phiw(i,j,k) might turns into negative values for some grid points
    											!causing overflow and seg. fault, so  
    											!I set soil water content equals 0 if it is negative  !Arseni
    			    
                
                si = soil_index(i,j,k)

However, it needs to be checked more precisely.

2. theta_norm** (1. / vg_m) becomes greater than 1, resulting in a rising of negative value to a fraction power in pure function calc_diffusivity_vg and pure function calc_conductivity_vg: (1. - theta_norm ** (1. / vg_m)).

Here, I replace theta_norm** (1. / vg_m) with 0.9 everytime when theta_norm** (1. / vg_m)>1. However, it should be fixed in a more prompt way, I guess:

pure function calc_diffusivity_vg( &
        theta_norm, vg_a, vg_l, vg_m, lambda_sat, theta_sat, theta_res) result(res)
    implicit none
    real, intent(in) :: theta_norm, vg_a, vg_l, vg_m, lambda_sat, theta_sat, theta_res
    real :: res, res_old
    real :: part1, part2, part3, part4, part5

    
    
    part1 = (1. - vg_m)
    part2 = lambda_sat / (vg_a * vg_m * (theta_sat - theta_res))
    part3 = theta_norm ** (vg_l - (1. / vg_m))
    
    if((theta_norm ** (1. / vg_m))>=1.0) then  	!some values above 1 are appeared and the part4 and part5 turns into 
    		  							      	!a complex number (due to rising a negative value to a fraction power) that couses segmentation fault:
    											!here I put 0.9 to avoid negative values !Arseni
    											
    	part4 = (1. - 0.9) ** (-vg_m)
    	part5 = (1. - 0.9) ** vg_m - 2. 
    	
    else
    
    	part4 = (1. - theta_norm ** (1. / vg_m)) ** (-vg_m)
    	part5 = (1. - theta_norm ** (1. / vg_m)) ** vg_m - 2.
    endif
    
    
    

    res = part1 * part2 * part3 * (part4 + part5)
            
    
end function calc_diffusivity_vg

!
! Calculate hydraulic conductivity using van Genuchten parameterisation.
!
pure function calc_conductivity_vg(theta_norm, vg_l, vg_m, gamma_sat) result(res)
    implicit none
    real, intent(in) :: theta_norm, vg_l, vg_m, gamma_sat
    real :: res

	
    if((theta_norm ** (1. / vg_m))>=1.0) then  	!some values above 1 are appeared and it turns into 
    		  							      	!a complex number (due to rising a negative value to a fraction power) that couses segmentation fault:
    											!here I put 0.9 to avoid negative values !Arseni
    
    	res = gamma_sat * theta_norm**vg_l * ( 1.- (1.-0.9)**vg_m )**2.
    else
    
    	res = gamma_sat * theta_norm**vg_l * ( 1.- (1.-theta_norm**(1./vg_m))**vg_m )**2.
    	
    endif

end function calc_conductivity_vg

3. I also changed the subroutine Calc_canopy_resistance_ags so as not to calculate anything above water, I put an if statement (because over the water, tskin is NaN here and the model crashes); set output variables to 0 in case when Am < 0.01; Also, the units in co2_abs calculation have been corrected, co2_abs seems uses mg/m3? units here (Ronda et al., 2001), because the CO2 scalar tracers in DALES are expressed in units of µg/g, to correctly calculate co2_abs you simply need to multiply the value of co2 in µg/g by density at bottom layer ( rhof(1), kg/m3):

  to_abs   =  rhof(1) !from ug/g to mg/m3
    from_abs = 1/to_abs !from mg/m3 to ug/g

    do j=2,j1
        do i=2,i1
       	 	

       	 	! Check if tskin(i, j) is NaN
        	if (isnan(tskin(i, j))) then
            	
            	! Surface resistances for moisture and carbon dioxide:
           		tile_lv%rs(i,j) = 0.
            	tile_hv%rs(i,j) = 0.
            	tile_bs%rs(i,j) = 0.
            	tile_ws%rs(i,j) = 0.
            	
            	!respiration and vegetation:
            	resp_co2(i,j) = 0.
            	an_co2(i,j) = 0.        	
            	        	
            	
       	 	else
       	 	
       	 	    ! Check if tskin(i, j) less than 200 or more than 350, and replace with 287.
        		!if ((tskin(i, j).lt.200.).OR.(tskin(i, j).gt.350.)) then
            		!tskin(i, j) = 287.
       	 		!endif
    			
 				!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            
            	si = soil_index(i, j, kmax_soil)
            

			
				Ts = tskin(i,j) * exnh(1) 
			

            	! Calculate the CO2 compensation concentration (IFS eq. 8.92)
            	! "The compensation point Γ is defined as the CO2 concentration at which the net CO2 assimilation of a fully lit leaf becomes zero."
            	! NOTE: The old DALES LSM used the atmospheric `thl`, IFS uses the surface temperature.
            	co2_comp = rhof(1) * co2_comp298 * Q10lambda**(0.1 * (Ts - 298.0)) ![g m-3 or ppm??]

           	 	! Calculate the mesophyll conductance (IFS eq. 8.93)
            	! "The mesophyll conductance gm describes the transport of CO2 from the substomatal cavities to the mesophyll cells where the carbon is fixed."
            	! NOTE: The old DALES LSM used the atmospheric `thl`, IFS uses the surface temperature.
            	gm = gm298 * Q10gm**(0.1 * (Ts - 298.0)) / ((1. + exp(0.3 * (T1gm - Ts))) * (1. + exp(0.3 * (Ts - T2gm)))) / 1000.
	
            	! Calculate CO2 concentration inside the leaf (ci)
            	! NOTE: Differs from IFS
            	fmin0 = gmin/nuco2q - (1./9.)*gm
            	fmin  = (-fmin0 + (fmin0**2 + 4*gmin/nuco2q*gm)**0.5) / (2.*gm)

            	! Calculate atmospheric moisture deficit
           	 	! NOTE: "Therefore Ci/Cs is specified as a function of atmospheric moisture deficit Ds at the leaf surface"
            	!       In DALES, this uses a mix between esat(surface) and e(atmosphere)
            	!       In IFS, this uses (in kg kg-1): qsat(Ts)-qs instead of qsat(Ts)-qa!
            	! NOTE: Old DALES LSM used the surface pressure in the calculation of `e`, not sure why...
            	esatsurf = 0.611e3 * exp(17.2694 * (Ts - 273.16) / (Ts - 35.86))
            	e = qt0(i,j,1) * presf(1) / 0.622
            	Ds = (esatsurf - e) / 1000.
             
            
            	! This seems to differ from IFS?
            	Dmax = (f0 - fmin) / ad

            	! Coupling factor (IFS eq. 8.101)
            	!cfrac = f0 * (1.0 - Ds/Dmax) + fmin * (Ds/Dmax)
            	cfrac = max(0.01, f0 * (1.0 - Ds/Dmax) + fmin * (Ds/Dmax))

            	! Absolute CO2 concentration [mg/m3?]:
            	if (l_emission) then
            		co2_abs =  svm(i,j,1,svco2sum)*to_abs !from ug/g to mg/m3 (svm scalar tasers for gases are in ug/g now)
			    	
            	else
                	co2_abs =  svm(i,j,1,co2_index)*to_abs !from ug/g to mg/m3 (svm scalar tasers for gases are in ug/g now)

				endif
				

            !if (lrelaxci) then
            !  if (ci_old_set) then
            !    ci_inf        = cfrac * (co2_abs - co2_comp) + co2_comp
            !    ci            = ci_old(i,j) + min(kci*rk3coef, 1.0) * (ci_inf - ci_old(i,j))
            !    if (rk3step  == 3) then
            !      ci_old(i,j) = ci
            !    endif
            !  else
            !    ci            = cfrac * (co2_abs - co2_comp) + co2_comp
            !    ci_old(i,j)   = ci
            !  endif
            !else
            !  ci              = cfrac * (co2_abs - co2_comp) + co2_comp
            !endif

            	! CO2 concentration in leaf (IFS eq. ???): [mg m−3 or ppm]
            	ci = cfrac * (co2_abs - co2_comp) + co2_comp

            	! Max gross primary production in high light conditions (Ag) (IFS eq. 8.94)
            	! NOTE: The old DALES LSM used the atmospheric `thl`, IFS uses the surface temperature.
            	Ammax = Ammax298 * Q10am**(0.1 * (Ts - 298.0)) / ((1.0 + exp(0.3 * (T1Am - Ts))) * (1. + exp(0.3 * (Ts - T2Am))))

            	! Effect of soil moisture stress on gross assimilation rate.
           	 	! NOTE: this seems to be different in IFS...
            	! NOTE: for now, this uses the relative soil moisture content from the low vegetation only!
            	fstr = max(1.0e-3, min(1.0, tile_lv%phiw_mean(i,j)))
            

            	! Gross assimilation rate (Am, IFS eq. 8.97)
            	Am = Ammax * (1 - exp(-(gm * (ci - co2_comp) / Ammax)))
            
            	if((Am.lt.0.01).OR.(isnan(Am))) then	!I fix Am.lt.0.01 as threshold 
            											!since Am below might cause problems later in the code !Arseni
            
            
            		! Surface resistances for moisture and carbon dioxide:
           			tile_lv%rs(i,j) = 0.
            		tile_hv%rs(i,j) = 0.
            		tile_bs%rs(i,j) = 0.
            		tile_ws%rs(i,j) = 0.
            	
            		!respiration and vegetation:
            		resp_co2(i,j) = 0.
            		an_co2(i,j) = 0.
            
           	 	else
            

            		! Autotrophic dark respiration (IFS eq. 8.99)
            		Rdark = Am/9.
            
            
.................
            		! NOTE:Combined assimilation for low and high veg, using only low veg as vegetation type:
            		resp_co2(i,j) = R10 * (1.-fw) * exp(Eact0 / (283.15 * 8.314) * (1.0 - 283.15 / t_mean))
            		! Scale with vegetation fraction, and translate to `ug/g m s-1`.
			
					resp_co2(i,j) = cveg * resp_co2(i,j) * from_abs

            		! Calculate net assimilation
            		! NOTE: this uses the aerodynamic resistance from low vegetation...
            		an_co2(i,j) = -(co2_abs - ci) / (tile_lv%ra(i,j) + rs_co2)
            		! Scale with vegetation + bare soil fraction, and translate to `ug/g m s-1`.
            		an_co2(i,j) = cland * an_co2(i,j) * from_abs

                        !Replace unreasonable values with 0:
            		if (an_co2(i,j)>0.) an_co2(i,j)=0.
		        if (resp_co2(i,j)<0.) resp_co2(i,j)=0.

			
				endif
			endif
			
			
			
			! Set CO2 fluxes:
			
            if (l_emission) then
                ! Total flux into the CO2 field holding the sum of all CO2 concentrations:
                svflux(i,j,svco2sum) = resp_co2(i,j) + an_co2(i,j) !Kinematic scalar flux [ug/g m/s]
				
                ! Respiration flux into the respiration specific field:
                svflux(i,j,svco2ags) = resp_co2(i,j)  !Kinematic scalar flux [ug/g m/s]

                ! Net update into the vegetation specific field:
                !svflux(i,j,svco2veg) = an_co2(i,j) !Kinematic scalar flux [ug/g m/s]
            else
            
            	! Respiration flux into the respiration specific field:
                svflux(i,j,svco2ags) = resp_co2(i,j)  !Kinematic scalar flux [ug/g m/s]

                svflux(i,j,co2_index) = resp_co2(i,j) + an_co2(i,j)  !Kinematic scalar flux [ug/g m/s] 
            endif

        end do
    end do

Let me know, if I was not right somewhere and if you have a more smart way to solve issues I mentioned.

Cheers,
Arseni

Edits by Fredrik: added code tags

[jchylik]

1. grid points may turn into highly negative values

That is strange. For what kind of grids is it happening?
It has been a while since I used the surface model the last time, but I am surprised it does things like this. Are you sure the file was not corrupted?

[adoyenne]

1. grid points may turn into highly negative values

That is strange. For what kind of grids is it happening? It has been a while since I used the surface model the last time, but I am surprised it does things like this. Are you sure the file was not corrupted?

I found that if simulation is warm start (starts from restart files), phiw(i,j,k) might turns into negative values for some grid points. This problem appears not at the beginning of the simulation but after about 10 hours, so it is probably resulted from some calculating issues before the crash. I used a large domain in that simulation itot = 1536, jtot = 768, dx = 156.25, dy = 156.25, for smaller domains, I did not see the same problem. Restart files are not corrupted and phiw(i,j,k) is correct at the simulation start. I did not dig deeper, just put the if statement and substitute all negative values to reasonable ones. May problems with advection somehow cause wrong values here? I found that advection scheme may cause some problems with CO2 values (too low at some grid points) that why I shifted to kappa advection scheme and problem was solved.

[fjansson]

1. phiw should obey 0 < phiw < 1, otherwise it's not physical and there will be numerical problems in the LSM in several places. I've added a range check (at least temporarily) on phiw in the time step in modlsm.f90, commit 86feea1

2. theta_norm > 1 can happen if phiw > theta_sat (i.e. even if phiw is within the general limit 0 < phiw < 1)
   I wonder if we can limit also from above, by repacing

theta_lim = max(phiw(i,j,k), 1.001*theta_res(si))

with

theta_lim = min(max(phiw(i,j,k), 1.001*theta_res(si)), 0.999*theta_sat(si))

here
@julietbravo what do you think?

I'm currently getting crashes near the surface caused by a too high temperature. This may or may not be related to these LSM issues.

Point 3 is a separate issue in my opinion. From the point of view of the DALES core, scalars such as CO2 are just advected around. You can give them any unit you want, but it has to be consistent in the different schemes that use the same scalar.

[bartvstratum]

Yes, I think that limiting theta as in your point 2. is fine. I think you can even safely limit it at theta_sat instead of 0.999 theta_sat.

[bartvstratum]

@jchylik I only found this issue now, so my response is a bit late. We added a new land-surface model in the last few years, which more closely follows (C)HTESSEL, and allows for realistic land-use. I guess that you were still using the old one?

[adoyenne]

@fjansson thanks for this suggestion, I will try to implement it.

About issue 3, what I found is that the model crashed with seg fault right at the first-time step before calculation of tile_xx%rs. As I found, It is resulted from NaNs in tskin. Interestingly, the crash doesn't occur when I run the debug compiled code (cmake .. -DCMAKE_BUILD_TYPE="DEBUG") only at normal and fast run NaNs are appeared. So, tskin is somehow sensitive to the compilation type. In addition, tskin is not sensitive to the isnan function (isnan(tskin(i,j))), so in my code, I use this approach to exclude NaNs value (I printed tskin here and the value of 287K seems an average value of tskin ):

(I used a copy to not change the original tskin)

! Check if tskin_copy (i, j) is less than 200 or more than 350, and replace with 287.
if ((tskin_copy (i, j) < 200.) .OR. (tskin_copy (i, j) > 350.)) then
tskin_copy (i, j) = 287.
Endif

I do not know, if NaNs in tskin are appeared everywhere or only in the subroutine Calc_canopy_resistance_ags, good to check.

About units of CO2 tracer, here could be some inconsistency with others, since Marko implemented GHG scalar tracers in ug/g unit, whereas in the A-Gs scheme, it is seem assumed that the unit of co2 input is in ppb.