Variables sent to atmosphere are not averaged properly #1

Issue description

From the main module output fields send to atmosphere are added with

!...
 DO j=1,MAX_VARS
    myname = name_send_t(j)
    IF (trim(myname) /= 'none') THEN  
      IF (send_to_atmos_t(j)) THEN
        CALL add_output_field(myname, 'A', 0, 1, num_surface_types, grid_size(1), .TRUE., val_flux_t(j), &
                              local_field, num_input_fields, input_field, num_output_fields, output_field)
							 
!...

Setting uniform=.TRUE. leads in the add_output_field routine field to

!...
	                    IF (surface_type==0) THEN
                        ! Yes, output for the entire cell
                        IF (uniform) THEN
                            ! Flux is the same everywhere so we only need to find one surface_type for which it is calculated
                            DO j=1,num_surface_types
                                IF (ASSOCIATED(local_field(j,which_grid)%var(i)%field) .AND. (.NOT. ASSOCIATED(local_field(0,which_grid)%var(i)%field))) THEN
                                    local_field(0,which_grid)%var(i)%field => local_field(j,which_grid)%var(i)%field
!...

which means that local_field(0,which_grid)%var(i)%field is not allocated and just points to the other elements, thus no allocated flag is set for this variable. Since in the routine average_across_surface_types this flag is checked

!...
		IF (local_field(0,1)%var(my_idx)%allocated) THEN
            WRITE (w_unit,*) "Really averaging"
            local_field(0,1)%var(my_idx)%field=0.0
            DO i=1,num_surface_types
                DO j=1,grid_size(which_grid)
                    local_field(0,1)%var(my_idx)%field(j) = local_field(0,1)%var(my_idx)%field(j) + &
                                                            local_field(i,1)%var(my_idx)%field(j)*local_field(i,1)%var(idx_FARE)%field(j)
                ENDDO
            ENDDO
        ENDIF
!...

an averaging did never actually happen!

# Issue description From the main module output fields send to atmosphere are added with ``` fortran !... DO j=1,MAX_VARS myname = name_send_t(j) IF (trim(myname) /= 'none') THEN IF (send_to_atmos_t(j)) THEN CALL add_output_field(myname, 'A', 0, 1, num_surface_types, grid_size(1), .TRUE., val_flux_t(j), & local_field, num_input_fields, input_field, num_output_fields, output_field) !... ``` Setting `uniform=.TRUE.` leads in the add_output_field routine field to ``` fortran !... IF (surface_type==0) THEN ! Yes, output for the entire cell IF (uniform) THEN ! Flux is the same everywhere so we only need to find one surface_type for which it is calculated DO j=1,num_surface_types IF (ASSOCIATED(local_field(j,which_grid)%var(i)%field) .AND. (.NOT. ASSOCIATED(local_field(0,which_grid)%var(i)%field))) THEN local_field(0,which_grid)%var(i)%field => local_field(j,which_grid)%var(i)%field !... ``` which means that `local_field(0,which_grid)%var(i)%field` is not allocated and just points to the other elements, thus no `allocated` flag is set for this variable. Since in the routine `average_across_surface_types` this flag is checked ``` fortran !... IF (local_field(0,1)%var(my_idx)%allocated) THEN WRITE (w_unit,*) "Really averaging" local_field(0,1)%var(my_idx)%field=0.0 DO i=1,num_surface_types DO j=1,grid_size(which_grid) local_field(0,1)%var(my_idx)%field(j) = local_field(0,1)%var(my_idx)%field(j) + & local_field(i,1)%var(my_idx)%field(j)*local_field(i,1)%var(idx_FARE)%field(j) ENDDO ENDDO ENDIF !... ``` **an averaging did never actually happen!**

First idea

If one would call add_output_field(..., .FALSE., ...) with uniform=.FALSE. one could achieve

!...
                            ! Flux will differ, so we will need a value for each surface_type and an area fraction of that surface_type
                            found_all_fluxes = .TRUE.
                            found_all_areas  = .TRUE.
                            DO j=1,num_surface_types
                                IF (.NOT. ASSOCIATED(local_field(j,which_grid)%var(i)%field))        found_all_fluxes=.FALSE.
                                IF (.NOT. ASSOCIATED(local_field(j,which_grid)%var(idx_FARE)%field)) found_all_areas=.FALSE.
                            ENDDO
                            IF (.NOT. found_all_areas) THEN
                                ! fractional area is not provided - this is fatal
                                WRITE (w_unit,*) "ERROR: Output field ",myname," has not been defined as uniform (flux_?_uniform=.FALSE.)."
                                WRITE (w_unit,*) "To calculate its average value across different surface_types, their fractional area (FARE) must be given but is missing."
                                CALL mpi_finalize(1)
                            ELSE 
                                IF (found_all_fluxes) THEN
                                    ! this is nice, we found everything we need, so just allocate
                                    IF (.NOT. ASSOCIATED(local_field(0,which_grid)%var(i)%field)) THEN
                                        ALLOCATE(local_field(0,which_grid)%var(i)%field(grid_size))
                                        local_field(0,which_grid)%var(i)%allocated=.TRUE.
                                    ENDIF
                                ENDIF
                            ENDIF
