Opened 7 weeks ago

Closed 10 days ago

#2665 closed Bug (fixed)

Various fixes required to run with key_qco or key_linssh

Reported by: hadcv Owned by: hadcv
Priority: normal Milestone: Unscheduled
Component: MULTIPLE Version: trunk
Severity: minor Keywords: QCO, key_qco, key_linssh
Cc: techene Branch review:
MP ready?: yes Task progress:

Description

Context

Running certain code with key_qco or key_linssh in ORCA2 (and possibly other configurations) results in segfaults.

Analysis

Workarounds and/or #include domzgr_substitute.h90 directives are required.

Fix

#include domzgr_substitute.h90 directives need to be added to diaobs.F90, obs_prep.F90, tradmp.F90, and zdfric.F90.

Additionally, in diaobs.F90:

   SUBROUTINE dia_obs( kstp, Kmm )
      !!----------------------------------------------------------------------
      !!                    ***  ROUTINE dia_obs  ***
      !!          
      !! ** Purpose : Call the observation operators on each time step
      !!
      !! ** Method  : Call the observation operators on each time step to
      !!              compute the model equivalent of the following data:
      !!               - Profile data, currently T/S or U/V
      !!               - Surface data, currently SST, SLA or sea-ice concentration.
      !!
      !! ** Action  :
      !!----------------------------------------------------------------------
      USE dom_oce, ONLY : gdept, gdept_1d     ! Ocean space domain variables (Kmm time-level only)
      USE phycst , ONLY : rday                ! Physical constants
      USE oce    , ONLY : ts, uu, vv, ssh     ! Ocean dynamics and tracers variables (Kmm time-level only)
      USE phycst , ONLY : rday                ! Physical constants
#if defined  key_si3
      USE ice    , ONLY : at_i                ! SI3 Ice model variables
#endif
#if defined key_cice
      USE sbc_oce, ONLY : fr_i     ! ice fraction
#endif

      IMPLICIT NONE

      !! * Arguments
      INTEGER, INTENT(IN) :: kstp  ! Current timestep
      INTEGER, INTENT(in) :: Kmm   ! ocean time level indices
      !! * Local declarations
      INTEGER :: idaystp           ! Number of timesteps per day
      INTEGER :: jtype             ! Data loop variable
      INTEGER :: jvar              ! Variable number
-     INTEGER :: ji, jj            ! Loop counters
+     INTEGER :: ji, jj, jk        ! Loop counters
      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
         & zprofvar                ! Model values for variables in a prof ob
      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
         & zprofmask               ! Mask associated with zprofvar
      REAL(wp), DIMENSION(jpi,jpj) :: &
         & zsurfvar, &             ! Model values equivalent to surface ob.
         & zsurfmask               ! Mask associated with surface variable
      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
         & zglam,    &             ! Model longitudes for prof variables
         & zgphi                   ! Model latitudes for prof variables
+     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zdept, zdepw

      !-----------------------------------------------------------------------

      IF(lwp) THEN
         WRITE(numout,*)
         WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
         WRITE(numout,*) '~~~~~~~'
      ENDIF

      idaystp = NINT( rday / rn_Dt )

      !-----------------------------------------------------------------------
      ! Call the profile and surface observation operators
      !-----------------------------------------------------------------------

      IF ( nproftypes > 0 ) THEN

