Opened 4 years ago
Closed 3 years 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 |
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)
| Changeset | Author | Time | ChangeLog |
|---|---|---|---|
| 14982 | hadcv | 2021-06-11T16:52:03+02:00 | #2665: Various fixes for code enabled with key_qco/key_linssh |
Change History (5)
comment:1 Changed 4 years ago by techene
comment:2 Changed 4 years ago by hadcv
Great, thanks Sibylle!
comment:3 Changed 3 years 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 3 years ago by hadcv
In 14982:
comment:5 Changed 3 years ago by hadcv
- Resolution set to fixed
- Status changed from new to closed
Note: See
TracTickets for help on using
tickets.
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.