From cc4a1d05b3c29df9325b2bf2f343b1764d43f4f9 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Tue, 27 Jun 2023 14:13:05 -0400 Subject: [PATCH 1/4] add u,v to vertical remapping --- src/soca/Fields/soca_fields_mod.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/soca/Fields/soca_fields_mod.F90 b/src/soca/Fields/soca_fields_mod.F90 index 98f8e42f0..d1790009d 100644 --- a/src/soca/Fields/soca_fields_mod.F90 +++ b/src/soca/Fields/soca_fields_mod.F90 @@ -361,7 +361,7 @@ subroutine soca_field_stencil_interp(self, geom, fromto) ! source point on land, skip if (masksrc_local(ij(1,sti), ij(2,sti)) == 0_kind_real) cycle - ! outcroping of layers, skip + ! outcroping of layers, skip if (abs(self%val(ij(1,sti), ij(2,sti),1)) > val_max) cycle ! store the valid neighbors @@ -1005,8 +1005,8 @@ subroutine soca_fields_read(self, f_conf, vdate) field => self%fields(n) select case(field%name) ! TODO remove hardcoded variable names here - ! TODO Add u and v. Remapping u and v will require interpolating h - case ('tocn','socn') + ! TODO Remapping u and v is only valid if they are on the tracer grid point. + case ('tocn','socn','uocn','vocn') if (associated(field%mask) .and. field%mask(i,j).eq.1) then varocn_ij = field%val(i,j,:) call remapping_core_h(remapCS, nz, h_common_ij, varocn_ij,& From 8cf2b6a4d27a9b8287a86b596d74a03637d2ac2b Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Tue, 27 Jun 2023 17:24:58 -0400 Subject: [PATCH 2/4] added biogeochem var to the vertical remapping --- src/soca/Fields/soca_fields_mod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/soca/Fields/soca_fields_mod.F90 b/src/soca/Fields/soca_fields_mod.F90 index d1790009d..61fc3ca44 100644 --- a/src/soca/Fields/soca_fields_mod.F90 +++ b/src/soca/Fields/soca_fields_mod.F90 @@ -1006,7 +1006,7 @@ subroutine soca_fields_read(self, f_conf, vdate) select case(field%name) ! TODO remove hardcoded variable names here ! TODO Remapping u and v is only valid if they are on the tracer grid point. - case ('tocn','socn','uocn','vocn') + case ('tocn','socn','uocn','vocn','chl','biop','po4') if (associated(field%mask) .and. field%mask(i,j).eq.1) then varocn_ij = field%val(i,j,:) call remapping_core_h(remapCS, nz, h_common_ij, varocn_ij,& From ed7bd11c583b5cdcce03b51ceafae569996ec767 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Thu, 6 Jul 2023 11:05:16 -0400 Subject: [PATCH 3/4] added vert_interp to fld metadata --- src/soca/Fields/soca_fields_metadata_mod.F90 | 10 ++++ src/soca/Fields/soca_fields_mod.F90 | 61 ++++++++++---------- test/Data/fields_metadata.yml | 2 + 3 files changed, 41 insertions(+), 32 deletions(-) diff --git a/src/soca/Fields/soca_fields_metadata_mod.F90 b/src/soca/Fields/soca_fields_metadata_mod.F90 index fd317241c..c6dc793e9 100644 --- a/src/soca/Fields/soca_fields_metadata_mod.F90 +++ b/src/soca/Fields/soca_fields_metadata_mod.F90 @@ -27,6 +27,7 @@ module soca_fields_metadata_mod character(len=:), allocatable :: io_file !< the restart file domain (ocn, sfc, ice) character(len=:), allocatable :: io_name !< the name use in the restart IO character(len=:), allocatable :: property !< physical property of the field, "none" or "positive_definite" + logical :: vert_interp !< true if the field can be vertically interpolated end type @@ -106,6 +107,15 @@ subroutine soca_fields_metadata_create(self, filename) if(.not. conf_list(i)%get("property", str)) str = "none" self%metadata(i)%property = str + if(.not. conf_list(i)%get("vert interp", bool)) then + if (self%metadata(i)%levels == "1" ) then + bool = .false. + else + bool = .true. + end if + end if + self%metadata(i)%vert_interp = bool + end do ! check for duplicates diff --git a/src/soca/Fields/soca_fields_mod.F90 b/src/soca/Fields/soca_fields_mod.F90 index 61fc3ca44..8eac072b5 100644 --- a/src/soca/Fields/soca_fields_mod.F90 +++ b/src/soca/Fields/soca_fields_mod.F90 @@ -969,6 +969,35 @@ subroutine soca_fields_read(self, f_conf, vdate) isc = self%geom%isc ; iec = self%geom%iec jsc = self%geom%jsc ; jec = self%geom%jec + ! Remap layers if needed + if (vert_remap) then + allocate(h_common_ij(nz), hocn_ij(nz), varocn_ij(nz), varocn2_ij(nz)) + call initialize_remapping(remapCS,'PCM') + do i = isc, iec + do j = jsc, jec + h_common_ij = h_common(i,j,:) + hocn_ij = hocn%val(i,j,:) + + do n=1,size(self%fields) + field => self%fields(n) + ! TODO Vertical remapping is only valid if the field is on the tracer grid point. + if (.not. field%metadata%vert_interp) cycle + if (associated(field%mask) .and. field%mask(i,j).eq.1) then + varocn_ij = field%val(i,j,:) + call remapping_core_h(remapCS, nz, h_common_ij, varocn_ij,& + &nz, hocn_ij, varocn2_ij) + field%val(i,j,:) = varocn2_ij + else + field%val(i,j,:) = 0.0_kind_real + end if + end do + end do + end do + hocn%val = h_common + deallocate(h_common_ij, hocn_ij, varocn_ij, varocn2_ij) + call end_remapping(remapCS) + end if + ! Initialize mid-layer depth from layer thickness if (self%has("layer_depth")) then call self%get("layer_depth", layer_depth) @@ -992,38 +1021,6 @@ subroutine soca_fields_read(self, f_conf, vdate) end do end if - ! Remap layers if needed - if (vert_remap) then - allocate(h_common_ij(nz), hocn_ij(nz), varocn_ij(nz), varocn2_ij(nz)) - call initialize_remapping(remapCS,'PCM') - do i = isc, iec - do j = jsc, jec - h_common_ij = h_common(i,j,:) - hocn_ij = hocn%val(i,j,:) - - do n=1,size(self%fields) - field => self%fields(n) - select case(field%name) - ! TODO remove hardcoded variable names here - ! TODO Remapping u and v is only valid if they are on the tracer grid point. - case ('tocn','socn','uocn','vocn','chl','biop','po4') - if (associated(field%mask) .and. field%mask(i,j).eq.1) then - varocn_ij = field%val(i,j,:) - call remapping_core_h(remapCS, nz, h_common_ij, varocn_ij,& - &nz, hocn_ij, varocn2_ij) - field%val(i,j,:) = varocn2_ij - else - field%val(i,j,:) = 0.0_kind_real - end if - end select - end do - end do - end do - hocn%val = h_common - deallocate(h_common_ij, hocn_ij, varocn_ij, varocn2_ij) - call end_remapping(remapCS) - end if - ! Update halo do n=1,size(self%fields) field => self%fields(n) diff --git a/test/Data/fields_metadata.yml b/test/Data/fields_metadata.yml index 371bed01d..e2807aeab 100644 --- a/test/Data/fields_metadata.yml +++ b/test/Data/fields_metadata.yml @@ -51,6 +51,7 @@ getval name: sea_water_cell_thickness io file: ocn io name: h + vert interp: false - name: ssh getval name: sea_surface_height_above_geoid @@ -149,6 +150,7 @@ - name: layer_depth levels: full_ocn + vert interp: false - name: mesoscale_representation_error From 2cbc7236f48689cc9b028569c2085afe13532875 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Thu, 6 Jul 2023 16:10:26 -0400 Subject: [PATCH 4/4] output log of interpolated fields --- src/soca/Fields/soca_fields_mod.F90 | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/soca/Fields/soca_fields_mod.F90 b/src/soca/Fields/soca_fields_mod.F90 index 8eac072b5..49218cbbf 100644 --- a/src/soca/Fields/soca_fields_mod.F90 +++ b/src/soca/Fields/soca_fields_mod.F90 @@ -18,6 +18,7 @@ module soca_fields_mod datetime_create, datetime_diff use duration_mod, only: duration, duration_to_string use fckit_configuration_module, only: fckit_configuration +use fckit_log_module, only: fckit_log use fckit_mpi_module, only: fckit_mpi_min, fckit_mpi_max, fckit_mpi_sum use kinds, only: kind_real use oops_variables_mod, only: oops_variables @@ -971,6 +972,13 @@ subroutine soca_fields_read(self, f_conf, vdate) ! Remap layers if needed if (vert_remap) then + + ! output log of what fields are going to be interpolated vertically + do n=1,size(self%fields) + if (.not. self%fields(n)%metadata%vert_interp) cycle + call fckit_log%info("vertically remapping "//trim(self%fields(n)%name)) + end do + allocate(h_common_ij(nz), hocn_ij(nz), varocn_ij(nz), varocn2_ij(nz)) call initialize_remapping(remapCS,'PCM') do i = isc, iec