+        ALLOCATE( zdept(jpi,jpj,jpk), zdepw(jpi,jpj,jpk) )
+        DO jk = 1, jpk
+           zdept(:,:,jk) = gdept(:,:,jk,Kmm)
+           zdepw(:,:,jk) = gdepw(:,:,jk,Kmm)
+        END DO

         DO jtype = 1, nproftypes

            ! Allocate local work arrays
            ALLOCATE( zprofvar (jpi, jpj, jpk, profdataqc(jtype)%nvar) )
            ALLOCATE( zprofmask(jpi, jpj, jpk, profdataqc(jtype)%nvar) )
            ALLOCATE( zglam    (jpi, jpj,      profdataqc(jtype)%nvar) )
            ALLOCATE( zgphi    (jpi, jpj,      profdataqc(jtype)%nvar) )  
                              
            ! Defaults which might change
            DO jvar = 1, profdataqc(jtype)%nvar
               zprofmask(:,:,:,jvar) = tmask(:,:,:)
               zglam(:,:,jvar)       = glamt(:,:)
               zgphi(:,:,jvar)       = gphit(:,:)
            END DO

            SELECT CASE ( TRIM(cobstypesprof(jtype)) )
            CASE('prof')
               zprofvar(:,:,:,1) = ts(:,:,:,jp_tem,Kmm)
               zprofvar(:,:,:,2) = ts(:,:,:,jp_sal,Kmm)
            CASE('vel')
               zprofvar(:,:,:,1) = uu(:,:,:,Kmm)
               zprofvar(:,:,:,2) = vv(:,:,:,Kmm)
               zprofmask(:,:,:,1) = umask(:,:,:)
               zprofmask(:,:,:,2) = vmask(:,:,:)
               zglam(:,:,1) = glamu(:,:)
               zglam(:,:,2) = glamv(:,:)
               zgphi(:,:,1) = gphiu(:,:)
               zgphi(:,:,2) = gphiv(:,:)
            CASE DEFAULT
               CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' )
            END SELECT

            DO jvar = 1, profdataqc(jtype)%nvar
               CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  &
                  &               nit000, idaystp, jvar,                   &
                  &               zprofvar(:,:,:,jvar),                    &
-                 &               gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm),      &
+                 &               zdept(:,:,:), zdepw(:,:,:),      &
                  &               zprofmask(:,:,:,jvar),                   &
                  &               zglam(:,:,jvar), zgphi(:,:,jvar),        &
                  &               nn_1dint, nn_2dint_default,              &
                  &               kdailyavtypes = nn_profdavtypes )
            END DO
            
            DEALLOCATE( zprofvar, zprofmask, zglam, zgphi )

         END DO

+        DEALLOCATE( zdept, zdepw )
+
      ENDIF

and in obs_prep.F90:

   SUBROUTINE obs_coo_spc_3d( kprofno, kobsno,  kpstart, kpend, &
      &                       kpi,     kpj,     kpk,            &
      &                       kobsi,   kobsj,   kobsk,          &
      &                       pobslam, pobsphi, pobsdep,        &
      &                       plam,    pphi,    pdep,    pmask, &
      &                       kpobsqc, kobsqc,  kosdobs,        &
      &                       klanobs, knlaobs, ld_nea,         &
      &                       kbdyobs, ld_bound_reject,         &
      &                       kqc_cutoff,       Kmm             )
      !!----------------------------------------------------------------------
      !!                    ***  ROUTINE obs_coo_spc_3d  ***
      !!
      !! ** Purpose : Check for points outside the domain and land points
      !!              Reset depth of observation above highest model level
      !!              to the value of highest model level
      !!
      !! ** Method  : Remove the observations that are outside the model space
      !!              and time domain or located within model land cells.
      !!
      !!              NB. T and S profile observations lying between the ocean
      !!              surface and the depth of the first model T point are 
      !!              assigned a depth equal to that of the first model T pt.
      !!
      !! ** Action  : 
      !!   
      !! History :
      !!        !  2007-01  (K. Mogensen) Rewrite of parts of obs_scr
      !!        !  2007-06  (K. Mogensen et al) Reject obs. near land.
      !!----------------------------------------------------------------------
      !! * Modules used
      USE dom_oce, ONLY : &       ! Geographical information 
         & gdepw_1d,      &
         & gdepw_0,       &                       
