Support mapping of mixing ratio profiles (ppb to ppb). Take care of undefined... (ea9c8e3f) · Commits · CAMS / CSO

oper/src/cso_pixels.F90

+142 −38

Original line number	Diff line number	Diff line
		@@ -33,6 +33,9 @@
		! Support kernel convolution 'xa + A ( x - xa )'.
		!
		! 2026-07, Arjo Segers
		! Support mapping of mixing ratio profiles (ppb to ppb).
		! Take care of undefined layers.
		! Support `ya + A ( x - xa )` formula for kernel (profile) convolution.
		! Updated logging messages.
		!
		!###############################################################################
		@@ -4677,12 +4680,15 @@ contains
		real :: fill_value

		integer :: nlayer
		integer :: ilayer
		integer :: nz
		real, pointer :: hp_data(:,:) ! (nlayer,npix)
		character(len=64) :: hp_units
		integer :: hp_isfc
		integer :: hp_iblh
		integer :: hp_itop
		real :: hp_fill_value
		real :: hp_sfc
		real, pointer :: mod_hp_data(:,:) ! (nz+1,npix)
		character(len=64) :: mod_hp_units
		integer :: mod_hp_isfc
		@@ -4695,6 +4701,7 @@ contains
		type(T_ProfileMapping) :: ProfileMapping
		real, allocatable :: mod_hpx(:) ! (0:nz)
		real, allocatable :: mod_g(:) ! (nz)
		real, allocatable :: hpx(:) ! (0:nlayer)
		real, allocatable :: hx(:) ! (nlayer)

		real, pointer :: A_data(:,:,:) ! (nretr,nlayer,npix)
		@@ -4715,6 +4722,8 @@ contains
		character(len=64) :: nla_units
		real, pointer :: y_data(:,:) ! (*,npix)
		character(len=64) :: y_units
		real, pointer :: ya_data(:,:) ! (*,npix)
		character(len=64) :: ya_units
		real, pointer :: y_m_data(:,:) ! (nlayer,npix)
		character(len=64) :: y_m_units
		real, pointer :: R_data(:,:,:) ! (nretr,nretr,npix)
		@@ -4756,7 +4765,7 @@ contains
		IF_NOT_OK_RETURN(status=1)
		! get pointers to source arrays:
		call self%GetFormulaData( p%formula_terms, 'hp', status, pd=pd, data1=hp_data, &
		units=hp_units, isfc_pressure=hp_isfc )
		units=hp_units, isfc_pressure=hp_isfc, fill_value=hp_fill_value )
		IF_NOT_OK_RETURN(status=1)
		call self%GetFormulaData( p%formula_terms, 'mod_hp', status, pd=pd, data1=mod_hp_data, &
		units=mod_hp_units, isfc_pressure=mod_hp_isfc )
		@@ -4783,19 +4792,22 @@ contains
		IF_NOT_OK_RETURN(status=1)

		! storage:
		allocate( mod_hpx(0:nz), stat=status )
		allocate( mod_hpx(1:nz+1), stat=status )
		IF_NOT_OK_RETURN(status=1)
		allocate( hx(nlayer), stat=status )
		allocate( hpx(1:nlayer+1), stat=status )
		IF_NOT_OK_RETURN(status=1)
		allocate( hx(1:nlayer), stat=status )
		IF_NOT_OK_RETURN(status=1)

		! select on unit conversion:
		select case ( trim(mod_conc_units)//' -> '//trim(y_units) )
		!~ mass mixing ratio to column density:

		!~ volume mixing ratio to column density:
		case ( 'ppb -> mol m-2' )
		! loop over pixels:
		do ipix = 1, npix
		! skip if no-data:
		if ( mod_hp_data(1,ipix) == p%fill_value ) cycle
		if ( all( mod_hp_data(:,ipix) == p%fill_value ) ) cycle
		! scale model pressures to have same surface pressure as pixel,
		! slightly lower surface to avoid tiny mismatches:
		mod_hpx = mod_hp_data(:,ipix) / mod_hp_data(mod_hp_isfc,ipix) * hp_data(hp_isfc,ipix)*1.0001
		@@ -4822,6 +4834,39 @@ contains
		/ xm_air & ! / ((kg air)/(mole air))
		* 1.0e-9 ! * (mole tr)/(mole air)/ppb
		end do ! ipix

		!~ volume mixing ratio "interpolation":
		case ( 'ppb -> ppb', 'ppmv -> ppmv' )
		! loop over pixels:
		do ipix = 1, npix
		! skip if no-data:
		if ( all( mod_hp_data(:,ipix) == p%fill_value ) ) cycle
		! skip if zero contirbution from this domain:
		if ( mod_hp_data(mod_hp_isfc,ipix) <= 0.0 ) cycle

		! largest defined layer in hp_data provides the surface pressure:
		hp_sfc = maxval( hp_data(:,ipix), mask=(hp_data(:,ipix) /= hp_fill_value) )
		! copy, fill surface for undefined layers:
		where ( hp_data(:,ipix) == hp_fill_value )
		hpx = hp_sfc
		elsewhere
		hpx = hp_data(:,ipix)
		endwhere

		!~ copy model pressures
		mod_hpx = mod_hp_data(:,ipix)
		!~ reset model surface to have model at least as low as target:
		mod_hpx(mod_hp_isfc) = max( hp_sfc, mod_hpx(mod_hp_isfc) )
		!~ weights:
		call ProfileMapping%Setup( mod_hpx, hpx, status )
		IF_NOT_OK_RETURN(status=1)
		! pressure (air mass) weighted mean:
		! y(j) = [ sum_i c(i) dp(i) ]/dp(j)
		! ppb Pa / Pa
		call ProfileMapping%Apply_WeightedMean( mod_conc_data(:,ipix), y_data(:,ipix), status )
		IF_NOT_OK_RETURN(status=1)
		end do ! ipix

		!~
		case default
		write (csol,'("unsupported conversion `",a,"`")') &
		@@ -4832,6 +4877,8 @@ contains
		! clear:
		deallocate( mod_hpx, stat=status )
		IF_NOT_OK_RETURN(status=1)
		deallocate( hpx, stat=status )
		IF_NOT_OK_RETURN(status=1)
		deallocate( hx, stat=status )
		IF_NOT_OK_RETURN(status=1)

		@@ -4996,7 +5043,7 @@ contains
		end do ! ipix

		!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		!~ kernel convolution:
		!~ kernel matrix convolution:
		case ( 'xa + A ( x - xa )' )
		!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		@@ -5035,7 +5082,7 @@ contains
		! apply:
		do ipix = 1, npix
		! filter on no-data ..
		if ( x_data(1,ipix) == fill_value ) cycle
		if ( all( x_data(:,ipix) == fill_value ) ) cycle
		! with some compilers problem using "matmul"; instead,
		! loop over target layers:
		do iretr = 1, size(y_data,1)
		@@ -5052,6 +5099,63 @@ contains
		end do ! iretr
		end do ! ipix

		!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		!~ kernel profile convolution:
		case ( 'ya + A ( x - xa )' )
		!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		! pointer to target array:
		call self%GetData( status, id=id, data1=y_data, units=y_units )
		IF_NOT_OK_RETURN(status=1)
		! get pointers to source arrays:
		call self%GetFormulaData( p%formula_terms, 'A', status, pd=pd, data2=A_data, units=A_units )
		IF_NOT_OK_RETURN(status=1)
		call self%GetFormulaData( p%formula_terms, 'xa', status, pd=pd, data1=xa_data, units=xa_units, fill_value=fill_value )
		IF_NOT_OK_RETURN(status=1)
		call self%GetFormulaData( p%formula_terms, 'x', status, pd=pd, data1=x_data, units=x_units, fill_value=fill_value )
		IF_NOT_OK_RETURN(status=1)
		call self%GetFormulaData( p%formula_terms, 'ya', status, pd=pd, data1=ya_data, units=ya_units, fill_value=fill_value )
		IF_NOT_OK_RETURN(status=1)
		! check units:
		if ( trim(A_units) /= '1' ) then
		write (csol,'("A units `",a,"` should be `1`")') trim(A_units); call csoErr
		write (csol,'(" formula : ",a)') trim(p%formula); call csoErr
		write (csol,'(" formula_terms : ",a)') trim(p%formula_terms); call csoErr
		write (csol,'(" variable : ",a)') trim(p%name); call csoErr
		TRACEBACK; status=1; return
		end if
		if ( trim(y_units) /= trim(ya_units) ) then
		write (csol,'("output units `",a,"` should be equal to ya units `",a,"`")') trim(y_units), trim(ya_units); call csoErr
		write (csol,'(" formula : ",a)') trim(p%formula); call csoErr
		write (csol,'(" formula_terms : ",a)') trim(p%formula_terms); call csoErr
		write (csol,'(" variable : ",a)') trim(p%name); call csoErr
		TRACEBACK; status=1; return
		end if
		if ( trim(x_units) /= trim(x_units) ) then
		write (csol,'("input units `",a,"` should be equal to xa units `",a,"`")') trim(x_units), trim(xa_units); call csoErr
		write (csol,'(" formula : ",a)') trim(p%formula); call csoErr
		write (csol,'(" formula_terms : ",a)') trim(p%formula_terms); call csoErr
		write (csol,'(" variable : ",a)') trim(p%name); call csoErr
		TRACEBACK; status=1; return
		end if
		! apply:
		do ipix = 1, npix
		! filter on no-data ..
		if ( all( x_data(:,ipix) == fill_value ) ) cycle
		! with some compilers problem using "matmul"; instead,
		! loop over target layers:
		do iretr = 1, size(y_data,1)
		! fill with dot product:
		do ilayer = 1, size(xa_data,1)
		! filter:
		if ( xa_data(ilayer,ipix) == fill_value ) cycle
		! add contribution:
		y_data(iretr,ipix) = ya_data(iretr,ipix) + &
		A_data(iretr,ilayer,ipix) * ( x_data(ilayer,ipix) - xa_data(ilayer,ipix) )
		end do ! ilayer
		end do ! iretr
		end do ! ipix

		!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		!~ kernel convolution, no apriori
		case ( 'A x' )