!...

i.e. the allocated flag is set. However, only if ASSOCIATED(local_field(j,which_grid)%var(idx_FARE)%field) is set, which is not true for the u- and v-grid.

Here it could be checked if only ASSOCIATED(local_field(j,1)%var(idx_FARE)%field) (for t-grid). However, since average_across_surface_types only checks for local_field(0,1)%var(my_idx)%allocated this still doesn't help u- and v-grid variables.

Probably the average_across_surface_types and add_output_field have to be strikingly modified.

# First idea If one would call `add_output_field(..., .FALSE., ...)` with `uniform=.FALSE.` one could achieve ``` fortran !... ! Flux will differ, so we will need a value for each surface_type and an area fraction of that surface_type found_all_fluxes = .TRUE. found_all_areas = .TRUE. DO j=1,num_surface_types IF (.NOT. ASSOCIATED(local_field(j,which_grid)%var(i)%field)) found_all_fluxes=.FALSE. IF (.NOT. ASSOCIATED(local_field(j,which_grid)%var(idx_FARE)%field)) found_all_areas=.FALSE. ENDDO IF (.NOT. found_all_areas) THEN ! fractional area is not provided - this is fatal WRITE (w_unit,*) "ERROR: Output field ",myname," has not been defined as uniform (flux_?_uniform=.FALSE.)." WRITE (w_unit,*) "To calculate its average value across different surface_types, their fractional area (FARE) must be given but is missing." CALL mpi_finalize(1) ELSE IF (found_all_fluxes) THEN ! this is nice, we found everything we need, so just allocate IF (.NOT. ASSOCIATED(local_field(0,which_grid)%var(i)%field)) THEN ALLOCATE(local_field(0,which_grid)%var(i)%field(grid_size)) local_field(0,which_grid)%var(i)%allocated=.TRUE. ENDIF ENDIF ENDIF !... ``` i.e. the allocated flag is set. However, only if `ASSOCIATED(local_field(j,which_grid)%var(idx_FARE)%field)` is set, which is not true for the u- and v-grid. Here it could be checked if only `ASSOCIATED(local_field(j,1)%var(idx_FARE)%field)` (for t-grid). However, since `average_across_surface_types` only checks for `local_field(0,1)%var(my_idx)%allocated` this still doesn't help u- and v-grid variables. **Probably the `average_across_surface_types` and `add_output_field` have to be strikingly modified.**

Not too much needs to be changed, as far as I see.

"First idea" should be correct: Calling add_output_field with uniform=.FALSE. Even better would be to check the existing namelist option send_uniform_t(i) and provide either .TRUE. or .FALSE. depending on that.

The code shown under "First idea" should then, at least for the t grid, do the correct thing and allocate the field for surface type zero, rather than just setting a pointer.

This allocation then should already lead to a correct averaging, as the code in the original issue shows:

!...
	IF (local_field(0,1)%var(my_idx)%allocated) THEN        
    	WRITE (w_unit,*) "Really averaging"            
		local_field(0,1)%var(my_idx)%field=0.0            
		DO i=1,num_surface_types            
			DO j=1,grid_size(which_grid)                
				local_field(0,1)%var(my_idx)%field(j) = local_field(0,1)%var(my_idx)%field(j) + &                                                     			local_field(i,1)%var(my_idx)%field(j)*local_field(i,1)%var(idx_FARE)%field(j)                 
			ENDDO
		ENDDO
	ENDIF
!...

This averaging is only possible for the t-grid, since (at least for the MOM ocean model) the fractional area FARE only exists for the t-grid. All variables on the u-grid and v-grid which are sent to the atmosmodel must have send_uniform_u(i)=.TRUE. and send_uniform_v(i)=.TRUE.

This averaging is only possible for the t-grid, since (at least for the MOM ocean model) the fractional area FARE only exists for the t-grid. All variables on the u-grid and v-grid which are sent to the atmosmodel must have send_uniform_u(i)=.TRUE. and send_uniform_v(i)=.TRUE.

But what about the UMOM00 and VMOM00 fluxes passed to the atmosphere? Shouldn't they be the averages of UMOM01..06 and VMOM01..06, respectively?

Even better would be to check the existing namelist option send_uniform_t(i) and provide either .TRUE. or .FALSE. depending on that.

Unfortuantely, there is not send_uniform_t(i). There is only a send_uniform_t(i,j) where the first index corresponds to the bottom model. So this quantity is not yet a property of the variable sent to the atmosphere. I guess the original idea was that variables sent to the atmosphere are anyway always uniform. However, they might not be always uniform but rather always averaged?

Indeed UMOM and VMOM need to be averaged. So probably I was wrong and we will need a regridding of FARE from the t-grid to the u-grid and v-grid to enable this averaging.

It is correct that the send_uniform_t property can have different values for different bottom models. But since every instance of the flux calculator is responsible for only one bottom model, this should be no problem.