-        & gdepw,         &
+        & gdepw, r3t,    &
         & gdept,         &
         & ln_zco,        &
         & ln_zps             

      !! * Arguments
      INTEGER, INTENT(IN) :: kprofno      ! Number of profiles
      INTEGER, INTENT(IN) :: kobsno       ! Total number of observations
      INTEGER, INTENT(IN) :: kpi          ! Number of grid points in (i,j,k)
      INTEGER, INTENT(IN) :: kpj
      INTEGER, INTENT(IN) :: kpk    
      INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
         & kpstart, &         ! Start of individual profiles
         & kpend              ! End of individual profiles
      INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
         & kobsi, &           ! Observation (i,j) coordinates
         & kobsj
      INTEGER, DIMENSION(kobsno), INTENT(IN) :: &
         & kobsk              ! Observation k coordinate
      REAL(KIND=wp), DIMENSION(kprofno), INTENT(IN) :: &
         & pobslam, &         ! Observation (lon,lat) coordinates
         & pobsphi
      REAL(KIND=wp), DIMENSION(kobsno), INTENT(INOUT) :: &
         & pobsdep            ! Observation depths  
      REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
         & plam, pphi         ! Model (lon,lat) coordinates
      REAL(KIND=wp), DIMENSION(kpk), INTENT(IN) :: &
         & pdep               ! Model depth coordinates
      REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) :: &
         & pmask              ! Land mask array
      INTEGER, DIMENSION(kprofno), INTENT(INOUT) :: &
         & kpobsqc             ! Profile quality control
      INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
         & kobsqc             ! Observation quality control
      INTEGER, INTENT(INOUT) :: kosdobs     ! Observations outside space domain
      INTEGER, INTENT(INOUT) :: klanobs     ! Observations within a model land cell
      INTEGER, INTENT(INOUT) :: knlaobs     ! Observations near land
      INTEGER, INTENT(INOUT) :: kbdyobs     ! Observations near boundary
      LOGICAL, INTENT(IN) :: ld_nea         ! Flag observations near land
      LOGICAL, INTENT(IN) :: ld_bound_reject  ! Flag observations near open boundary
      INTEGER, INTENT(IN) :: kqc_cutoff     ! Cutoff QC value
      INTEGER, INTENT(IN) :: Kmm            ! time-level index

      !! * Local declarations
      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
         & zgmsk              ! Grid mask
      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
         & zbmsk              ! Boundary mask
      REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask
      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
         & zgdepw
      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
         & zglam, &           ! Model longitude at grid points
         & zgphi              ! Model latitude at grid points
+     REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdepw
      INTEGER, DIMENSION(2,2,kprofno) :: &
         & igrdi, &           ! Grid i,j
         & igrdj
      LOGICAL :: lgridobs           ! Is observation on a model grid point.
      LOGICAL :: ll_next_to_land    ! Is a profile next to land 
      INTEGER :: iig, ijg           ! i,j of observation on model grid point.
      INTEGER :: jobs, jobsp, jk, ji, jj
      !!----------------------------------------------------------------------

      ! Get grid point indices
      
      DO jobs = 1, kprofno

         ! For invalid points use 2,2

         IF ( kpobsqc(jobs) >= kqc_cutoff ) THEN
            
            igrdi(1,1,jobs) = 1
            igrdj(1,1,jobs) = 1
            igrdi(1,2,jobs) = 1
            igrdj(1,2,jobs) = 2
            igrdi(2,1,jobs) = 2
            igrdj(2,1,jobs) = 1
            igrdi(2,2,jobs) = 2
            igrdj(2,2,jobs) = 2
            
         ELSE
            
            igrdi(1,1,jobs) = kobsi(jobs)-1
            igrdj(1,1,jobs) = kobsj(jobs)-1
            igrdi(1,2,jobs) = kobsi(jobs)-1
            igrdj(1,2,jobs) = kobsj(jobs)
            igrdi(2,1,jobs) = kobsi(jobs)
            igrdj(2,1,jobs) = kobsj(jobs)-1
            igrdi(2,2,jobs) = kobsi(jobs)
            igrdj(2,2,jobs) = kobsj(jobs)
            
         ENDIF
         
      END DO

      IF (ln_bdy) THEN
        ! Create a mask grid points in boundary rim
        IF (ld_bound_reject) THEN           
           zbdymask(:,:) = 1.0_wp
           DO ji = 1, nb_bdy
              DO jj = 1, idx_bdy(ji)%nblen(1)
                 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp
              ENDDO
           ENDDO
        ENDIF
   
        CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk )
      ENDIF
      
      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk )
      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam )
      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi )
-     CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw(:,:,:,Kmm), &
-       &                     zgdepw )
+     DO jk = 1, jpk
+        zdepw(:,:,jk) = gdepw(:,:,jk,Kmm)
+     END DO
+     CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, zdepw(:,:,:), zgdepw )

      DO jobs = 1, kprofno

Commit History (1)

ChangesetAuthorTimeChangeLog
14982hadcv2021-06-11T16:52:03+02:00

#2665: Various fixes for code enabled with key_qco/key_linssh

Change History (5)

comment:1 Changed 7 weeks ago by techene

Yes, you got it exactly. Your fix is correct.
With qco, most of vertical grid variables no longer exist in array form.
Each of them is replaced by an expression defined in domzgr_substitute.h90.

comment:2 Changed 7 weeks ago by hadcv

Great, thanks Sibylle!

Last edited 7 weeks ago by hadcv (previous) (diff)

comment:3 Changed 10 days ago by hadcv

Additional fixes are needed in tradmp.F90 for the diagnostics added by [14718]:

      INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
      REAL(wp), DIMENSION(A2D(nn_hls),jpk,jpts)     ::  zts_dta
+      REAL(wp), DIMENSION(:,:,:)  , ALLOCATABLE ::  zwrk
      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  ztrdts
      IF( l_trdtra .OR. iom_use('hflx_dmp_cea') .OR. iom_use('sflx_dmp_cea') ) THEN   !* Save ta and sa trends
-         ALLOCATE( ztrdts(jpi,jpj,jpk,jpts) )
+         ALLOCATE( ztrdts(A2D(nn_hls),jpk,jpts) )
-         ztrdts(:,:,:,:) = pts(:,:,:,:,Krhs)
+         DO jn = 1, jpts
+            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
+               ztrdts(ji,jj,jk,jn) = pts(ji,jj,jk,jn,Krhs)
+            END_3D
+         END DO
      ENDIF
      ! outputs (clem trunk)
-      IF( iom_use('hflx_dmp_cea') )       &
-         &   CALL iom_put('hflx_dmp_cea', &
-         &   SUM( ( pts(:,:,:,jp_tem,Krhs) - ztrdts(:,:,:,jp_tem) ) * e3t(:,:,:,Kmm), dim=3 ) * rcp * rho0 ) ! W/m2
-      IF( iom_use('sflx_dmp_cea') )       &
-         &   CALL iom_put('sflx_dmp_cea', &
-         &   SUM( ( pts(:,:,:,jp_sal,Krhs) - ztrdts(:,:,:,jp_sal) ) * e3t(:,:,:,Kmm), dim=3 ) * rho0 )       ! g/m2/s
+      IF( iom_use('hflx_dmp_cea') .OR. iom_use('sflx_dmp_cea') ) THEN
+         ALLOCATE( zwrk(A2D(nn_hls),jpk) )          ! Needed to handle expressions containing e3t when using key_qco or key_linssh
+         zwrk(:,:,:) = 0._wp

+         IF( iom_use('hflx_dmp_cea') ) THEN
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               zwrk(ji,jj,jk) = ( pts(ji,jj,jk,jp_tem,Krhs) - ztrdts(ji,jj,jk,jp_tem) ) * e3t(ji,jj,jk,Kmm)
+            END_3D
+            CALL iom_put('hflx_dmp_cea', SUM( zwrk(:,:,:), dim=3 ) * rcp * rho0 ) ! W/m2
+         ENDIF
+         IF( iom_use('sflx_dmp_cea') ) THEN
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               zwrk(ji,jj,jk) = ( pts(ji,jj,jk,jp_sal,Krhs) - ztrdts(ji,jj,jk,jp_sal) ) * e3t(ji,jj,jk,Kmm)
+            END_3D
+            CALL iom_put('sflx_dmp_cea', SUM( zwrk(:,:,:), dim=3 ) * rho0 )       ! g/m2/s
+         ENDIF
+
+         DEALLOCATE( zwrk )
+      ENDIF

comment:4 Changed 10 days ago by hadcv

In 14982:

#2665: Various fixes for code enabled with key_qco/key_linssh

comment:5 Changed 10 days ago by hadcv

  • Resolution set to fixed
  • Status changed from new to closed
Note: See TracTickets for help on using tickets.