Changeset 11361 for NEMO/branches/UKMO
- Timestamp:
- 2019-07-29T11:26:23+02:00 (5 years ago)
- Location:
- NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 2 added
- 13 edited
-
ASM/asminc.F90 (modified) (1 diff)
-
OBS/diaobs.F90 (modified) (26 diffs)
-
OBS/obs_averg_h2d.F90 (added)
-
OBS/obs_mpp.F90 (modified) (5 diffs)
-
OBS/obs_oper.F90 (modified) (24 diffs)
-
OBS/obs_prep.F90 (modified) (54 diffs)
-
OBS/obs_profiles_def.F90 (modified) (1 diff)
-
OBS/obs_read_altbias.F90 (modified) (1 diff)
-
OBS/obs_read_prof.F90 (modified) (23 diffs)
-
OBS/obs_read_surf.F90 (modified) (13 diffs)
-
OBS/obs_readmdt.F90 (modified) (3 diffs)
-
OBS/obs_surf_def.F90 (modified) (14 diffs)
-
OBS/obs_write.F90 (modified) (16 diffs)
-
OBS/obsaverg_h2d.h90 (added)
-
OBS/obsinter_h2d.h90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
r11350 r11361 71 71 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkginc, v_bkginc !: Increment to the u- & v-components 72 72 REAL(wp), PUBLIC, DIMENSION(:) , ALLOCATABLE :: wgtiau !: IAU weights for each time step 73 #if defined key_asminc 74 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ssh_iau !: IAU-weighted sea surface height increment 75 #endif 73 REAL(wp), PUBLIC, DIMENSION(:,:) , ALLOCATABLE :: ssh_iau !: IAU-weighted sea surface height increment 76 74 ! !!! time steps relative to the cycle interval [0,nitend-nit000-1] 77 75 INTEGER , PUBLIC :: nitbkg !: Time step of the background state used in the Jb term -
NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
r11350 r11361 46 46 47 47 !! * Module variables 48 LOGICAL, PUBLIC :: ln_diaobs !: Logical switch for the obs operator 49 LOGICAL :: ln_sstnight !: Logical switch for night mean SST obs 50 51 INTEGER :: nn_1dint !: Vertical interpolation method 52 INTEGER :: nn_2dint !: Horizontal interpolation method 48 LOGICAL, PUBLIC :: & 49 & ln_diaobs !: Logical switch for the obs operator 50 LOGICAL :: ln_sstnight !: Logical switch for night mean SST obs 51 LOGICAL :: ln_default_fp_indegs !: T=> Default obs footprint size specified in degrees, F=> in metres 52 LOGICAL :: ln_sla_fp_indegs !: T=> SLA obs footprint size specified in degrees, F=> in metres 53 LOGICAL :: ln_sst_fp_indegs !: T=> SST obs footprint size specified in degrees, F=> in metres 54 LOGICAL :: ln_sss_fp_indegs !: T=> SSS obs footprint size specified in degrees, F=> in metres 55 LOGICAL :: ln_sic_fp_indegs !: T=> sea-ice obs footprint size specified in degrees, F=> in metres 56 57 REAL(wp) :: rn_default_avglamscl !: Default E/W diameter of observation footprint 58 REAL(wp) :: rn_default_avgphiscl !: Default N/S diameter of observation footprint 59 REAL(wp) :: rn_sla_avglamscl !: E/W diameter of SLA observation footprint 60 REAL(wp) :: rn_sla_avgphiscl !: N/S diameter of SLA observation footprint 61 REAL(wp) :: rn_sst_avglamscl !: E/W diameter of SST observation footprint 62 REAL(wp) :: rn_sst_avgphiscl !: N/S diameter of SST observation footprint 63 REAL(wp) :: rn_sss_avglamscl !: E/W diameter of SSS observation footprint 64 REAL(wp) :: rn_sss_avgphiscl !: N/S diameter of SSS observation footprint 65 REAL(wp) :: rn_sic_avglamscl !: E/W diameter of sea-ice observation footprint 66 REAL(wp) :: rn_sic_avgphiscl !: N/S diameter of sea-ice observation footprint 67 68 INTEGER :: nn_1dint !: Vertical interpolation method 69 INTEGER :: nn_2dint_default !: Default horizontal interpolation method 70 INTEGER :: nn_2dint_sla !: SLA horizontal interpolation method (-1 = default) 71 INTEGER :: nn_2dint_sst !: SST horizontal interpolation method (-1 = default) 72 INTEGER :: nn_2dint_sss !: SSS horizontal interpolation method (-1 = default) 73 INTEGER :: nn_2dint_sic !: Seaice horizontal interpolation method (-1 = default) 74 53 75 INTEGER, DIMENSION(imaxavtypes) :: & 54 76 & nn_profdavtypes !: Profile data types representing a daily average … … 62 84 & nextrsurf !: Number of surface extra variables 63 85 INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: sstbias_type !SST bias type 86 INTEGER, DIMENSION(:), ALLOCATABLE :: & 87 & n2dintsurf !: Interpolation option for surface variables 88 REAL(wp), DIMENSION(:), ALLOCATABLE :: & 89 & ravglamscl, & !: E/W diameter of averaging footprint for surface variables 90 & ravgphiscl !: N/S diameter of averaging footprint for surface variables 91 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 92 & lfpindegs, & !: T=> surface obs footprint size specified in degrees, F=> in metres 93 & llnightav !: Logical for calculating night-time averages 64 94 TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: & 65 95 & surfdata, & !: Initial surface data … … 69 99 & profdataqc !: Profile data after quality control 70 100 71 CHARACTER(len= 6), PUBLIC, DIMENSION(:), ALLOCATABLE :: &101 CHARACTER(len=8), PUBLIC, DIMENSION(:), ALLOCATABLE :: & 72 102 & cobstypesprof, & !: Profile obs types 73 103 & cobstypessurf !: Surface obs types … … 96 126 !! ! 06-10 (A. Weaver) Cleaning and add controls 97 127 !! ! 07-03 (K. Mogensen) General handling of profiles 98 !! ! 14-08 (J.While) Incorporated SST bias correction 128 !! ! 14-08 (J.While) Incorporated SST bias correction 99 129 !! ! 15-02 (M. Martin) Simplification of namelist and code 100 130 !!---------------------------------------------------------------------- … … 112 142 INTEGER :: jvar ! Counter for variables 113 143 INTEGER :: jfile ! Counter for files 144 INTEGER :: jnumsstbias ! Number of SST bias files to read and apply 145 INTEGER :: n2dint_type ! Local version of nn_2dint* 114 146 115 147 CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: & 116 & cn_profbfiles, & ! T/S profile input filenames 117 & cn_sstfbfiles, & ! Sea surface temperature input filenames 118 & cn_slafbfiles, & ! Sea level anomaly input filenames 119 & cn_sicfbfiles, & ! Seaice concentration input filenames 120 & cn_velfbfiles, & ! Velocity profile input filenames 121 & cn_sstbias_files ! SST bias input filenames 148 & cn_profbfiles, & ! T/S profile input filenames 149 & cn_sstfbfiles, & ! Sea surface temperature input filenames 150 & cn_slafbfiles, & ! Sea level anomaly input filenames 151 & cn_sicfbfiles, & ! Seaice concentration input filenames 152 & cn_velfbfiles, & ! Velocity profile input filenames 153 & cn_sssfbfiles, & ! Sea surface salinity input filenames 154 & cn_slchltotfbfiles, & ! Surface total log10(chlorophyll) input filenames 155 & cn_slchldiafbfiles, & ! Surface diatom log10(chlorophyll) input filenames 156 & cn_slchlnonfbfiles, & ! Surface non-diatom log10(chlorophyll) input filenames 157 & cn_slchldinfbfiles, & ! Surface dinoflagellate log10(chlorophyll) input filenames 158 & cn_slchlmicfbfiles, & ! Surface microphytoplankton log10(chlorophyll) input filenames 159 & cn_slchlnanfbfiles, & ! Surface nanophytoplankton log10(chlorophyll) input filenames 160 & cn_slchlpicfbfiles, & ! Surface picophytoplankton log10(chlorophyll) input filenames 161 & cn_schltotfbfiles, & ! Surface total chlorophyll input filenames 162 & cn_slphytotfbfiles, & ! Surface total log10(phytoplankton carbon) input filenames 163 & cn_slphydiafbfiles, & ! Surface diatom log10(phytoplankton carbon) input filenames 164 & cn_slphynonfbfiles, & ! Surface non-diatom log10(phytoplankton carbon) input filenames 165 & cn_sspmfbfiles, & ! Surface suspended particulate matter input filenames 166 & cn_sfco2fbfiles, & ! Surface fugacity of carbon dioxide input filenames 167 & cn_spco2fbfiles, & ! Surface partial pressure of carbon dioxide input filenames 168 & cn_plchltotfbfiles, & ! Profile total log10(chlorophyll) input filenames 169 & cn_pchltotfbfiles, & ! Profile total chlorophyll input filenames 170 & cn_pno3fbfiles, & ! Profile nitrate input filenames 171 & cn_psi4fbfiles, & ! Profile silicate input filenames 172 & cn_ppo4fbfiles, & ! Profile phosphate input filenames 173 & cn_pdicfbfiles, & ! Profile dissolved inorganic carbon input filenames 174 & cn_palkfbfiles, & ! Profile alkalinity input filenames 175 & cn_pphfbfiles, & ! Profile pH input filenames 176 & cn_po2fbfiles, & ! Profile dissolved oxygen input filenames 177 & cn_sstbiasfiles ! SST bias input filenames 178 122 179 CHARACTER(LEN=128) :: & 123 180 & cn_altbiasfile ! Altimeter bias input filename 124 CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: &125 & clproffiles, & ! Profile filenames126 & clsurffiles ! Surface filenames127 181 128 182 LOGICAL :: ln_t3d ! Logical switch for temperature profiles … … 131 185 LOGICAL :: ln_sst ! Logical switch for sea surface temperature 132 186 LOGICAL :: ln_sic ! Logical switch for sea ice concentration 187 LOGICAL :: ln_sss ! Logical switch for sea surface salinity obs 133 188 LOGICAL :: ln_vel3d ! Logical switch for velocity (u,v) obs 189 LOGICAL :: ln_slchltot ! Logical switch for surface total log10(chlorophyll) obs 190 LOGICAL :: ln_slchldia ! Logical switch for surface diatom log10(chlorophyll) obs 191 LOGICAL :: ln_slchlnon ! Logical switch for surface non-diatom log10(chlorophyll) obs 192 LOGICAL :: ln_slchldin ! Logical switch for surface dinoflagellate log10(chlorophyll) obs 193 LOGICAL :: ln_slchlmic ! Logical switch for surface microphytoplankton log10(chlorophyll) obs 194 LOGICAL :: ln_slchlnan ! Logical switch for surface nanophytoplankton log10(chlorophyll) obs 195 LOGICAL :: ln_slchlpic ! Logical switch for surface picophytoplankton log10(chlorophyll) obs 196 LOGICAL :: ln_schltot ! Logical switch for surface total chlorophyll obs 197 LOGICAL :: ln_slphytot ! Logical switch for surface total log10(phytoplankton carbon) obs 198 LOGICAL :: ln_slphydia ! Logical switch for surface diatom log10(phytoplankton carbon) obs 199 LOGICAL :: ln_slphynon ! Logical switch for surface non-diatom log10(phytoplankton carbon) obs 200 LOGICAL :: ln_sspm ! Logical switch for surface suspended particulate matter obs 201 LOGICAL :: ln_sfco2 ! Logical switch for surface fugacity of carbon dioxide obs 202 LOGICAL :: ln_spco2 ! Logical switch for surface partial pressure of carbon dioxide obs 203 LOGICAL :: ln_plchltot ! Logical switch for profile total log10(chlorophyll) obs 204 LOGICAL :: ln_pchltot ! Logical switch for profile total chlorophyll obs 205 LOGICAL :: ln_pno3 ! Logical switch for profile nitrate obs 206 LOGICAL :: ln_psi4 ! Logical switch for profile silicate obs 207 LOGICAL :: ln_ppo4 ! Logical switch for profile phosphate obs 208 LOGICAL :: ln_pdic ! Logical switch for profile dissolved inorganic carbon obs 209 LOGICAL :: ln_palk ! Logical switch for profile alkalinity obs 210 LOGICAL :: ln_pph ! Logical switch for profile pH obs 211 LOGICAL :: ln_po2 ! Logical switch for profile dissolved oxygen obs 134 212 LOGICAL :: ln_nea ! Logical switch to remove obs near land 135 213 LOGICAL :: ln_altbias ! Logical switch for altimeter bias 136 LOGICAL :: ln_sstbias !: Logical switch for bias corection of SST214 LOGICAL :: ln_sstbias ! Logical switch for bias correction of SST 137 215 LOGICAL :: ln_ignmis ! Logical switch for ignoring missing files 138 216 LOGICAL :: ln_s_at_t ! Logical switch to compute model S at T obs 139 LOGICAL :: llvar1 ! Logical for profile variable 1 140 LOGICAL :: llvar2 ! Logical for profile variable 1 141 LOGICAL :: llnightav ! Logical for calculating night-time averages 217 LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary 142 218 LOGICAL, DIMENSION(jpmaxnfiles) :: lmask ! Used for finding number of sstbias files 143 219 144 220 REAL(dp) :: rn_dobsini ! Obs window start date YYYYMMDD.HHMMSS 145 221 REAL(dp) :: rn_dobsend ! Obs window end date YYYYMMDD.HHMMSS 146 REAL(wp), POINTER, DIMENSION(:,:) :: & 147 & zglam1, & ! Model longitudes for profile variable 1 148 & zglam2 ! Model longitudes for profile variable 2 149 REAL(wp), POINTER, DIMENSION(:,:) :: & 150 & zgphi1, & ! Model latitudes for profile variable 1 151 & zgphi2 ! Model latitudes for profile variable 2 222 223 REAL(wp) :: ztype_avglamscl ! Local version of rn_*_avglamscl 224 REAL(wp) :: ztype_avgphiscl ! Local version of rn_*_avgphiscl 225 226 CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: & 227 & clproffiles, & ! Profile filenames 228 & clsurffiles ! Surface filenames 229 230 LOGICAL, DIMENSION(:), ALLOCATABLE :: llvar ! Logical for profile variable read 231 LOGICAL :: ltype_fp_indegs ! Local version of ln_*_fp_indegs 232 LOGICAL :: ltype_night ! Local version of ln_sstnight (false for other variables) 233 152 234 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 153 & zmask1, & ! Model land/sea mask associated with variable 1 154 & zmask2 ! Model land/sea mask associated with variable 2 235 & zglam ! Model longitudes for profile variables 236 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 237 & zgphi ! Model latitudes for profile variables 238 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 239 & zmask ! Model land/sea mask associated with variables 240 155 241 156 242 NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla, & 157 & ln_sst, ln_sic, ln_vel3d, & 158 & ln_altbias, ln_nea, ln_grid_global, & 159 & ln_grid_search_lookup, & 160 & ln_ignmis, ln_s_at_t, ln_sstnight, & 243 & ln_sst, ln_sic, ln_sss, ln_vel3d, & 244 & ln_slchltot, ln_slchldia, ln_slchlnon, & 245 & ln_slchldin, ln_slchlmic, ln_slchlnan, & 246 & ln_slchlpic, ln_schltot, & 247 & ln_slphytot, ln_slphydia, ln_slphynon, & 248 & ln_sspm, ln_sfco2, ln_spco2, & 249 & ln_plchltot, ln_pchltot, ln_pno3, & 250 & ln_psi4, ln_ppo4, ln_pdic, & 251 & ln_palk, ln_pph, ln_po2, & 252 & ln_altbias, ln_sstbias, ln_nea, & 253 & ln_grid_global, ln_grid_search_lookup, & 254 & ln_ignmis, ln_s_at_t, ln_bound_reject, & 255 & ln_sstnight, ln_default_fp_indegs, & 256 & ln_sla_fp_indegs, ln_sst_fp_indegs, & 257 & ln_sss_fp_indegs, ln_sic_fp_indegs, & 161 258 & cn_profbfiles, cn_slafbfiles, & 162 259 & cn_sstfbfiles, cn_sicfbfiles, & 163 & cn_velfbfiles, cn_altbiasfile, & 260 & cn_velfbfiles, cn_sssfbfiles, & 261 & cn_slchltotfbfiles, cn_slchldiafbfiles, & 262 & cn_slchlnonfbfiles, cn_slchldinfbfiles, & 263 & cn_slchlmicfbfiles, cn_slchlnanfbfiles, & 264 & cn_slchlpicfbfiles, cn_schltotfbfiles, & 265 & cn_slphytotfbfiles, cn_slphydiafbfiles, & 266 & cn_slphynonfbfiles, cn_sspmfbfiles, & 267 & cn_sfco2fbfiles, cn_spco2fbfiles, & 268 & cn_plchltotfbfiles, cn_pchltotfbfiles, & 269 & cn_pno3fbfiles, cn_psi4fbfiles, cn_ppo4fbfiles, & 270 & cn_pdicfbfiles, cn_palkfbfiles, cn_pphfbfiles, & 271 & cn_po2fbfiles, & 272 & cn_sstbiasfiles, cn_altbiasfile, & 164 273 & cn_gridsearchfile, rn_gridsearchres, & 165 & rn_dobsini, rn_dobsend, nn_1dint, nn_2dint, & 274 & rn_dobsini, rn_dobsend, & 275 & rn_default_avglamscl, rn_default_avgphiscl, & 276 & rn_sla_avglamscl, rn_sla_avgphiscl, & 277 & rn_sst_avglamscl, rn_sst_avgphiscl, & 278 & rn_sss_avglamscl, rn_sss_avgphiscl, & 279 & rn_sic_avglamscl, rn_sic_avgphiscl, & 280 & nn_1dint, nn_2dint_default, & 281 & nn_2dint_sla, nn_2dint_sst, & 282 & nn_2dint_sss, nn_2dint_sic, & 166 283 & nn_msshc, rn_mdtcorr, rn_mdtcutoff, & 167 & nn_profdavtypes , ln_sstbias, cn_sstbias_files284 & nn_profdavtypes 168 285 169 286 INTEGER :: jnumsstbias 170 CALL wrk_alloc( jpi, jpj, zglam1 ) 171 CALL wrk_alloc( jpi, jpj, zglam2 ) 172 CALL wrk_alloc( jpi, jpj, zgphi1 ) 173 CALL wrk_alloc( jpi, jpj, zgphi2 ) 174 CALL wrk_alloc( jpi, jpj, jpk, zmask1 ) 175 CALL wrk_alloc( jpi, jpj, jpk, zmask2 ) 287 ALLOCATE(sstbias_type(jpmaxnfiles)) 176 288 177 289 !----------------------------------------------------------------------- 178 290 ! Read namelist parameters 179 291 !----------------------------------------------------------------------- 180 181 !Initalise all values in namelist arrays 182 ALLOCATE(sstbias_type(jpmaxnfiles)) 292 183 293 ! Some namelist arrays need initialising 184 cn_profbfiles(:) = '' 185 cn_slafbfiles(:) = '' 186 cn_sstfbfiles(:) = '' 187 cn_sicfbfiles(:) = '' 188 cn_velfbfiles(:) = '' 189 cn_sstbias_files(:) = '' 190 nn_profdavtypes(:) = -1 294 cn_profbfiles(:) = '' 295 cn_slafbfiles(:) = '' 296 cn_sstfbfiles(:) = '' 297 cn_sicfbfiles(:) = '' 298 cn_velfbfiles(:) = '' 299 cn_sssfbfiles(:) = '' 300 cn_slchltotfbfiles(:) = '' 301 cn_slchldiafbfiles(:) = '' 302 cn_slchlnonfbfiles(:) = '' 303 cn_slchldinfbfiles(:) = '' 304 cn_slchlmicfbfiles(:) = '' 305 cn_slchlnanfbfiles(:) = '' 306 cn_slchlpicfbfiles(:) = '' 307 cn_schltotfbfiles(:) = '' 308 cn_slphytotfbfiles(:) = '' 309 cn_slphydiafbfiles(:) = '' 310 cn_slphynonfbfiles(:) = '' 311 cn_sspmfbfiles(:) = '' 312 cn_sfco2fbfiles(:) = '' 313 cn_spco2fbfiles(:) = '' 314 cn_plchltotfbfiles(:) = '' 315 cn_pchltotfbfiles(:) = '' 316 cn_pno3fbfiles(:) = '' 317 cn_psi4fbfiles(:) = '' 318 cn_ppo4fbfiles(:) = '' 319 cn_pdicfbfiles(:) = '' 320 cn_palkfbfiles(:) = '' 321 cn_pphfbfiles(:) = '' 322 cn_po2fbfiles(:) = '' 323 cn_sstbiasfiles(:) = '' 324 nn_profdavtypes(:) = -1 191 325 192 326 CALL ini_date( rn_dobsini ) … … 205 339 IF ( .NOT. ln_diaobs ) THEN 206 340 IF(lwp) WRITE(numout,cform_war) 207 IF(lwp) WRITE(numout,*)' ln_diaobs is set to false so not calling dia_obs'341 IF(lwp) WRITE(numout,*)' ln_diaobs is set to false, so not calling dia_obs' 208 342 RETURN 209 ENDIF210 211 !-----------------------------------------------------------------------212 ! Set up list of observation types to be used213 ! and the files associated with each type214 !-----------------------------------------------------------------------215 216 nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) )217 nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic /) )218 219 IF (ln_sstbias) THEN220 lmask(:) = .FALSE.221 WHERE (cn_sstbias_files(:) /= '') lmask(:) = .TRUE.222 jnumsstbias = COUNT(lmask)223 lmask(:) = .FALSE.224 ENDIF225 226 IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN227 IF(lwp) WRITE(numout,cform_war)228 IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', &229 & ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', &230 & ' are set to .FALSE. so turning off calls to dia_obs'231 nwarn = nwarn + 1232 ln_diaobs = .FALSE.233 RETURN234 ENDIF235 236 IF ( nproftypes > 0 ) THEN237 238 ALLOCATE( cobstypesprof(nproftypes) )239 ALLOCATE( ifilesprof(nproftypes) )240 ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) )241 242 jtype = 0243 IF (ln_t3d .OR. ln_s3d) THEN244 jtype = jtype + 1245 clproffiles(jtype,:) = cn_profbfiles(:)246 cobstypesprof(jtype) = 'prof '247 ifilesprof(jtype) = 0248 DO jfile = 1, jpmaxnfiles249 IF ( trim(clproffiles(jtype,jfile)) /= '' ) &250 ifilesprof(jtype) = ifilesprof(jtype) + 1251 END DO252 ENDIF253 IF (ln_vel3d) THEN254 jtype = jtype + 1255 clproffiles(jtype,:) = cn_velfbfiles(:)256 cobstypesprof(jtype) = 'vel '257 ifilesprof(jtype) = 0258 DO jfile = 1, jpmaxnfiles259 IF ( trim(clproffiles(jtype,jfile)) /= '' ) &260 ifilesprof(jtype) = ifilesprof(jtype) + 1261 END DO262 ENDIF263 264 ENDIF265 266 IF ( nsurftypes > 0 ) THEN267 268 ALLOCATE( cobstypessurf(nsurftypes) )269 ALLOCATE( ifilessurf(nsurftypes) )270 ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) )271 272 jtype = 0273 IF (ln_sla) THEN274 jtype = jtype + 1275 clsurffiles(jtype,:) = cn_slafbfiles(:)276 cobstypessurf(jtype) = 'sla '277 ifilessurf(jtype) = 0278 DO jfile = 1, jpmaxnfiles279 IF ( trim(clsurffiles(jtype,jfile)) /= '' ) &280 ifilessurf(jtype) = ifilessurf(jtype) + 1281 END DO282 ENDIF283 IF (ln_sst) THEN284 jtype = jtype + 1285 clsurffiles(jtype,:) = cn_sstfbfiles(:)286 cobstypessurf(jtype) = 'sst '287 ifilessurf(jtype) = 0288 DO jfile = 1, jpmaxnfiles289 IF ( trim(clsurffiles(jtype,jfile)) /= '' ) &290 ifilessurf(jtype) = ifilessurf(jtype) + 1291 END DO292 ENDIF293 #if defined key_lim2 || defined key_lim3294 IF (ln_sic) THEN295 jtype = jtype + 1296 clsurffiles(jtype,:) = cn_sicfbfiles(:)297 cobstypessurf(jtype) = 'sic '298 ifilessurf(jtype) = 0299 DO jfile = 1, jpmaxnfiles300 IF ( trim(clsurffiles(jtype,jfile)) /= '' ) &301 ifilessurf(jtype) = ifilessurf(jtype) + 1302 END DO303 ENDIF304 #endif305 306 343 ENDIF 307 344 … … 318 355 WRITE(numout,*) ' Logical switch for Sea Ice observations ln_sic = ', ln_sic 319 356 WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d 320 WRITE(numout,*) ' Global distribution of observations ln_grid_global = ',ln_grid_global 321 WRITE(numout,*) ' Logical switch for SST bias correction ln_sstbias = ', ln_sstbias 322 WRITE(numout,*) ' Logical switch for obs grid search lookup ln_grid_search_lookup = ',ln_grid_search_lookup 357 WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss 358 WRITE(numout,*) ' Logical switch for surface total logchl obs ln_slchltot = ', ln_slchltot 359 WRITE(numout,*) ' Logical switch for surface diatom logchl obs ln_slchldia = ', ln_slchldia 360 WRITE(numout,*) ' Logical switch for surface non-diatom logchl obs ln_slchlnon = ', ln_slchlnon 361 WRITE(numout,*) ' Logical switch for surface dino logchl obs ln_slchldin = ', ln_slchldin 362 WRITE(numout,*) ' Logical switch for surface micro logchl obs ln_slchlmic = ', ln_slchlmic 363 WRITE(numout,*) ' Logical switch for surface nano logchl obs ln_slchlnan = ', ln_slchlnan 364 WRITE(numout,*) ' Logical switch for surface pico logchl obs ln_slchlpic = ', ln_slchlpic 365 WRITE(numout,*) ' Logical switch for surface total chl obs ln_schltot = ', ln_schltot 366 WRITE(numout,*) ' Logical switch for surface total log(phyC) obs ln_slphytot = ', ln_slphytot 367 WRITE(numout,*) ' Logical switch for surface diatom log(phyC) obs ln_slphydia = ', ln_slphydia 368 WRITE(numout,*) ' Logical switch for surface non-diatom log(phyC) obs ln_slphynon = ', ln_slphynon 369 WRITE(numout,*) ' Logical switch for surface SPM observations ln_sspm = ', ln_sspm 370 WRITE(numout,*) ' Logical switch for surface fCO2 observations ln_sfco2 = ', ln_sfco2 371 WRITE(numout,*) ' Logical switch for surface pCO2 observations ln_spco2 = ', ln_spco2 372 WRITE(numout,*) ' Logical switch for profile total logchl obs ln_plchltot = ', ln_plchltot 373 WRITE(numout,*) ' Logical switch for profile total chl obs ln_pchltot = ', ln_pchltot 374 WRITE(numout,*) ' Logical switch for profile nitrate obs ln_pno3 = ', ln_pno3 375 WRITE(numout,*) ' Logical switch for profile silicate obs ln_psi4 = ', ln_psi4 376 WRITE(numout,*) ' Logical switch for profile phosphate obs ln_ppo4 = ', ln_ppo4 377 WRITE(numout,*) ' Logical switch for profile DIC obs ln_pdic = ', ln_pdic 378 WRITE(numout,*) ' Logical switch for profile alkalinity obs ln_palk = ', ln_palk 379 WRITE(numout,*) ' Logical switch for profile pH obs ln_pph = ', ln_pph 380 WRITE(numout,*) ' Logical switch for profile oxygen obs ln_po2 = ', ln_po2 381 WRITE(numout,*) ' Global distribution of observations ln_grid_global = ', ln_grid_global 382 WRITE(numout,*) ' Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup 323 383 IF (ln_grid_search_lookup) & 324 384 WRITE(numout,*) ' Grid search lookup file header cn_gridsearchfile = ', cn_gridsearchfile … … 326 386 WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS rn_dobsend = ', rn_dobsend 327 387 WRITE(numout,*) ' Type of vertical interpolation method nn_1dint = ', nn_1dint 328 WRITE(numout,*) ' Type of horizontal interpolation method nn_2dint = ', nn_2dint 388 WRITE(numout,*) ' Default horizontal interpolation method nn_2dint_default = ', nn_2dint_default 389 WRITE(numout,*) ' Type of horizontal interpolation method for SLA nn_2dint_sla = ', nn_2dint_sla 390 WRITE(numout,*) ' Type of horizontal interpolation method for SST nn_2dint_sst = ', nn_2dint_sst 391 WRITE(numout,*) ' Type of horizontal interpolation method for SSS nn_2dint_sss = ', nn_2dint_sss 392 WRITE(numout,*) ' Type of horizontal interpolation method for SIC nn_2dint_sic = ', nn_2dint_sic 393 WRITE(numout,*) ' Default E/W diameter of obs footprint rn_default_avglamscl = ', rn_default_avglamscl 394 WRITE(numout,*) ' Default N/S diameter of obs footprint rn_default_avgphiscl = ', rn_default_avgphiscl 395 WRITE(numout,*) ' Default obs footprint in deg [T] or m [F] ln_default_fp_indegs = ', ln_default_fp_indegs 396 WRITE(numout,*) ' SLA E/W diameter of obs footprint rn_sla_avglamscl = ', rn_sla_avglamscl 397 WRITE(numout,*) ' SLA N/S diameter of obs footprint rn_sla_avgphiscl = ', rn_sla_avgphiscl 398 WRITE(numout,*) ' SLA obs footprint in deg [T] or m [F] ln_sla_fp_indegs = ', ln_sla_fp_indegs 399 WRITE(numout,*) ' SST E/W diameter of obs footprint rn_sst_avglamscl = ', rn_sst_avglamscl 400 WRITE(numout,*) ' SST N/S diameter of obs footprint rn_sst_avgphiscl = ', rn_sst_avgphiscl 401 WRITE(numout,*) ' SST obs footprint in deg [T] or m [F] ln_sst_fp_indegs = ', ln_sst_fp_indegs 402 WRITE(numout,*) ' SIC E/W diameter of obs footprint rn_sic_avglamscl = ', rn_sic_avglamscl 403 WRITE(numout,*) ' SIC N/S diameter of obs footprint rn_sic_avgphiscl = ', rn_sic_avgphiscl 404 WRITE(numout,*) ' SIC obs footprint in deg [T] or m [F] ln_sic_fp_indegs = ', ln_sic_fp_indegs 329 405 WRITE(numout,*) ' Rejection of observations near land switch ln_nea = ', ln_nea 406 WRITE(numout,*) ' Rejection of obs near open bdys ln_bound_reject = ', ln_bound_reject 330 407 WRITE(numout,*) ' MSSH correction scheme nn_msshc = ', nn_msshc 331 408 WRITE(numout,*) ' MDT correction rn_mdtcorr = ', rn_mdtcorr 332 409 WRITE(numout,*) ' MDT cutoff for computed correction rn_mdtcutoff = ', rn_mdtcutoff 333 410 WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias 411 WRITE(numout,*) ' Logical switch for sst bias ln_sstbias = ', ln_sstbias 334 412 WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis 335 413 WRITE(numout,*) ' Daily average types nn_profdavtypes = ', nn_profdavtypes 336 414 WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight 337 WRITE(numout,*) ' Number of profile obs types: ',nproftypes 338 339 IF ( nproftypes > 0 ) THEN 340 DO jtype = 1, nproftypes 341 DO jfile = 1, ifilesprof(jtype) 342 WRITE(numout,'(1X,2A)') ' '//cobstypesprof(jtype)//' input observation file names = ', & 343 TRIM(clproffiles(jtype,jfile)) 344 END DO 345 END DO 346 ENDIF 347 348 WRITE(numout,*)' Number of surface obs types: ',nsurftypes 349 IF ( nsurftypes > 0 ) THEN 350 DO jtype = 1, nsurftypes 351 DO jfile = 1, ifilessurf(jtype) 352 WRITE(numout,'(1X,2A)') ' '//cobstypessurf(jtype)//' input observation file names = ', & 353 TRIM(clsurffiles(jtype,jfile)) 354 END DO 355 END DO 356 ENDIF 357 WRITE(numout,*) '~~~~~~~~~~~~' 358 359 ENDIF 415 ENDIF 416 !----------------------------------------------------------------------- 417 ! Set up list of observation types to be used 418 ! and the files associated with each type 419 !----------------------------------------------------------------------- 420 421 nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d, ln_plchltot, & 422 & ln_pchltot, ln_pno3, ln_psi4, ln_ppo4, & 423 & ln_pdic, ln_palk, ln_pph, ln_po2 /) ) 424 nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss, & 425 & ln_slchltot, ln_slchldia, ln_slchlnon, ln_slchldin, & 426 & ln_slchlmic, ln_slchlnan, ln_slchlpic, ln_schltot, & 427 & ln_slphytot, ln_slphydia, ln_slphynon, ln_sspm, & 428 & ln_sfco2, ln_spco2 /) ) 429 430 IF (ln_sstbias) THEN 431 lmask(:) = .FALSE. 432 WHERE (cn_sstbiasfiles(:) /= '') lmask(:) = .TRUE. 433 jnumsstbias = COUNT(lmask) 434 lmask(:) = .FALSE. 435 ENDIF 436 437 IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 438 IF(lwp) WRITE(numout,cform_war) 439 IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 440 & ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', & 441 & ' are set to .FALSE. so turning off calls to dia_obs' 442 nwarn = nwarn + 1 443 ln_diaobs = .FALSE. 444 RETURN 445 ENDIF 446 447 IF(lwp) WRITE(numout,*) ' Number of profile obs types: ',nproftypes 448 IF ( nproftypes > 0 ) THEN 449 450 ALLOCATE( cobstypesprof(nproftypes) ) 451 ALLOCATE( ifilesprof(nproftypes) ) 452 ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 453 454 jtype = 0 455 IF (ln_t3d .OR. ln_s3d) THEN 456 jtype = jtype + 1 457 cobstypesprof(jtype) = 'prof' 458 clproffiles(jtype,:) = cn_profbfiles 459 ENDIF 460 IF (ln_vel3d) THEN 461 jtype = jtype + 1 462 cobstypesprof(jtype) = 'vel' 463 clproffiles(jtype,:) = cn_velfbfiles 464 ENDIF 465 IF (ln_plchltot) THEN 466 jtype = jtype + 1 467 cobstypesprof(jtype) = 'plchltot' 468 clproffiles(jtype,:) = cn_plchltotfbfiles 469 ENDIF 470 IF (ln_pchltot) THEN 471 jtype = jtype + 1 472 cobstypesprof(jtype) = 'pchltot' 473 clproffiles(jtype,:) = cn_pchltotfbfiles 474 ENDIF 475 IF (ln_pno3) THEN 476 jtype = jtype + 1 477 cobstypesprof(jtype) = 'pno3' 478 clproffiles(jtype,:) = cn_pno3fbfiles 479 ENDIF 480 IF (ln_psi4) THEN 481 jtype = jtype + 1 482 cobstypesprof(jtype) = 'psi4' 483 clproffiles(jtype,:) = cn_psi4fbfiles 484 ENDIF 485 IF (ln_ppo4) THEN 486 jtype = jtype + 1 487 cobstypesprof(jtype) = 'ppo4' 488 clproffiles(jtype,:) = cn_ppo4fbfiles 489 ENDIF 490 IF (ln_pdic) THEN 491 jtype = jtype + 1 492 cobstypesprof(jtype) = 'pdic' 493 clproffiles(jtype,:) = cn_pdicfbfiles 494 ENDIF 495 IF (ln_palk) THEN 496 jtype = jtype + 1 497 cobstypesprof(jtype) = 'palk' 498 clproffiles(jtype,:) = cn_palkfbfiles 499 ENDIF 500 IF (ln_pph) THEN 501 jtype = jtype + 1 502 cobstypesprof(jtype) = 'pph' 503 clproffiles(jtype,:) = cn_pphfbfiles 504 ENDIF 505 IF (ln_po2) THEN 506 jtype = jtype + 1 507 cobstypesprof(jtype) = 'po2' 508 clproffiles(jtype,:) = cn_po2fbfiles 509 ENDIF 510 511 CALL obs_settypefiles( nproftypes, jpmaxnfiles, ifilesprof, cobstypesprof, clproffiles ) 512 513 ENDIF 514 515 IF(lwp) WRITE(numout,*)' Number of surface obs types: ',nsurftypes 516 IF ( nsurftypes > 0 ) THEN 517 518 ALLOCATE( cobstypessurf(nsurftypes) ) 519 ALLOCATE( ifilessurf(nsurftypes) ) 520 ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 521 ALLOCATE(n2dintsurf(nsurftypes)) 522 ALLOCATE(ravglamscl(nsurftypes)) 523 ALLOCATE(ravgphiscl(nsurftypes)) 524 ALLOCATE(lfpindegs(nsurftypes)) 525 ALLOCATE(llnightav(nsurftypes)) 526 527 jtype = 0 528 IF (ln_sla) THEN 529 jtype = jtype + 1 530 cobstypessurf(jtype) = 'sla' 531 clsurffiles(jtype,:) = cn_slafbfiles 532 ENDIF 533 IF (ln_sst) THEN 534 jtype = jtype + 1 535 cobstypessurf(jtype) = 'sst' 536 clsurffiles(jtype,:) = cn_sstfbfiles 537 ENDIF 538 #if defined key_lim2 || defined key_lim3 539 IF (ln_sic) THEN 540 jtype = jtype + 1 541 cobstypessurf(jtype) = 'sic' 542 clsurffiles(jtype,:) = cn_sicfbfiles 543 ENDIF 544 #endif 545 IF (ln_sss) THEN 546 jtype = jtype + 1 547 cobstypessurf(jtype) = 'sss' 548 clsurffiles(jtype,:) = cn_sssfbfiles 549 ENDIF 550 IF (ln_slchltot) THEN 551 jtype = jtype + 1 552 cobstypessurf(jtype) = 'slchltot' 553 clsurffiles(jtype,:) = cn_slchltotfbfiles 554 ENDIF 555 IF (ln_slchldia) THEN 556 jtype = jtype + 1 557 cobstypessurf(jtype) = 'slchldia' 558 clsurffiles(jtype,:) = cn_slchldiafbfiles 559 ENDIF 560 IF (ln_slchlnon) THEN 561 jtype = jtype + 1 562 cobstypessurf(jtype) = 'slchlnon' 563 clsurffiles(jtype,:) = cn_slchlnonfbfiles 564 ENDIF 565 IF (ln_slchldin) THEN 566 jtype = jtype + 1 567 cobstypessurf(jtype) = 'slchldin' 568 clsurffiles(jtype,:) = cn_slchldinfbfiles 569 ENDIF 570 IF (ln_slchlmic) THEN 571 jtype = jtype + 1 572 cobstypessurf(jtype) = 'slchlmic' 573 clsurffiles(jtype,:) = cn_slchlmicfbfiles 574 ENDIF 575 IF (ln_slchlnan) THEN 576 jtype = jtype + 1 577 cobstypessurf(jtype) = 'slchlnan' 578 clsurffiles(jtype,:) = cn_slchlnanfbfiles 579 ENDIF 580 IF (ln_slchlpic) THEN 581 jtype = jtype + 1 582 cobstypessurf(jtype) = 'slchlpic' 583 clsurffiles(jtype,:) = cn_slchlpicfbfiles 584 ENDIF 585 IF (ln_schltot) THEN 586 jtype = jtype + 1 587 cobstypessurf(jtype) = 'schltot' 588 clsurffiles(jtype,:) = cn_schltotfbfiles 589 ENDIF 590 IF (ln_slphytot) THEN 591 jtype = jtype + 1 592 cobstypessurf(jtype) = 'slphytot' 593 clsurffiles(jtype,:) = cn_slphytotfbfiles 594 ENDIF 595 IF (ln_slphydia) THEN 596 jtype = jtype + 1 597 cobstypessurf(jtype) = 'slphydia' 598 clsurffiles(jtype,:) = cn_slphydiafbfiles 599 ENDIF 600 IF (ln_slphynon) THEN 601 jtype = jtype + 1 602 cobstypessurf(jtype) = 'slphynon' 603 clsurffiles(jtype,:) = cn_slphynonfbfiles 604 ENDIF 605 IF (ln_sspm) THEN 606 jtype = jtype + 1 607 cobstypessurf(jtype) = 'sspm' 608 clsurffiles(jtype,:) = cn_sspmfbfiles 609 ENDIF 610 IF (ln_sfco2) THEN 611 jtype = jtype + 1 612 cobstypessurf(jtype) = 'sfco2' 613 clsurffiles(jtype,:) = cn_sfco2fbfiles 614 ENDIF 615 IF (ln_spco2) THEN 616 jtype = jtype + 1 617 cobstypessurf(jtype) = 'spco2' 618 clsurffiles(jtype,:) = cn_spco2fbfiles 619 ENDIF 620 621 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, ifilessurf, cobstypessurf, clsurffiles ) 622 623 DO jtype = 1, nsurftypes 624 625 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 626 IF ( nn_2dint_sla == -1 ) THEN 627 n2dint_type = nn_2dint_default 628 ELSE 629 n2dint_type = nn_2dint_sla 630 ENDIF 631 ztype_avglamscl = rn_sla_avglamscl 632 ztype_avgphiscl = rn_sla_avgphiscl 633 ltype_fp_indegs = ln_sla_fp_indegs 634 ltype_night = .FALSE. 635 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN 636 IF ( nn_2dint_sst == -1 ) THEN 637 n2dint_type = nn_2dint_default 638 ELSE 639 n2dint_type = nn_2dint_sst 640 ENDIF 641 ztype_avglamscl = rn_sst_avglamscl 642 ztype_avgphiscl = rn_sst_avgphiscl 643 ltype_fp_indegs = ln_sst_fp_indegs 644 ltype_night = ln_sstnight 645 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN 646 IF ( nn_2dint_sic == -1 ) THEN 647 n2dint_type = nn_2dint_default 648 ELSE 649 n2dint_type = nn_2dint_sic 650 ENDIF 651 ztype_avglamscl = rn_sic_avglamscl 652 ztype_avgphiscl = rn_sic_avgphiscl 653 ltype_fp_indegs = ln_sic_fp_indegs 654 ltype_night = .FALSE. 655 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN 656 IF ( nn_2dint_sss == -1 ) THEN 657 n2dint_type = nn_2dint_default 658 ELSE 659 n2dint_type = nn_2dint_sss 660 ENDIF 661 ztype_avglamscl = rn_sss_avglamscl 662 ztype_avgphiscl = rn_sss_avgphiscl 663 ltype_fp_indegs = ln_sss_fp_indegs 664 ltype_night = .FALSE. 665 ELSE 666 n2dint_type = nn_2dint_default 667 ztype_avglamscl = rn_default_avglamscl 668 ztype_avgphiscl = rn_default_avgphiscl 669 ltype_fp_indegs = ln_default_fp_indegs 670 ltype_night = .FALSE. 671 ENDIF 672 673 CALL obs_setinterpopts( nsurftypes, jtype, TRIM(cobstypessurf(jtype)), & 674 & nn_2dint_default, n2dint_type, & 675 & ztype_avglamscl, ztype_avgphiscl, & 676 & ltype_fp_indegs, ltype_night, & 677 & n2dintsurf, ravglamscl, ravgphiscl, & 678 & lfpindegs, llnightav ) 679 680 END DO 681 682 ENDIF 683 684 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 685 360 686 361 687 !----------------------------------------------------------------------- … … 377 703 ENDIF 378 704 379 IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 4) ) THEN380 CALL ctl_stop(' Choice of horizontal (2D) interpolation method', &705 IF ( ( nn_2dint_default < 0 ) .OR. ( nn_2dint_default > 6 ) ) THEN 706 CALL ctl_stop(' Choice of default horizontal (2D) interpolation method', & 381 707 & ' is not available') 382 708 ENDIF … … 402 728 DO jtype = 1, nproftypes 403 729 404 nvarsprof(jtype) = 2405 730 IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN 731 nvarsprof(jtype) = 2 406 732 nextrprof(jtype) = 1 407 llvar1 = ln_t3d 408 llvar2 = ln_s3d 409 zglam1 = glamt 410 zgphi1 = gphit 411 zmask1 = tmask 412 zglam2 = glamt 413 zgphi2 = gphit 414 zmask2 = tmask 415 ENDIF 416 IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN 733 ALLOCATE(llvar(nvarsprof(jtype))) 734 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zglam ) 735 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zgphi ) 736 CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) 737 llvar(1) = ln_t3d 738 llvar(2) = ln_s3d 739 zglam(:,:,1) = glamt(:,:) 740 zglam(:,:,2) = glamt(:,:) 741 zgphi(:,:,1) = gphit(:,:) 742 zgphi(:,:,2) = gphit(:,:) 743 zmask(:,:,:,1) = tmask(:,:,:) 744 zmask(:,:,:,2) = tmask(:,:,:) 745 ELSE IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN 746 nvarsprof(jtype) = 2 417 747 nextrprof(jtype) = 2 418 llvar1 = ln_vel3d 419 llvar2 = ln_vel3d 420 zglam1 = glamu 421 zgphi1 = gphiu 422 zmask1 = umask 423 zglam2 = glamv 424 zgphi2 = gphiv 425 zmask2 = vmask 748 ALLOCATE(llvar(nvarsprof(jtype))) 749 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zglam ) 750 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zgphi ) 751 CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) 752 llvar(1) = ln_vel3d 753 llvar(2) = ln_vel3d 754 zglam(:,:,1) = glamu(:,:) 755 zglam(:,:,2) = glamv(:,:) 756 zgphi(:,:,1) = gphiu(:,:) 757 zgphi(:,:,2) = gphiv(:,:) 758 zmask(:,:,:,1) = umask(:,:,:) 759 zmask(:,:,:,2) = vmask(:,:,:) 760 ELSE 761 nvarsprof(jtype) = 1 762 nextrprof(jtype) = 0 763 ALLOCATE(llvar(nvarsprof(jtype))) 764 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zglam ) 765 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zgphi ) 766 CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) 767 llvar(1) = .TRUE. 768 zglam(:,:,1) = glamt(:,:) 769 zgphi(:,:,1) = gphit(:,:) 770 zmask(:,:,:,1) = tmask(:,:,:) 426 771 ENDIF 427 772 … … 430 775 & clproffiles(jtype,1:ifilesprof(jtype)), & 431 776 & nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, & 432 & rn_dobsini, rn_dobsend, llvar 1, llvar2, &777 & rn_dobsini, rn_dobsend, llvar, & 433 778 & ln_ignmis, ln_s_at_t, .FALSE., & 434 779 & kdailyavtypes = nn_profdavtypes ) … … 439 784 440 785 CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), & 441 & llvar 1, llvar2, &786 & llvar, & 442 787 & jpi, jpj, jpk, & 443 & zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2, & 444 & ln_nea, kdailyavtypes = nn_profdavtypes ) 788 & zmask, zglam, zgphi, & 789 & ln_nea, ln_bound_reject, & 790 & kdailyavtypes = nn_profdavtypes ) 791 792 DEALLOCATE( llvar ) 793 CALL wrk_dealloc( jpi, jpj, nvarsprof(jtype), zglam ) 794 CALL wrk_dealloc( jpi, jpj, nvarsprof(jtype), zgphi ) 795 CALL wrk_dealloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) 445 796 446 797 END DO … … 461 812 nvarssurf(jtype) = 1 462 813 nextrsurf(jtype) = 0 463 llnightav = .FALSE.464 814 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2 465 IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav = ln_sstnight815 IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav(jtype) = ln_sstnight 466 816 467 817 !Read in surface obs types … … 469 819 & clsurffiles(jtype,1:ifilessurf(jtype)), & 470 820 & nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 471 & rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav ) 472 473 474 CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea ) 821 & rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) ) 822 823 CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject ) 475 824 476 825 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 477 CALL obs_rea_mdt( surfdataqc(jtype), nn_2dint ) 478 IF ( ln_altbias ) CALL obs_rea_altbias ( surfdataqc(jtype), nn_2dint, cn_altbiasfile ) 826 CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) ) 827 IF ( ln_altbias ) & 828 & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile ) 479 829 ENDIF 480 830 IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN 481 !Read in bias field and correct SST. 482 IF ( jnumsstbias == 0 ) CALL ctl_stop("ln_sstbias set,"// & 483 " but no bias"// & 484 " files to read in") 485 CALL obs_app_sstbias( surfdataqc(jtype), nn_2dint, & 486 jnumsstbias, cn_sstbias_files(1:jnumsstbias) ) 831 jnumsstbias = 0 832 DO jfile = 1, jpmaxnfiles 833 IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) & 834 & jnumsstbias = jnumsstbias + 1 835 END DO 836 IF ( jnumsstbias == 0 ) THEN 837 CALL ctl_stop("ln_sstbias set but no bias files to read in") 838 ENDIF 839 840 CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), & 841 & jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) ) 842 487 843 ENDIF 488 844 END DO … … 491 847 492 848 ENDIF 493 494 CALL wrk_dealloc( jpi, jpj, zglam1 )495 CALL wrk_dealloc( jpi, jpj, zglam2 )496 CALL wrk_dealloc( jpi, jpj, zgphi1 )497 CALL wrk_dealloc( jpi, jpj, zgphi2 )498 CALL wrk_dealloc( jpi, jpj, jpk, zmask1 )499 CALL wrk_dealloc( jpi, jpj, jpk, zmask2 )500 849 501 850 END SUBROUTINE dia_obs_init … … 526 875 !!---------------------------------------------------------------------- 527 876 !! * Modules used 528 USE dom_oce, ONLY : & ! Ocean space and time domain variables529 & gdept_n, &530 & gdept_1d531 USE phycst, ONLY : & ! Physical constants532 & rday533 USE oce, ONLY : & ! Ocean dynamics and tracers variables534 & tsn, &535 & un, vn, &536 & sshn537 877 USE phycst, ONLY : & ! Physical constants 878 #if defined key_fabm 879 & rt0, & 880 #endif 538 881 & rday 882 USE oce, ONLY : & ! Ocean dynamics and tracers variables 883 & tsn, & 884 & un, & 885 & vn, & 886 & sshn 539 887 #if defined key_lim3 540 888 USE ice, ONLY : & ! LIM3 Ice model variables … … 545 893 & frld 546 894 #endif 895 #if defined key_cice 896 USE sbc_oce, ONLY : fr_i ! ice fraction 897 #endif 898 #if defined key_top 899 USE trc, ONLY : & ! Biogeochemical state variables 900 & trn 901 #endif 902 #if defined key_hadocc 903 USE par_hadocc ! HadOCC parameters 904 USE trc, ONLY : & 905 & HADOCC_CHL, & 906 & HADOCC_FCO2, & 907 & HADOCC_PCO2, & 908 & HADOCC_FILL_FLT 909 USE had_bgc_const, ONLY: c2n_p 910 #elif defined key_medusa 911 USE par_medusa ! MEDUSA parameters 912 USE sms_medusa, ONLY: & 913 & xthetapn, & 914 & xthetapd 915 #if defined key_roam 916 USE sms_medusa, ONLY: & 917 & f2_pco2w, & 918 & f2_fco2w, & 919 & f3_pH 920 #endif 921 #elif defined key_fabm 922 USE par_fabm ! FABM parameters 923 USE fabm, ONLY: & 924 & fabm_get_interior_diagnostic_data 925 #endif 926 #if defined key_spm 927 USE par_spm, ONLY: & ! Sediment parameters 928 & jp_spm 929 #endif 930 547 931 IMPLICIT NONE 548 932 … … 553 937 INTEGER :: jtype ! Data loop variable 554 938 INTEGER :: jvar ! Variable number 555 INTEGER :: ji, jj ! Loop counters 939 INTEGER :: ji, jj, jk ! Loop counters 940 REAL(wp) :: tiny ! small number 941 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 942 & zprofvar ! Model values for variables in a prof ob 943 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 944 & zprofmask ! Mask associated with zprofvar 945 REAL(wp), POINTER, DIMENSION(:,:) :: & 946 & zsurfvar, & ! Model values equivalent to surface ob. 947 & zsurfmask ! Mask associated with surface variable 556 948 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 557 & zprofvar1, & ! Model values for 1st variable in a prof ob 558 & zprofvar2 ! Model values for 2nd variable in a prof ob 949 & zglam, & ! Model longitudes for prof variables 950 & zgphi ! Model latitudes for prof variables 951 LOGICAL :: llog10 ! Perform log10 transform of variable 952 #if defined key_fabm 559 953 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 560 & zprofmask1, & ! Mask associated with zprofvar1 561 & zprofmask2 ! Mask associated with zprofvar2 562 REAL(wp), POINTER, DIMENSION(:,:) :: & 563 & zsurfvar ! Model values equivalent to surface ob. 564 REAL(wp), POINTER, DIMENSION(:,:) :: & 565 & zglam1, & ! Model longitudes for prof variable 1 566 & zglam2, & ! Model longitudes for prof variable 2 567 & zgphi1, & ! Model latitudes for prof variable 1 568 & zgphi2 ! Model latitudes for prof variable 2 569 #if ! defined key_lim2 && ! defined key_lim3 570 REAL(wp), POINTER, DIMENSION(:,:) :: frld 571 #endif 572 LOGICAL :: llnightav ! Logical for calculating night-time average 573 574 !Allocate local work arrays 575 CALL wrk_alloc( jpi, jpj, jpk, zprofvar1 ) 576 CALL wrk_alloc( jpi, jpj, jpk, zprofvar2 ) 577 CALL wrk_alloc( jpi, jpj, jpk, zprofmask1 ) 578 CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 ) 579 CALL wrk_alloc( jpi, jpj, zsurfvar ) 580 CALL wrk_alloc( jpi, jpj, zglam1 ) 581 CALL wrk_alloc( jpi, jpj, zglam2 ) 582 CALL wrk_alloc( jpi, jpj, zgphi1 ) 583 CALL wrk_alloc( jpi, jpj, zgphi2 ) 584 #if ! defined key_lim2 && ! defined key_lim3 585 CALL wrk_alloc(jpi,jpj,frld) 954 & pco2_3d ! 3D pCO2 from FABM 586 955 #endif 587 956 … … 595 964 596 965 !----------------------------------------------------------------------- 597 ! No LIM => frld == 0.0_wp598 !-----------------------------------------------------------------------599 #if ! defined key_lim2 && ! defined key_lim3600 frld(:,:) = 0.0_wp601 #endif602 !-----------------------------------------------------------------------603 966 ! Call the profile and surface observation operators 604 967 !----------------------------------------------------------------------- … … 608 971 DO jtype = 1, nproftypes 609 972 973 ! Allocate local work arrays 974 CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar ) 975 CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask ) 976 CALL wrk_alloc( jpi, jpj, profdataqc(jtype)%nvar, zglam ) 977 CALL wrk_alloc( jpi, jpj, profdataqc(jtype)%nvar, zgphi ) 978 979 ! Defaults which might change 980 DO jvar = 1, profdataqc(jtype)%nvar 981 zprofmask(:,:,:,jvar) = tmask(:,:,:) 982 zglam(:,:,jvar) = glamt(:,:) 983 zgphi(:,:,jvar) = gphit(:,:) 984 END DO 985 610 986 SELECT CASE ( TRIM(cobstypesprof(jtype)) ) 611 987 CASE('prof') 612 zprofvar1(:,:,:) = tsn(:,:,:,jp_tem) 613 zprofvar2(:,:,:) = tsn(:,:,:,jp_sal) 614 zprofmask1(:,:,:) = tmask(:,:,:) 615 zprofmask2(:,:,:) = tmask(:,:,:) 616 zglam1(:,:) = glamt(:,:) 617 zglam2(:,:) = glamt(:,:) 618 zgphi1(:,:) = gphit(:,:) 619 zgphi2(:,:) = gphit(:,:) 988 zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem) 989 zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal) 990 620 991 CASE('vel') 621 zprofvar1(:,:,:) = un(:,:,:) 622 zprofvar2(:,:,:) = vn(:,:,:) 623 zprofmask1(:,:,:) = umask(:,:,:) 624 zprofmask2(:,:,:) = vmask(:,:,:) 625 zglam1(:,:) = glamu(:,:) 626 zglam2(:,:) = glamv(:,:) 627 zgphi1(:,:) = gphiu(:,:) 628 zgphi2(:,:) = gphiv(:,:) 992 zprofvar(:,:,:,1) = un(:,:,:) 993 zprofvar(:,:,:,2) = vn(:,:,:) 994 zprofmask(:,:,:,1) = umask(:,:,:) 995 zprofmask(:,:,:,2) = vmask(:,:,:) 996 zglam(:,:,1) = glamu(:,:) 997 zglam(:,:,2) = glamv(:,:) 998 zgphi(:,:,1) = gphiu(:,:) 999 zgphi(:,:,2) = gphiv(:,:) 1000 1001 CASE('plchltot') 1002 #if defined key_hadocc 1003 ! Chlorophyll from HadOCC 1004 zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) 1005 #elif defined key_medusa 1006 ! Add non-diatom and diatom chlorophyll from MEDUSA 1007 zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 1008 #elif defined key_fabm 1009 ! Add all chlorophyll groups from ERSEM 1010 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + & 1011 & trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4) 1012 #else 1013 CALL ctl_stop( ' Trying to run plchltot observation operator', & 1014 & ' but no biogeochemical model appears to have been defined' ) 1015 #endif 1016 ! Take the log10 where we can, otherwise exclude 1017 tiny = 1.0e-20 1018 WHERE(zprofvar(:,:,:,:) > tiny .AND. zprofvar(:,:,:,:) /= obfillflt ) 1019 zprofvar(:,:,:,:) = LOG10(zprofvar(:,:,:,:)) 1020 ELSEWHERE 1021 zprofvar(:,:,:,:) = obfillflt 1022 zprofmask(:,:,:,:) = 0 1023 END WHERE 1024 ! Mask out model below any excluded values, 1025 ! to avoid interpolation issues 1026 DO jvar = 1, profdataqc(jtype)%nvar 1027 DO jj = 1, jpj 1028 DO ji = 1, jpi 1029 depth_loop: DO jk = 1, jpk 1030 IF ( zprofmask(ji,jj,jk,jvar) == 0 ) THEN 1031 zprofmask(ji,jj,jk:jpk,jvar) = 0 1032 EXIT depth_loop 1033 ENDIF 1034 END DO depth_loop 1035 END DO 1036 END DO 1037 END DO 1038 1039 CASE('pchltot') 1040 #if defined key_hadocc 1041 ! Chlorophyll from HadOCC 1042 zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) 1043 #elif defined key_medusa 1044 ! Add non-diatom and diatom chlorophyll from MEDUSA 1045 zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 1046 #elif defined key_fabm 1047 ! Add all chlorophyll groups from ERSEM 1048 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + & 1049 & trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4) 1050 #else 1051 CALL ctl_stop( ' Trying to run pchltot observation operator', & 1052 & ' but no biogeochemical model appears to have been defined' ) 1053 #endif 1054 1055 CASE('pno3') 1056 #if defined key_hadocc 1057 ! Dissolved inorganic nitrogen from HadOCC 1058 zprofvar(:,:,:,1) = trn(:,:,:,jp_had_nut) 1059 #elif defined key_medusa 1060 ! Dissolved inorganic nitrogen from MEDUSA 1061 zprofvar(:,:,:,1) = trn(:,:,:,jpdin) 1062 #elif defined key_fabm 1063 ! Nitrate from ERSEM 1064 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n) 1065 #else 1066 CALL ctl_stop( ' Trying to run pno3 observation operator', & 1067 & ' but no biogeochemical model appears to have been defined' ) 1068 #endif 1069 1070 CASE('psi4') 1071 #if defined key_hadocc 1072 CALL ctl_stop( ' Trying to run psi4 observation operator', & 1073 & ' but HadOCC does not simulate silicate' ) 1074 #elif defined key_medusa 1075 ! Silicate from MEDUSA 1076 zprofvar(:,:,:,1) = trn(:,:,:,jpsil) 1077 #elif defined key_fabm 1078 ! Silicate from ERSEM 1079 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n5s) 1080 #else 1081 CALL ctl_stop( ' Trying to run psi4 observation operator', & 1082 & ' but no biogeochemical model appears to have been defined' ) 1083 #endif 1084 1085 CASE('ppo4') 1086 #if defined key_hadocc 1087 CALL ctl_stop( ' Trying to run ppo4 observation operator', & 1088 & ' but HadOCC does not simulate phosphate' ) 1089 #elif defined key_medusa 1090 CALL ctl_stop( ' Trying to run ppo4 observation operator', & 1091 & ' but MEDUSA does not simulate phosphate' ) 1092 #elif defined key_fabm 1093 ! Phosphate from ERSEM 1094 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n1p) 1095 #else 1096 CALL ctl_stop( ' Trying to run ppo4 observation operator', & 1097 & ' but no biogeochemical model appears to have been defined' ) 1098 #endif 1099 1100 CASE('pdic') 1101 #if defined key_hadocc 1102 ! Dissolved inorganic carbon from HadOCC 1103 zprofvar(:,:,:,1) = trn(:,:,:,jp_had_dic) 1104 #elif defined key_medusa 1105 ! Dissolved inorganic carbon from MEDUSA 1106 zprofvar(:,:,:,1) = trn(:,:,:,jpdic) 1107 #elif defined key_fabm 1108 ! Dissolved inorganic carbon from ERSEM 1109 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3c) 1110 #else 1111 CALL ctl_stop( ' Trying to run pdic observation operator', & 1112 & ' but no biogeochemical model appears to have been defined' ) 1113 #endif 1114 1115 CASE('palk') 1116 #if defined key_hadocc 1117 ! Alkalinity from HadOCC 1118 zprofvar(:,:,:,1) = trn(:,:,:,jp_had_alk) 1119 #elif defined key_medusa 1120 ! Alkalinity from MEDUSA 1121 zprofvar(:,:,:,1) = trn(:,:,:,jpalk) 1122 #elif defined key_fabm 1123 ! Alkalinity from ERSEM 1124 zprofvar(:,:,:,1) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta) 1125 #else 1126 CALL ctl_stop( ' Trying to run palk observation operator', & 1127 & ' but no biogeochemical model appears to have been defined' ) 1128 #endif 1129 1130 CASE('pph') 1131 #if defined key_hadocc 1132 CALL ctl_stop( ' Trying to run pph observation operator', & 1133 & ' but HadOCC has no pH diagnostic defined' ) 1134 #elif defined key_medusa && defined key_roam 1135 ! pH from MEDUSA 1136 zprofvar(:,:,:,1) = f3_pH(:,:,:) 1137 #elif defined key_fabm 1138 ! pH from ERSEM 1139 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3ph) 1140 #else 1141 CALL ctl_stop( ' Trying to run pph observation operator', & 1142 & ' but no biogeochemical model appears to have been defined' ) 1143 #endif 1144 1145 CASE('po2') 1146 #if defined key_hadocc 1147 CALL ctl_stop( ' Trying to run po2 observation operator', & 1148 & ' but HadOCC does not simulate oxygen' ) 1149 #elif defined key_medusa 1150 ! Oxygen from MEDUSA 1151 zprofvar(:,:,:,1) = trn(:,:,:,jpoxy) 1152 #elif defined key_fabm 1153 ! Oxygen from ERSEM 1154 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o2o) 1155 #else 1156 CALL ctl_stop( ' Trying to run po2 observation operator', & 1157 & ' but no biogeochemical model appears to have been defined' ) 1158 #endif 1159 1160 CASE DEFAULT 1161 CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) 1162 629 1163 END SELECT 630 1164 631 IF( ln_zco .OR. ln_zps ) THEN1165 DO jvar = 1, profdataqc(jtype)%nvar 632 1166 CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk, & 633 & nit000, idaystp, & 634 & zprofvar1, zprofvar2, & 635 & gdept_1d, zprofmask1, zprofmask2, & 636 & zglam1, zglam2, zgphi1, zgphi2, & 637 & nn_1dint, nn_2dint, & 1167 & nit000, idaystp, jvar, & 1168 & zprofvar(:,:,:,jvar), & 1169 & gdept_n(:,:,:), gdepw_n(:,:,:), & 1170 & zprofmask(:,:,:,jvar), & 1171 & zglam(:,:,jvar), zgphi(:,:,jvar), & 1172 & nn_1dint, nn_2dint_default, & 638 1173 & kdailyavtypes = nn_profdavtypes ) 639 ELSE IF(TRIM(cobstypesprof(jtype)) == 'prof') THEN 640 !TG - THIS NEEDS MODIFICATION TO MATCH SIMPLIFICATION 641 CALL obs_pro_sco_opt( profdataqc(jtype), & 642 & kstp, jpi, jpj, jpk, nit000, idaystp, & 643 & zprofvar1, zprofvar2, & 644 & gdept_n(:,:,:), gdepw_n(:,:,:), & 645 & tmask, nn_1dint, nn_2dint, & 646 & kdailyavtypes = nn_profdavtypes ) 647 ELSE 648 CALL ctl_stop('DIA_OBS: Generalised vertical interpolation not'// & 649 'yet working for velocity data (turn off velocity observations') 650 ENDIF 1174 END DO 1175 1176 CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar ) 1177 CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask ) 1178 CALL wrk_dealloc( jpi, jpj, profdataqc(jtype)%nvar, zglam ) 1179 CALL wrk_dealloc( jpi, jpj, profdataqc(jtype)%nvar, zgphi ) 651 1180 652 1181 END DO … … 656 1185 IF ( nsurftypes > 0 ) THEN 657 1186 1187 !Allocate local work arrays 1188 CALL wrk_alloc( jpi, jpj, zsurfvar ) 1189 CALL wrk_alloc( jpi, jpj, zsurfmask ) 1190 #if defined key_fabm 1191 CALL wrk_alloc( jpi, jpj, jpk, pco2_3d ) 1192 #endif 1193 658 1194 DO jtype = 1, nsurftypes 1195 1196 !Defaults which might be changed 1197 zsurfmask(:,:) = tmask(:,:,1) 1198 llog10 = .FALSE. 659 1199 660 1200 SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 661 1201 CASE('sst') 662 1202 zsurfvar(:,:) = tsn(:,:,1,jp_tem) 663 llnightav = ln_sstnight1203 llnightav(jtype) = ln_sstnight 664 1204 CASE('sla') 665 1205 zsurfvar(:,:) = sshn(:,:) 666 llnightav = .FALSE.667 #if defined key_lim2 || defined key_lim3 1206 CASE('sss') 1207 zsurfvar(:,:) = tsn(:,:,1,jp_sal) 668 1208 CASE('sic') 669 1209 IF ( kstp == 0 ) THEN … … 678 1218 CYCLE 679 1219 ELSE 1220 #if defined key_cice 1221 zsurfvar(:,:) = fr_i(:,:) 1222 #elif defined key_lim2 || defined key_lim3 680 1223 zsurfvar(:,:) = 1._wp - frld(:,:) 1224 #else 1225 CALL ctl_stop( ' Trying to run sea-ice observation operator', & 1226 & ' but no sea-ice model appears to have been defined' ) 1227 #endif 681 1228 ENDIF 682 683 llnightav = .FALSE. 684 #endif 1229 CASE('slchltot') 1230 #if defined key_hadocc 1231 ! Surface chlorophyll from HadOCC 1232 zsurfvar(:,:) = HADOCC_CHL(:,:,1) 1233 #elif defined key_medusa 1234 ! Add non-diatom and diatom surface chlorophyll from MEDUSA 1235 zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 1236 #elif defined key_fabm 1237 ! Add all surface chlorophyll groups from ERSEM 1238 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 1239 & trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1240 #else 1241 CALL ctl_stop( ' Trying to run slchltot observation operator', & 1242 & ' but no biogeochemical model appears to have been defined' ) 1243 #endif 1244 llog10 = .TRUE. 1245 1246 CASE('slchldia') 1247 #if defined key_hadocc 1248 CALL ctl_stop( ' Trying to run slchldia observation operator', & 1249 & ' but HadOCC does not explicitly simulate diatoms' ) 1250 #elif defined key_medusa 1251 ! Diatom surface chlorophyll from MEDUSA 1252 zsurfvar(:,:) = trn(:,:,1,jpchd) 1253 #elif defined key_fabm 1254 ! Diatom surface chlorophyll from ERSEM 1255 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) 1256 #else 1257 CALL ctl_stop( ' Trying to run slchldia observation operator', & 1258 & ' but no biogeochemical model appears to have been defined' ) 1259 #endif 1260 llog10 = .TRUE. 1261 1262 CASE('slchlnon') 1263 #if defined key_hadocc 1264 CALL ctl_stop( ' Trying to run slchlnon observation operator', & 1265 & ' but HadOCC does not explicitly simulate non-diatoms' ) 1266 #elif defined key_medusa 1267 ! Non-diatom surface chlorophyll from MEDUSA 1268 zsurfvar(:,:) = trn(:,:,1,jpchn) 1269 #elif defined key_fabm 1270 ! Add all non-diatom surface chlorophyll groups from ERSEM 1271 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 1272 & trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1273 #else 1274 CALL ctl_stop( ' Trying to run slchlnon observation operator', & 1275 & ' but no biogeochemical model appears to have been defined' ) 1276 #endif 1277 llog10 = .TRUE. 1278 1279 CASE('slchldin') 1280 #if defined key_hadocc 1281 CALL ctl_stop( ' Trying to run slchldin observation operator', & 1282 & ' but HadOCC does not explicitly simulate dinoflagellates' ) 1283 #elif defined key_medusa 1284 CALL ctl_stop( ' Trying to run slchldin observation operator', & 1285 & ' but MEDUSA does not explicitly simulate dinoflagellates' ) 1286 #elif defined key_fabm 1287 ! Dinoflagellate surface chlorophyll from ERSEM 1288 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1289 #else 1290 CALL ctl_stop( ' Trying to run slchldin observation operator', & 1291 & ' but no biogeochemical model appears to have been defined' ) 1292 #endif 1293 llog10 = .TRUE. 1294 1295 CASE('slchlmic') 1296 #if defined key_hadocc 1297 CALL ctl_stop( ' Trying to run slchlmic observation operator', & 1298 & ' but HadOCC does not explicitly simulate microphytoplankton' ) 1299 #elif defined key_medusa 1300 CALL ctl_stop( ' Trying to run slchlmic observation operator', & 1301 & ' but MEDUSA does not explicitly simulate microphytoplankton' ) 1302 #elif defined key_fabm 1303 ! Add diatom and dinoflagellate surface chlorophyll from ERSEM 1304 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1305 #else 1306 CALL ctl_stop( ' Trying to run slchlmic observation operator', & 1307 & ' but no biogeochemical model appears to have been defined' ) 1308 #endif 1309 llog10 = .TRUE. 1310 1311 CASE('slchlnan') 1312 #if defined key_hadocc 1313 CALL ctl_stop( ' Trying to run slchlnan observation operator', & 1314 & ' but HadOCC does not explicitly simulate nanophytoplankton' ) 1315 #elif defined key_medusa 1316 CALL ctl_stop( ' Trying to run slchlnan observation operator', & 1317 & ' but MEDUSA does not explicitly simulate nanophytoplankton' ) 1318 #elif defined key_fabm 1319 ! Nanophytoplankton surface chlorophyll from ERSEM 1320 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) 1321 #else 1322 CALL ctl_stop( ' Trying to run slchlnan observation operator', & 1323 & ' but no biogeochemical model appears to have been defined' ) 1324 #endif 1325 llog10 = .TRUE. 1326 1327 CASE('slchlpic') 1328 #if defined key_hadocc 1329 CALL ctl_stop( ' Trying to run slchlpic observation operator', & 1330 & ' but HadOCC does not explicitly simulate picophytoplankton' ) 1331 #elif defined key_medusa 1332 CALL ctl_stop( ' Trying to run slchlpic observation operator', & 1333 & ' but MEDUSA does not explicitly simulate picophytoplankton' ) 1334 #elif defined key_fabm 1335 ! Picophytoplankton surface chlorophyll from ERSEM 1336 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) 1337 #else 1338 CALL ctl_stop( ' Trying to run slchlpic observation operator', & 1339 & ' but no biogeochemical model appears to have been defined' ) 1340 #endif 1341 llog10 = .TRUE. 1342 1343 CASE('schltot') 1344 #if defined key_hadocc 1345 ! Surface chlorophyll from HadOCC 1346 zsurfvar(:,:) = HADOCC_CHL(:,:,1) 1347 #elif defined key_medusa 1348 ! Add non-diatom and diatom surface chlorophyll from MEDUSA 1349 zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 1350 #elif defined key_fabm 1351 ! Add all surface chlorophyll groups from ERSEM 1352 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 1353 & trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1354 #else 1355 CALL ctl_stop( ' Trying to run schltot observation operator', & 1356 & ' but no biogeochemical model appears to have been defined' ) 1357 #endif 1358 1359 CASE('slphytot') 1360 #if defined key_hadocc 1361 ! Surface phytoplankton nitrogen from HadOCC multiplied by C:N ratio 1362 zsurfvar(:,:) = trn(:,:,1,jp_had_phy) * c2n_p 1363 #elif defined key_medusa 1364 ! Add non-diatom and diatom surface phytoplankton nitrogen from MEDUSA 1365 ! multiplied by C:N ratio for each 1366 zsurfvar(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd) 1367 #elif defined key_fabm 1368 ! Add all surface phytoplankton carbon groups from ERSEM 1369 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 1370 & trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c) 1371 #else 1372 CALL ctl_stop( ' Trying to run slphytot observation operator', & 1373 & ' but no biogeochemical model appears to have been defined' ) 1374 #endif 1375 llog10 = .TRUE. 1376 1377 CASE('slphydia') 1378 #if defined key_hadocc 1379 CALL ctl_stop( ' Trying to run slphydia observation operator', & 1380 & ' but HadOCC does not explicitly simulate diatoms' ) 1381 #elif defined key_medusa 1382 ! Diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio 1383 zsurfvar(:,:) = trn(:,:,1,jpphd) * xthetapd 1384 #elif defined key_fabm 1385 ! Diatom surface phytoplankton carbon from ERSEM 1386 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) 1387 #else 1388 CALL ctl_stop( ' Trying to run slphydia observation operator', & 1389 & ' but no biogeochemical model appears to have been defined' ) 1390 #endif 1391 llog10 = .TRUE. 1392 1393 CASE('slphynon') 1394 #if defined key_hadocc 1395 CALL ctl_stop( ' Trying to run slphynon observation operator', & 1396 & ' but HadOCC does not explicitly simulate non-diatoms' ) 1397 #elif defined key_medusa 1398 ! Non-diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio 1399 zsurfvar(:,:) = trn(:,:,1,jpphn) * xthetapn 1400 #elif defined key_fabm 1401 ! Add all non-diatom surface phytoplankton carbon groups from ERSEM 1402 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 1403 & trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c) 1404 #else 1405 CALL ctl_stop( ' Trying to run slphynon observation operator', & 1406 & ' but no biogeochemical model appears to have been defined' ) 1407 #endif 1408 llog10 = .TRUE. 1409 1410 CASE('sspm') 1411 #if defined key_spm 1412 zsurfvar(:,:) = 0.0 1413 DO jn = 1, jp_spm 1414 zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn) ! sum SPM sizes 1415 END DO 1416 #else 1417 CALL ctl_stop( ' Trying to run sspm observation operator', & 1418 & ' but no spm model appears to have been defined' ) 1419 #endif 1420 1421 CASE('sfco2') 1422 #if defined key_hadocc 1423 zsurfvar(:,:) = HADOCC_FCO2(:,:) ! fCO2 from HadOCC 1424 IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. & 1425 & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN 1426 zsurfvar(:,:) = obfillflt 1427 zsurfmask(:,:) = 0 1428 CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', & 1429 & ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' ) 1430 ENDIF 1431 #elif defined key_medusa && defined key_roam 1432 zsurfvar(:,:) = f2_fco2w(:,:) 1433 #elif defined key_fabm 1434 ! First, get pCO2 from FABM 1435 pco2_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc) 1436 zsurfvar(:,:) = pco2_3d(:,:,1) 1437 ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of: 1438 ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems 1439 ! and data reduction routines, Deep-Sea Research II, 56: 512-522. 1440 ! and 1441 ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas, 1442 ! Marine Chemistry, 2: 203-215. 1443 ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so 1444 ! not explicitly included - atmospheric pressure is not necessarily available so this is 1445 ! the best assumption. 1446 ! Further, the (1-xCO2)^2 term has been neglected. This is common practice 1447 ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes) 1448 ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal 1449 ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway. 1450 zsurfvar(:,:) = zsurfvar(:,:) * EXP((-1636.75 + & 1451 & 12.0408 * (tsn(:,:,1,jp_tem)+rt0) - & 1452 & 0.0327957 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & 1453 & 0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & 1454 & 2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0))) / & 1455 & (82.0578 * (tsn(:,:,1,jp_tem)+rt0))) 1456 #else 1457 CALL ctl_stop( ' Trying to run sfco2 observation operator', & 1458 & ' but no biogeochemical model appears to have been defined' ) 1459 #endif 1460 1461 CASE('spco2') 1462 #if defined key_hadocc 1463 zsurfvar(:,:) = HADOCC_PCO2(:,:) ! pCO2 from HadOCC 1464 IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. & 1465 & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN 1466 zsurfvar(:,:) = obfillflt 1467 zsurfmask(:,:) = 0 1468 CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', & 1469 & ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' ) 1470 ENDIF 1471 #elif defined key_medusa && defined key_roam 1472 zsurfvar(:,:) = f2_pco2w(:,:) 1473 #elif defined key_fabm 1474 pco2_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc) 1475 zsurfvar(:,:) = pco2_3d(:,:,1) 1476 #else 1477 CALL ctl_stop( ' Trying to run spco2 observation operator', & 1478 & ' but no biogeochemical model appears to have been defined' ) 1479 #endif 1480 1481 CASE DEFAULT 1482 1483 CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' ) 1484 685 1485 END SELECT 1486 1487 IF ( llog10 ) THEN 1488 ! Take the log10 where we can, otherwise exclude 1489 tiny = 1.0e-20 1490 WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt ) 1491 zsurfvar(:,:) = LOG10(zsurfvar(:,:)) 1492 ELSEWHERE 1493 zsurfvar(:,:) = obfillflt 1494 zsurfmask(:,:) = 0 1495 END WHERE 1496 ENDIF 686 1497 687 1498 CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj, & 688 & nit000, idaystp, zsurfvar, tmask(:,:,1), & 689 & nn_2dint, llnightav ) 1499 & nit000, idaystp, zsurfvar, zsurfmask, & 1500 & n2dintsurf(jtype), llnightav(jtype), & 1501 & ravglamscl(jtype), ravgphiscl(jtype), & 1502 & lfpindegs(jtype) ) 690 1503 691 1504 END DO 692 1505 693 ENDIF 694 695 CALL wrk_dealloc( jpi, jpj, jpk, zprofvar1 ) 696 CALL wrk_dealloc( jpi, jpj, jpk, zprofvar2 ) 697 CALL wrk_dealloc( jpi, jpj, jpk, zprofmask1 ) 698 CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 ) 699 CALL wrk_dealloc( jpi, jpj, zsurfvar ) 700 CALL wrk_dealloc( jpi, jpj, zglam1 ) 701 CALL wrk_dealloc( jpi, jpj, zglam2 ) 702 CALL wrk_dealloc( jpi, jpj, zgphi1 ) 703 CALL wrk_dealloc( jpi, jpj, zgphi2 ) 704 #if ! defined key_lim2 && ! defined key_lim3 705 CALL wrk_dealloc(jpi,jpj,frld) 706 #endif 1506 CALL wrk_dealloc( jpi, jpj, zsurfvar ) 1507 CALL wrk_dealloc( jpi, jpj, zsurfmask ) 1508 #if defined key_fabm 1509 CALL wrk_dealloc( jpi, jpj, jpk, pco2_3d ) 1510 #endif 1511 1512 ENDIF 707 1513 708 1514 END SUBROUTINE dia_obs … … 755 1561 & ) 756 1562 757 CALL obs_rotvel( profdataqc(jtype), nn_2dint , zu, zv )1563 CALL obs_rotvel( profdataqc(jtype), nn_2dint_default, zu, zv ) 758 1564 759 1565 DO jo = 1, profdataqc(jtype)%nprof … … 821 1627 822 1628 IF ( nsurftypes > 0 ) & 823 & DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf ) 1629 & DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, & 1630 & n2dintsurf, ravglamscl, ravgphiscl, lfpindegs, llnightav ) 824 1631 825 1632 END SUBROUTINE dia_obs_dealloc … … 970 1777 END SUBROUTINE fin_date 971 1778 1779 SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, ifiles, cobstypes, cfiles ) 1780 1781 INTEGER, INTENT(IN) :: ntypes ! Total number of obs types 1782 INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type 1783 INTEGER, DIMENSION(ntypes), INTENT(OUT) :: & 1784 & ifiles ! Out number of files for each type 1785 CHARACTER(len=8), DIMENSION(ntypes), INTENT(IN) :: & 1786 & cobstypes ! List of obs types 1787 CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(IN) :: & 1788 & cfiles ! List of files for all types 1789 1790 !Local variables 1791 INTEGER :: jfile 1792 INTEGER :: jtype 1793 1794 DO jtype = 1, ntypes 1795 1796 ifiles(jtype) = 0 1797 DO jfile = 1, jpmaxnfiles 1798 IF ( trim(cfiles(jtype,jfile)) /= '' ) & 1799 ifiles(jtype) = ifiles(jtype) + 1 1800 END DO 1801 1802 IF ( ifiles(jtype) == 0 ) THEN 1803 CALL ctl_stop( 'Logical for observation type '//TRIM(cobstypes(jtype))// & 1804 & ' set to true but no files available to read' ) 1805 ENDIF 1806 1807 IF(lwp) THEN 1808 WRITE(numout,*) ' '//cobstypes(jtype)//' input observation file names:' 1809 DO jfile = 1, ifiles(jtype) 1810 WRITE(numout,*) ' '//TRIM(cfiles(jtype,jfile)) 1811 END DO 1812 ENDIF 1813 1814 END DO 1815 1816 END SUBROUTINE obs_settypefiles 1817 1818 SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein, & 1819 & n2dint_default, n2dint_type, & 1820 & ravglamscl_type, ravgphiscl_type, & 1821 & lfp_indegs_type, lavnight_type, & 1822 & n2dint, ravglamscl, ravgphiscl, & 1823 & lfpindegs, lavnight ) 1824 1825 INTEGER, INTENT(IN) :: ntypes ! Total number of obs types 1826 INTEGER, INTENT(IN) :: jtype ! Index of the current type of obs 1827 INTEGER, INTENT(IN) :: n2dint_default ! Default option for interpolation type 1828 INTEGER, INTENT(IN) :: n2dint_type ! Option for interpolation type 1829 REAL(wp), INTENT(IN) :: & 1830 & ravglamscl_type, & !E/W diameter of obs footprint for this type 1831 & ravgphiscl_type !N/S diameter of obs footprint for this type 1832 LOGICAL, INTENT(IN) :: lfp_indegs_type !T=> footprint in degrees, F=> in metres 1833 LOGICAL, INTENT(IN) :: lavnight_type !T=> obs represent night time average 1834 CHARACTER(len=8), INTENT(IN) :: ctypein 1835 1836 INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 1837 & n2dint 1838 REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: & 1839 & ravglamscl, ravgphiscl 1840 LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: & 1841 & lfpindegs, lavnight 1842 1843 lavnight(jtype) = lavnight_type 1844 1845 IF ( (n2dint_type >= 0) .AND. (n2dint_type <= 6) ) THEN 1846 n2dint(jtype) = n2dint_type 1847 ELSE IF ( n2dint_type == -1 ) THEN 1848 n2dint(jtype) = n2dint_default 1849 ELSE 1850 CALL ctl_stop(' Choice of '//TRIM(ctypein)//' horizontal (2D) interpolation method', & 1851 & ' is not available') 1852 ENDIF 1853 1854 ! For averaging observation footprints set options for size of footprint 1855 IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN 1856 IF ( ravglamscl_type > 0._wp ) THEN 1857 ravglamscl(jtype) = ravglamscl_type 1858 ELSE 1859 CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 1860 'scale (ravglamscl) for observation type '//TRIM(ctypein) ) 1861 ENDIF 1862 1863 IF ( ravgphiscl_type > 0._wp ) THEN 1864 ravgphiscl(jtype) = ravgphiscl_type 1865 ELSE 1866 CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 1867 'scale (ravgphiscl) for observation type '//TRIM(ctypein) ) 1868 ENDIF 1869 1870 lfpindegs(jtype) = lfp_indegs_type 1871 1872 ENDIF 1873 1874 ! Write out info 1875 IF(lwp) THEN 1876 IF ( n2dint(jtype) <= 4 ) THEN 1877 WRITE(numout,*) ' '//TRIM(ctypein)// & 1878 & ' model counterparts will be interpolated horizontally' 1879 ELSE IF ( n2dint(jtype) <= 6 ) THEN 1880 WRITE(numout,*) ' '//TRIM(ctypein)// & 1881 & ' model counterparts will be averaged horizontally' 1882 WRITE(numout,*) ' '//' with E/W scale: ',ravglamscl(jtype) 1883 WRITE(numout,*) ' '//' with N/S scale: ',ravgphiscl(jtype) 1884 IF ( lfpindegs(jtype) ) THEN 1885 WRITE(numout,*) ' '//' (in degrees)' 1886 ELSE 1887 WRITE(numout,*) ' '//' (in metres)' 1888 ENDIF 1889 ENDIF 1890 ENDIF 1891 1892 END SUBROUTINE obs_setinterpopts 1893 972 1894 END MODULE diaobs -
NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_mpp.F90
r11350 r11361 98 98 ! 99 99 INTEGER :: ierr 100 INTEGER, DIMENSION(kno) :: ivals 101 ! 102 INCLUDE 'mpif.h' 103 !!---------------------------------------------------------------------- 100 INTEGER, DIMENSION(:), ALLOCATABLE :: ivals 101 ! 102 INCLUDE 'mpif.h' 103 !!---------------------------------------------------------------------- 104 105 ALLOCATE( ivals(kno) ) 104 106 105 107 ! Call the MPI library to find the maximum across processors … … 107 109 & mpi_max, mpi_comm_opa, ierr ) 108 110 kvals(:) = ivals(:) 111 112 DEALLOCATE( ivals ) 109 113 #else 110 114 ! no MPI: empty routine … … 138 142 ! 139 143 INTEGER :: ji, isum 140 INTEGER, DIMENSION(kno) :: iobsp 141 !! 142 !! 143 144 iobsp=kobsp 144 INTEGER, DIMENSION(:), ALLOCATABLE :: iobsp 145 !! 146 !! 147 148 ALLOCATE( iobsp(kno) ) 149 150 iobsp(:)=kobsp(:) 145 151 146 152 WHERE( iobsp(:) == -1 ) … … 148 154 END WHERE 149 155 150 iobsp =-1*iobsp156 iobsp(:)=-1*iobsp(:) 151 157 152 158 CALL obs_mpp_max_integer( iobsp, kno ) 153 159 154 kobsp =-1*iobsp160 kobsp(:)=-1*iobsp(:) 155 161 156 162 isum=0 … … 168 174 ENDIF 169 175 176 DEALLOCATE( iobsp ) 177 170 178 #else 171 179 ! no MPI: empty routine -
NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
r11350 r11361 9 9 !! obs_prof_opt : Compute the model counterpart of profile data 10 10 !! obs_surf_opt : Compute the model counterpart of surface data 11 !! obs_pro_sco_opt: Compute the model counterpart of temperature and12 !! salinity observations from profiles in generalised13 !! vertical coordinates14 11 !!---------------------------------------------------------------------- 15 12 … … 22 19 & obs_int_h2d, & 23 20 & obs_int_h2d_init 21 USE obs_averg_h2d, ONLY : & ! Horizontal averaging to the obs footprint 22 & obs_avg_h2d, & 23 & obs_avg_h2d_init, & 24 & obs_max_fpsize 24 25 USE obs_inter_z1d, ONLY : & ! Vertical interpolation to the obs pt 25 26 & obs_int_z1d, & 26 27 & obs_int_z1d_spl 27 USE obs_const, ONLY : &28 & obfillflt ! Fillvalue28 USE obs_const, ONLY : & ! Obs fill value 29 & obfillflt 29 30 USE dom_oce, ONLY : & 30 & glamt, glamu, glamv, & 31 & gphit, gphiu, gphiv, & 32 & gdept_n, gdept_0 33 USE lib_mpp, ONLY : & 31 & glamt, glamf, & 32 & gphit, gphif 33 USE lib_mpp, ONLY : & ! Warning and stopping routines 34 34 & ctl_warn, ctl_stop 35 USE sbcdcy, ONLY : & ! For calculation of where it is night-time 36 & sbc_dcy, nday_qsr 35 37 USE obs_grid, ONLY : & 36 38 & obs_level_search 37 USE sbcdcy, ONLY : & ! For calculation of where it is night-time38 & sbc_dcy, nday_qsr39 39 40 40 IMPLICIT NONE … … 44 44 45 45 PUBLIC obs_prof_opt, & ! Compute the model counterpart of profile obs 46 & obs_pro_sco_opt, & ! Compute the model counterpart of profile observations47 46 & obs_surf_opt ! Compute the model counterpart of surface obs 48 47 … … 58 57 CONTAINS 59 58 60 SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, & 61 & kit000, kdaystp, & 62 & pvar1, pvar2, pgdept, pmask1, pmask2, & 63 & plam1, plam2, pphi1, pphi2, & 59 60 SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, & 61 & kit000, kdaystp, kvar, & 62 & pvar, pgdept, pgdepw, & 63 & pmask, & 64 & plam, pphi, & 64 65 & k1dint, k2dint, kdailyavtypes ) 65 66 … … 111 112 !! ! 07-03 (K. Mogensen) General handling of profiles 112 113 !! ! 15-02 (M. Martin) Combined routine for all profile types 114 !! ! 17-02 (M. Martin) Include generalised vertical coordinate changes 113 115 !!----------------------------------------------------------------------- 114 116 … … 130 132 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 131 133 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 134 INTEGER, INTENT(IN) :: kvar ! Number of variable in prodatqc 132 135 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 133 & pvar1, & ! Model field 1 134 & pvar2, & ! Model field 2 135 & pmask1, & ! Land-sea mask 1 136 & pmask2 ! Land-sea mask 2 136 & pvar, & ! Model field for variable 137 & pmask ! Land-sea mask for variable 137 138 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 138 & plam1, & ! Model longitudes for variable 1 139 & plam2, & ! Model longitudes for variable 2 140 & pphi1, & ! Model latitudes for variable 1 141 & pphi2 ! Model latitudes for variable 2 142 REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 143 & pgdept ! Model array of depth levels 139 & plam, & ! Model longitudes for variable 140 & pphi ! Model latitudes for variable 141 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 142 & pgdept, & ! Model array of depth T levels 143 & pgdepw ! Model array of depth W levels 144 144 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 145 145 & kdailyavtypes ! Types for daily averages … … 156 156 INTEGER :: iend 157 157 INTEGER :: iobs 158 INTEGER :: iin, ijn, ikn, ik ! looping indices over interpolation nodes 159 INTEGER :: inum_obs 158 160 INTEGER, DIMENSION(imaxavtypes) :: & 159 161 & idailyavtypes 160 162 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 161 & igrdi 1, &162 & igrd i2, &163 & igrdj1, &164 & igrdj2 163 & igrdi, & 164 & igrdj 165 INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic 166 165 167 REAL(KIND=wp) :: zlam 166 168 REAL(KIND=wp) :: zphi 167 169 REAL(KIND=wp) :: zdaystp 168 170 REAL(KIND=wp), DIMENSION(kpk) :: & 169 & zobsmask1, &170 & zobsmask2, &171 171 & zobsk, & 172 172 & zobs2k 173 REAL(KIND=wp), DIMENSION(2,2, kpk) :: &173 REAL(KIND=wp), DIMENSION(2,2,1) :: & 174 174 & zweig1, & 175 & zweig 2175 & zweig 176 176 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 177 & zmask1, & 178 & zmask2, & 179 & zint1, & 180 & zint2, & 181 & zinm1, & 182 & zinm2 177 & zmask, & 178 & zint, & 179 & zinm, & 180 & zgdept, & 181 & zgdepw 183 182 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 184 & zglam1, & 185 & zglam2, & 186 & zgphi1, & 187 & zgphi2 183 & zglam, & 184 & zgphi 185 REAL(KIND=wp), DIMENSION(1) :: zmsk 186 REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 187 188 188 LOGICAL :: ld_dailyav 189 189 … … 215 215 DO jj = 1, jpj 216 216 DO ji = 1, jpi 217 prodatqc%vdmean(ji,jj,jk,1) = 0.0 218 prodatqc%vdmean(ji,jj,jk,2) = 0.0 217 prodatqc%vdmean(ji,jj,jk,kvar) = 0.0 219 218 END DO 220 219 END DO … … 225 224 DO jj = 1, jpj 226 225 DO ji = 1, jpi 227 ! Increment field 1 for computing daily mean 228 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 229 & + pvar1(ji,jj,jk) 230 ! Increment field 2 for computing daily mean 231 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 232 & + pvar2(ji,jj,jk) 226 ! Increment field for computing daily mean 227 prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) & 228 & + pvar(ji,jj,jk) 233 229 END DO 234 230 END DO … … 243 239 DO jj = 1, jpj 244 240 DO ji = 1, jpi 245 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 246 & * zdaystp 247 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 248 & * zdaystp 241 prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) & 242 & * zdaystp 249 243 END DO 250 244 END DO … … 256 250 ! Get the data for interpolation 257 251 ALLOCATE( & 258 & igrdi1(2,2,ipro), & 259 & igrdi2(2,2,ipro), & 260 & igrdj1(2,2,ipro), & 261 & igrdj2(2,2,ipro), & 262 & zglam1(2,2,ipro), & 263 & zglam2(2,2,ipro), & 264 & zgphi1(2,2,ipro), & 265 & zgphi2(2,2,ipro), & 266 & zmask1(2,2,kpk,ipro), & 267 & zmask2(2,2,kpk,ipro), & 268 & zint1(2,2,kpk,ipro), & 269 & zint2(2,2,kpk,ipro) & 252 & igrdi(2,2,ipro), & 253 & igrdj(2,2,ipro), & 254 & zglam(2,2,ipro), & 255 & zgphi(2,2,ipro), & 256 & zmask(2,2,kpk,ipro), & 257 & zint(2,2,kpk,ipro), & 258 & zgdept(2,2,kpk,ipro), & 259 & zgdepw(2,2,kpk,ipro) & 270 260 & ) 271 261 272 262 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 273 263 iobs = jobs - prodatqc%nprofup 274 igrdi1(1,1,iobs) = prodatqc%mi(jobs,1)-1 275 igrdj1(1,1,iobs) = prodatqc%mj(jobs,1)-1 276 igrdi1(1,2,iobs) = prodatqc%mi(jobs,1)-1 277 igrdj1(1,2,iobs) = prodatqc%mj(jobs,1) 278 igrdi1(2,1,iobs) = prodatqc%mi(jobs,1) 279 igrdj1(2,1,iobs) = prodatqc%mj(jobs,1)-1 280 igrdi1(2,2,iobs) = prodatqc%mi(jobs,1) 281 igrdj1(2,2,iobs) = prodatqc%mj(jobs,1) 282 igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1 283 igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1 284 igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1 285 igrdj2(1,2,iobs) = prodatqc%mj(jobs,2) 286 igrdi2(2,1,iobs) = prodatqc%mi(jobs,2) 287 igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1 288 igrdi2(2,2,iobs) = prodatqc%mi(jobs,2) 289 igrdj2(2,2,iobs) = prodatqc%mj(jobs,2) 264 igrdi(1,1,iobs) = prodatqc%mi(jobs,kvar)-1 265 igrdj(1,1,iobs) = prodatqc%mj(jobs,kvar)-1 266 igrdi(1,2,iobs) = prodatqc%mi(jobs,kvar)-1 267 igrdj(1,2,iobs) = prodatqc%mj(jobs,kvar) 268 igrdi(2,1,iobs) = prodatqc%mi(jobs,kvar) 269 igrdj(2,1,iobs) = prodatqc%mj(jobs,kvar)-1 270 igrdi(2,2,iobs) = prodatqc%mi(jobs,kvar) 271 igrdj(2,2,iobs) = prodatqc%mj(jobs,kvar) 290 272 END DO 291 273 292 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 ) 293 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 ) 294 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pmask1, zmask1 ) 295 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pvar1, zint1 ) 296 297 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, plam2, zglam2 ) 298 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, pphi2, zgphi2 ) 299 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pmask2, zmask2 ) 300 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2, zint2 ) 274 ! Initialise depth arrays 275 zgdept(:,:,:,:) = 0.0 276 zgdepw(:,:,:,:) = 0.0 277 278 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, plam, zglam ) 279 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, pphi, zgphi ) 280 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pmask, zmask ) 281 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pvar, zint ) 282 283 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept, zgdept ) 284 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw, zgdepw ) 301 285 302 286 ! At the end of the day also get interpolated means 303 287 IF ( ld_dailyav .AND. idayend == 0 ) THEN 304 288 305 ALLOCATE( & 306 & zinm1(2,2,kpk,ipro), & 307 & zinm2(2,2,kpk,ipro) & 308 & ) 309 310 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, & 311 & prodatqc%vdmean(:,:,:,1), zinm1 ) 312 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, & 313 & prodatqc%vdmean(:,:,:,2), zinm2 ) 289 ALLOCATE( zinm(2,2,kpk,ipro) ) 290 291 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, & 292 & prodatqc%vdmean(:,:,:,kvar), zinm ) 314 293 315 294 ENDIF 295 296 ! Return if no observations to process 297 ! Has to be done after comm commands to ensure processors 298 ! stay in sync 299 IF ( ipro == 0 ) RETURN 316 300 317 301 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro … … 339 323 zphi = prodatqc%rphi(jobs) 340 324 341 ! Horizontal weights and vertical mask 342 343 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 344 345 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, & 346 & zglam1(:,:,iobs), zgphi1(:,:,iobs), & 347 & zmask1(:,:,:,iobs), zweig1, zobsmask1 ) 348 349 ENDIF 350 351 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 352 353 CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi, & 354 & zglam2(:,:,iobs), zgphi2(:,:,iobs), & 355 & zmask2(:,:,:,iobs), zweig2, zobsmask2 ) 356 357 ENDIF 358 359 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 325 ! Horizontal weights 326 ! Masked values are calculated later. 327 IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN 328 329 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 330 & zglam(:,:,iobs), zgphi(:,:,iobs), & 331 & zmask(:,:,1,iobs), zweig1, zmsk ) 332 333 ENDIF 334 335 IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN 360 336 361 337 zobsk(:) = obfillflt … … 365 341 IF ( idayend == 0 ) THEN 366 342 ! Daily averaged data 367 CALL obs_int_h2d( kpk, kpk, & 368 & zweig1, zinm1(:,:,:,iobs), zobsk ) 369 370 ENDIF 371 372 ELSE 373 374 ! Point data 375 CALL obs_int_h2d( kpk, kpk, & 376 & zweig1, zint1(:,:,:,iobs), zobsk ) 377 378 ENDIF 379 380 !------------------------------------------------------------- 381 ! Compute vertical second-derivative of the interpolating 382 ! polynomial at obs points 383 !------------------------------------------------------------- 384 385 IF ( k1dint == 1 ) THEN 386 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 387 & pgdept, zobsmask1 ) 388 ENDIF 389 390 !----------------------------------------------------------------- 391 ! Vertical interpolation to the observation point 392 !----------------------------------------------------------------- 393 ista = prodatqc%npvsta(jobs,1) 394 iend = prodatqc%npvend(jobs,1) 395 CALL obs_int_z1d( kpk, & 396 & prodatqc%var(1)%mvk(ista:iend), & 397 & k1dint, iend - ista + 1, & 398 & prodatqc%var(1)%vdep(ista:iend), & 399 & zobsk, zobs2k, & 400 & prodatqc%var(1)%vmod(ista:iend), & 401 & pgdept, zobsmask1 ) 402 403 ENDIF 404 405 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 406 407 zobsk(:) = obfillflt 408 409 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 410 411 IF ( idayend == 0 ) THEN 412 413 ! Daily averaged data 414 CALL obs_int_h2d( kpk, kpk, & 415 & zweig2, zinm2(:,:,:,iobs), zobsk ) 416 417 ENDIF 418 419 ELSE 420 421 ! Point data 422 CALL obs_int_h2d( kpk, kpk, & 423 & zweig2, zint2(:,:,:,iobs), zobsk ) 424 425 ENDIF 426 427 428 !------------------------------------------------------------- 429 ! Compute vertical second-derivative of the interpolating 430 ! polynomial at obs points 431 !------------------------------------------------------------- 432 433 IF ( k1dint == 1 ) THEN 434 CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 435 & pgdept, zobsmask2 ) 436 ENDIF 437 438 !---------------------------------------------------------------- 439 ! Vertical interpolation to the observation point 440 !---------------------------------------------------------------- 441 ista = prodatqc%npvsta(jobs,2) 442 iend = prodatqc%npvend(jobs,2) 443 CALL obs_int_z1d( kpk, & 444 & prodatqc%var(2)%mvk(ista:iend),& 445 & k1dint, iend - ista + 1, & 446 & prodatqc%var(2)%vdep(ista:iend),& 447 & zobsk, zobs2k, & 448 & prodatqc%var(2)%vmod(ista:iend),& 449 & pgdept, zobsmask2 ) 450 451 ENDIF 452 453 END DO 454 455 ! Deallocate the data for interpolation 456 DEALLOCATE( & 457 & igrdi1, & 458 & igrdi2, & 459 & igrdj1, & 460 & igrdj2, & 461 & zglam1, & 462 & zglam2, & 463 & zgphi1, & 464 & zgphi2, & 465 & zmask1, & 466 & zmask2, & 467 & zint1, & 468 & zint2 & 469 & ) 470 471 ! At the end of the day also get interpolated means 472 IF ( ld_dailyav .AND. idayend == 0 ) THEN 473 DEALLOCATE( & 474 & zinm1, & 475 & zinm2 & 476 & ) 477 ENDIF 478 479 prodatqc%nprofup = prodatqc%nprofup + ipro 480 481 END SUBROUTINE obs_prof_opt 482 483 SUBROUTINE obs_pro_sco_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 484 & ptn, psn, pgdept, pgdepw, ptmask, k1dint, k2dint, & 485 & kdailyavtypes ) 486 !!----------------------------------------------------------------------- 487 !! 488 !! *** ROUTINE obs_pro_opt *** 489 !! 490 !! ** Purpose : Compute the model counterpart of profiles 491 !! data by interpolating from the model grid to the 492 !! observation point. Generalised vertical coordinate version 493 !! 494 !! ** Method : Linearly interpolate to each observation point using 495 !! the model values at the corners of the surrounding grid box. 496 !! 497 !! First, model values on the model grid are interpolated vertically to the 498 !! Depths of the profile observations. Two vertical interpolation schemes are 499 !! available: 500 !! - linear (k1dint = 0) 501 !! - Cubic spline (k1dint = 1) 502 !! 503 !! 504 !! Secondly the interpolated values are interpolated horizontally to the 505 !! obs (lon, lat) point. 506 !! Several horizontal interpolation schemes are available: 507 !! - distance-weighted (great circle) (k2dint = 0) 508 !! - distance-weighted (small angle) (k2dint = 1) 509 !! - bilinear (geographical grid) (k2dint = 2) 510 !! - bilinear (quadrilateral grid) (k2dint = 3) 511 !! - polynomial (quadrilateral grid) (k2dint = 4) 512 !! 513 !! For the cubic spline the 2nd derivative of the interpolating 514 !! polynomial is computed before entering the vertical interpolation 515 !! routine. 516 !! 517 !! For ENACT moored buoy data (e.g., TAO), the model equivalent is 518 !! a daily mean model temperature field. So, we first compute 519 !! the mean, then interpolate only at the end of the day. 520 !! 521 !! This is the procedure to be used with generalised vertical model 522 !! coordinates (ie s-coordinates. It is ~4x slower than the equivalent 523 !! horizontal then vertical interpolation algorithm, but can deal with situations 524 !! where the model levels are not flat. 525 !! ONLY PERFORMED if ln_sco=.TRUE. 526 !! 527 !! Note: the in situ temperature observations must be converted 528 !! to potential temperature (the model variable) prior to 529 !! assimilation. 530 !!?????????????????????????????????????????????????????????????? 531 !! INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR??? 532 !!?????????????????????????????????????????????????????????????? 533 !! 534 !! ** Action : 535 !! 536 !! History : 537 !! ! 2014-08 (J. While) Adapted from obs_pro_opt to handel generalised 538 !! vertical coordinates 539 !!----------------------------------------------------------------------- 540 541 !! * Modules used 542 USE obs_profiles_def ! Definition of storage space for profile obs. 543 544 IMPLICIT NONE 545 546 !! * Arguments 547 TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening 548 INTEGER, INTENT(IN) :: kt ! Time step 549 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 550 INTEGER, INTENT(IN) :: kpj 551 INTEGER, INTENT(IN) :: kpk 552 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 553 ! (kit000-1 = restart time) 554 INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header) 555 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 556 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 557 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 558 & ptn, & ! Model temperature field 559 & psn, & ! Model salinity field 560 & ptmask ! Land-sea mask 561 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 562 & pgdept, & ! Model array of depth T levels 563 & pgdepw ! Model array of depth W levels 564 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 565 & kdailyavtypes ! Types for daily averages 566 567 !! * Local declarations 568 INTEGER :: ji 569 INTEGER :: jj 570 INTEGER :: jk 571 INTEGER :: iico, ijco 572 INTEGER :: jobs 573 INTEGER :: inrc 574 INTEGER :: ipro 575 INTEGER :: idayend 576 INTEGER :: ista 577 INTEGER :: iend 578 INTEGER :: iobs 579 INTEGER :: iin, ijn, ikn, ik ! looping indices over interpolation nodes 580 INTEGER, DIMENSION(imaxavtypes) :: & 581 & idailyavtypes 582 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 583 & igrdi, & 584 & igrdj 585 INTEGER :: & 586 & inum_obs 587 INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic 588 REAL(KIND=wp) :: zlam 589 REAL(KIND=wp) :: zphi 590 REAL(KIND=wp) :: zdaystp 591 REAL(KIND=wp), DIMENSION(kpk) :: & 592 & zobsmask, & 593 & zobsk, & 594 & zobs2k 595 REAL(KIND=wp), DIMENSION(2,2,1) :: & 596 & zweig, & 597 & l_zweig 598 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 599 & zmask, & 600 & zintt, & 601 & zints, & 602 & zinmt, & 603 & zgdept,& 604 & zgdepw,& 605 & zinms 606 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 607 & zglam, & 608 & zgphi 609 REAL(KIND=wp), DIMENSION(1) :: zmsk_1 610 REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 611 612 !------------------------------------------------------------------------ 613 ! Local initialization 614 !------------------------------------------------------------------------ 615 ! ... Record and data counters 616 inrc = kt - kit000 + 2 617 ipro = prodatqc%npstp(inrc) 618 619 ! Daily average types 620 IF ( PRESENT(kdailyavtypes) ) THEN 621 idailyavtypes(:) = kdailyavtypes(:) 622 ELSE 623 idailyavtypes(:) = -1 624 ENDIF 625 626 ! Initialize daily mean for first time-step 627 idayend = MOD( kt - kit000 + 1, kdaystp ) 628 629 ! Added kt == 0 test to catch restart case 630 IF ( idayend == 1 .OR. kt == 0) THEN 631 632 IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt 633 DO jk = 1, jpk 634 DO jj = 1, jpj 635 DO ji = 1, jpi 636 prodatqc%vdmean(ji,jj,jk,1) = 0.0 637 prodatqc%vdmean(ji,jj,jk,2) = 0.0 638 END DO 639 END DO 640 END DO 641 642 ENDIF 643 644 DO jk = 1, jpk 645 DO jj = 1, jpj 646 DO ji = 1, jpi 647 ! Increment the temperature field for computing daily mean 648 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 649 & + ptn(ji,jj,jk) 650 ! Increment the salinity field for computing daily mean 651 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 652 & + psn(ji,jj,jk) 653 END DO 654 END DO 655 END DO 656 657 ! Compute the daily mean at the end of day 658 zdaystp = 1.0 / REAL( kdaystp ) 659 IF ( idayend == 0 ) THEN 660 DO jk = 1, jpk 661 DO jj = 1, jpj 662 DO ji = 1, jpi 663 prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 664 & * zdaystp 665 prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 666 & * zdaystp 667 END DO 668 END DO 669 END DO 670 ENDIF 671 672 ! Get the data for interpolation 673 ALLOCATE( & 674 & igrdi(2,2,ipro), & 675 & igrdj(2,2,ipro), & 676 & zglam(2,2,ipro), & 677 & zgphi(2,2,ipro), & 678 & zmask(2,2,kpk,ipro), & 679 & zintt(2,2,kpk,ipro), & 680 & zints(2,2,kpk,ipro), & 681 & zgdept(2,2,kpk,ipro), & 682 & zgdepw(2,2,kpk,ipro) & 683 & ) 684 685 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 686 iobs = jobs - prodatqc%nprofup 687 igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1 688 igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1 689 igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1 690 igrdj(1,2,iobs) = prodatqc%mj(jobs,1) 691 igrdi(2,1,iobs) = prodatqc%mi(jobs,1) 692 igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1 693 igrdi(2,2,iobs) = prodatqc%mi(jobs,1) 694 igrdj(2,2,iobs) = prodatqc%mj(jobs,1) 695 END DO 696 697 ! Initialise depth arrays 698 zgdept = 0.0 699 zgdepw = 0.0 700 701 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, glamt, zglam ) 702 CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, gphit, zgphi ) 703 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, ptmask,zmask ) 704 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, ptn, zintt ) 705 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, psn, zints ) 706 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept(:,:,:), & 707 & zgdept ) 708 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw(:,:,:), & 709 & zgdepw ) 710 711 ! At the end of the day also get interpolated means 712 IF ( idayend == 0 ) THEN 713 714 ALLOCATE( & 715 & zinmt(2,2,kpk,ipro), & 716 & zinms(2,2,kpk,ipro) & 717 & ) 718 719 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, & 720 & prodatqc%vdmean(:,:,:,1), zinmt ) 721 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, & 722 & prodatqc%vdmean(:,:,:,2), zinms ) 723 724 ENDIF 725 726 ! Return if no observations to process 727 ! Has to be done after comm commands to ensure processors 728 ! stay in sync 729 IF ( ipro == 0 ) RETURN 730 731 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 732 733 iobs = jobs - prodatqc%nprofup 734 735 IF ( kt /= prodatqc%mstp(jobs) ) THEN 736 737 IF(lwp) THEN 738 WRITE(numout,*) 739 WRITE(numout,*) ' E R R O R : Observation', & 740 & ' time step is not consistent with the', & 741 & ' model time step' 742 WRITE(numout,*) ' =========' 743 WRITE(numout,*) 744 WRITE(numout,*) ' Record = ', jobs, & 745 & ' kt = ', kt, & 746 & ' mstp = ', prodatqc%mstp(jobs), & 747 & ' ntyp = ', prodatqc%ntyp(jobs) 748 ENDIF 749 CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' ) 750 ENDIF 751 752 zlam = prodatqc%rlam(jobs) 753 zphi = prodatqc%rphi(jobs) 754 755 ! Horizontal weights 756 ! Only calculated once, for both T and S. 757 ! Masked values are calculated later. 758 759 IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. & 760 & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN 761 762 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 763 & zglam(:,:,iobs), zgphi(:,:,iobs), & 764 & zmask(:,:,1,iobs), zweig, zmsk_1 ) 765 766 ENDIF 767 768 ! IF zmsk_1 = 0; then ob is on land 769 IF (zmsk_1(1) < 0.1) THEN 770 WRITE(numout,*) 'WARNING (obs_oper) :- profile found within landmask' 771 772 ELSE 773 774 ! Temperature 775 776 IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 777 778 zobsk(:) = obfillflt 779 780 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 781 782 IF ( idayend == 0 ) THEN 783 784 ! Daily averaged moored buoy (MRB) data 785 786 ! vertically interpolate all 4 corners 787 ista = prodatqc%npvsta(jobs,1) 788 iend = prodatqc%npvend(jobs,1) 789 inum_obs = iend - ista + 1 790 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 791 792 DO iin=1,2 793 DO ijn=1,2 794 795 796 797 IF ( k1dint == 1 ) THEN 798 CALL obs_int_z1d_spl( kpk, & 799 & zinmt(iin,ijn,:,iobs), & 800 & zobs2k, zgdept(iin,ijn,:,iobs), & 801 & zmask(iin,ijn,:,iobs)) 802 ENDIF 803 804 CALL obs_level_search(kpk, & 805 & zgdept(iin,ijn,:,iobs), & 806 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 807 & iv_indic) 808 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 809 & prodatqc%var(1)%vdep(ista:iend), & 810 & zinmt(iin,ijn,:,iobs), & 811 & zobs2k, interp_corner(iin,ijn,:), & 812 & zgdept(iin,ijn,:,iobs), & 813 & zmask(iin,ijn,:,iobs)) 814 815 ENDDO 816 ENDDO 817 818 819 ELSE 820 821 CALL ctl_stop( ' A nonzero' // & 822 & ' number of profile T BUOY data should' // & 823 & ' only occur at the end of a given day' ) 824 825 ENDIF 826 827 ELSE 828 829 ! Point data 830 343 831 344 ! vertically interpolate all 4 corners 832 ista = prodatqc%npvsta(jobs, 1)833 iend = prodatqc%npvend(jobs, 1)345 ista = prodatqc%npvsta(jobs,kvar) 346 iend = prodatqc%npvend(jobs,kvar) 834 347 inum_obs = iend - ista + 1 835 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 836 DO iin=1,2 348 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 349 350 DO iin=1,2 837 351 DO ijn=1,2 838 839 352 840 353 IF ( k1dint == 1 ) THEN 841 354 CALL obs_int_z1d_spl( kpk, & 842 & zintt(iin,ijn,:,iobs),& 843 & zobs2k, zgdept(iin,ijn,:,iobs), & 844 & zmask(iin,ijn,:,iobs)) 845 355 & zinm(iin,ijn,:,iobs), & 356 & zobs2k, zgdept(iin,ijn,:,iobs), & 357 & zmask(iin,ijn,:,iobs)) 846 358 ENDIF 847 359 848 360 CALL obs_level_search(kpk, & 849 & zgdept(iin,ijn,:,iobs),& 850 & inum_obs, prodatqc%var(1)%vdep(ista:iend), & 851 & iv_indic) 852 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 853 & prodatqc%var(1)%vdep(ista:iend), & 854 & zintt(iin,ijn,:,iobs), & 855 & zobs2k,interp_corner(iin,ijn,:), & 856 & zgdept(iin,ijn,:,iobs), & 857 & zmask(iin,ijn,:,iobs) ) 858 361 & zgdept(iin,ijn,:,iobs), & 362 & inum_obs, prodatqc%var(kvar)%vdep(ista:iend), & 363 & iv_indic) 364 365 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 366 & prodatqc%var(kvar)%vdep(ista:iend), & 367 & zinm(iin,ijn,:,iobs), & 368 & zobs2k, interp_corner(iin,ijn,:), & 369 & zgdept(iin,ijn,:,iobs), & 370 & zmask(iin,ijn,:,iobs)) 371 859 372 ENDDO 860 373 ENDDO 374 375 ENDIF !idayend 376 377 ELSE 378 379 ! Point data 380 381 ! vertically interpolate all 4 corners 382 ista = prodatqc%npvsta(jobs,kvar) 383 iend = prodatqc%npvend(jobs,kvar) 384 inum_obs = iend - ista + 1 385 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 386 DO iin=1,2 387 DO ijn=1,2 388 389 IF ( k1dint == 1 ) THEN 390 CALL obs_int_z1d_spl( kpk, & 391 & zint(iin,ijn,:,iobs),& 392 & zobs2k, zgdept(iin,ijn,:,iobs), & 393 & zmask(iin,ijn,:,iobs)) 394 395 ENDIF 396 397 CALL obs_level_search(kpk, & 398 & zgdept(iin,ijn,:,iobs),& 399 & inum_obs, prodatqc%var(kvar)%vdep(ista:iend), & 400 & iv_indic) 401 402 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 403 & prodatqc%var(kvar)%vdep(ista:iend), & 404 & zint(iin,ijn,:,iobs), & 405 & zobs2k,interp_corner(iin,ijn,:), & 406 & zgdept(iin,ijn,:,iobs), & 407 & zmask(iin,ijn,:,iobs) ) 408 409 ENDDO 410 ENDDO 861 411 862 863 864 865 866 412 ENDIF 413 414 !------------------------------------------------------------- 415 ! Compute the horizontal interpolation for every profile level 416 !------------------------------------------------------------- 867 417 868 869 iend=ista+ikn-1870 871 l_zweig(:,:,1) = 0._wp872 873 874 875 876 ! the mask. This is important for observations arenear877 878 879 418 DO ikn=1,inum_obs 419 iend=ista+ikn-1 420 421 zweig(:,:,1) = 0._wp 422 423 ! This code forces the horizontal weights to be 424 ! zero IF the observation is below the bottom of the 425 ! corners of the interpolation nodes, Or if it is in 426 ! the mask. This is important for observations near 427 ! steep bathymetry 428 DO iin=1,2 429 DO ijn=1,2 880 430 881 depth_loop1: DO ik=kpk,2,-1882 431 depth_loop: DO ik=kpk,2,-1 432 IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN 883 433 884 l_zweig(iin,ijn,1) = &885 & zweig(iin,ijn,1) * &886 887 & - prodatqc%var(1)%vdep(iend)),0._wp)434 zweig(iin,ijn,1) = & 435 & zweig1(iin,ijn,1) * & 436 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 437 & - prodatqc%var(kvar)%vdep(iend)),0._wp) 888 438 889 EXIT depth_loop1 890 ENDIF 891 ENDDO depth_loop1 439 EXIT depth_loop 440 441 ENDIF 442 443 ENDDO depth_loop 892 444 893 ENDDO894 445 ENDDO 446 ENDDO 895 447 896 CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), & 897 & prodatqc%var(1)%vmod(iend:iend) ) 448 CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 449 & prodatqc%var(kvar)%vmod(iend:iend) ) 450 451 ! Set QC flag for any observations found below the bottom 452 ! needed as the check here is more strict than that in obs_prep 453 IF (sum(zweig) == 0.0_wp) prodatqc%var(kvar)%nvqc(iend:iend)=4 898 454 899 455 ENDDO 900 456 901 902 DEALLOCATE(interp_corner,iv_indic) 457 DEALLOCATE(interp_corner,iv_indic) 903 458 904 ENDIF 905 906 907 ! Salinity 908 909 IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 910 911 zobsk(:) = obfillflt 912 913 IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 914 915 IF ( idayend == 0 ) THEN 916 917 ! Daily averaged moored buoy (MRB) data 918 919 ! vertically interpolate all 4 corners 920 ista = prodatqc%npvsta(iobs,2) 921 iend = prodatqc%npvend(iobs,2) 922 inum_obs = iend - ista + 1 923 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 924 925 DO iin=1,2 926 DO ijn=1,2 927 928 929 930 IF ( k1dint == 1 ) THEN 931 CALL obs_int_z1d_spl( kpk, & 932 & zinms(iin,ijn,:,iobs), & 933 & zobs2k, zgdept(iin,ijn,:,iobs), & 934 & zmask(iin,ijn,:,iobs)) 935 ENDIF 936 937 CALL obs_level_search(kpk, & 938 & zgdept(iin,ijn,:,iobs), & 939 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 940 & iv_indic) 941 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 942 & prodatqc%var(2)%vdep(ista:iend), & 943 & zinms(iin,ijn,:,iobs), & 944 & zobs2k, interp_corner(iin,ijn,:), & 945 & zgdept(iin,ijn,:,iobs), & 946 & zmask(iin,ijn,:,iobs)) 947 948 ENDDO 949 ENDDO 950 951 952 ELSE 953 954 CALL ctl_stop( ' A nonzero' // & 955 & ' number of profile T BUOY data should' // & 956 & ' only occur at the end of a given day' ) 957 958 ENDIF 959 960 ELSE 961 962 ! Point data 963 964 ! vertically interpolate all 4 corners 965 ista = prodatqc%npvsta(jobs,2) 966 iend = prodatqc%npvend(jobs,2) 967 inum_obs = iend - ista + 1 968 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 969 970 DO iin=1,2 971 DO ijn=1,2 972 973 974 IF ( k1dint == 1 ) THEN 975 CALL obs_int_z1d_spl( kpk, & 976 & zints(iin,ijn,:,iobs),& 977 & zobs2k, zgdept(iin,ijn,:,iobs), & 978 & zmask(iin,ijn,:,iobs)) 979 980 ENDIF 981 982 CALL obs_level_search(kpk, & 983 & zgdept(iin,ijn,:,iobs),& 984 & inum_obs, prodatqc%var(2)%vdep(ista:iend), & 985 & iv_indic) 986 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 987 & prodatqc%var(2)%vdep(ista:iend), & 988 & zints(iin,ijn,:,iobs), & 989 & zobs2k,interp_corner(iin,ijn,:), & 990 & zgdept(iin,ijn,:,iobs), & 991 & zmask(iin,ijn,:,iobs) ) 992 993 ENDDO 994 ENDDO 995 996 ENDIF 997 998 !------------------------------------------------------------- 999 ! Compute the horizontal interpolation for every profile level 1000 !------------------------------------------------------------- 1001 1002 DO ikn=1,inum_obs 1003 iend=ista+ikn-1 1004 1005 l_zweig(:,:,1) = 0._wp 1006 1007 ! This code forces the horizontal weights to be 1008 ! zero IF the observation is below the bottom of the 1009 ! corners of the interpolation nodes, Or if it is in 1010 ! the mask. This is important for observations are near 1011 ! steep bathymetry 1012 DO iin=1,2 1013 DO ijn=1,2 1014 1015 depth_loop2: DO ik=kpk,2,-1 1016 IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN 1017 1018 l_zweig(iin,ijn,1) = & 1019 & zweig(iin,ijn,1) * & 1020 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 1021 & - prodatqc%var(2)%vdep(iend)),0._wp) 1022 1023 EXIT depth_loop2 1024 ENDIF 1025 ENDDO depth_loop2 1026 1027 ENDDO 1028 ENDDO 1029 1030 CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), & 1031 & prodatqc%var(2)%vmod(iend:iend) ) 1032 1033 ENDDO 1034 1035 1036 DEALLOCATE(interp_corner,iv_indic) 1037 1038 ENDIF 1039 1040 ENDIF 1041 1042 END DO 1043 1044 ! Deallocate the data for interpolation 1045 DEALLOCATE( & 1046 & igrdi, & 1047 & igrdj, & 1048 & zglam, & 1049 & zgphi, & 1050 & zmask, & 1051 & zintt, & 1052 & zints, & 1053 & zgdept,& 1054 & zgdepw & 1055 & ) 1056 ! At the end of the day also get interpolated means 1057 IF ( idayend == 0 ) THEN 1058 DEALLOCATE( & 1059 & zinmt, & 1060 & zinms & 1061 & ) 1062 ENDIF 1063 1064 prodatqc%nprofup = prodatqc%nprofup + ipro 1065 1066 END SUBROUTINE obs_pro_sco_opt 1067 1068 SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, & 1069 & kit000, kdaystp, psurf, psurfmask, & 1070 & k2dint, ldnightav ) 459 ENDIF 460 461 ENDDO 462 463 ! Deallocate the data for interpolation 464 DEALLOCATE( & 465 & igrdi, & 466 & igrdj, & 467 & zglam, & 468 & zgphi, & 469 & zmask, & 470 & zint, & 471 & zgdept, & 472 & zgdepw & 473 & ) 474 475 ! At the end of the day also get interpolated means 476 IF ( ld_dailyav .AND. idayend == 0 ) THEN 477 DEALLOCATE( zinm ) 478 ENDIF 479 480 IF ( kvar == prodatqc%nvar ) THEN 481 prodatqc%nprofup = prodatqc%nprofup + ipro 482 ENDIF 483 484 END SUBROUTINE obs_prof_opt 485 486 SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, & 487 & kit000, kdaystp, psurf, psurfmask, & 488 & k2dint, ldnightav, plamscl, pphiscl, & 489 & lindegrees ) 1071 490 1072 491 !!----------------------------------------------------------------------- … … 1090 509 !! - polynomial (quadrilateral grid) (k2dint = 4) 1091 510 !! 511 !! Two horizontal averaging schemes are also available: 512 !! - weighted radial footprint (k2dint = 5) 513 !! - weighted rectangular footprint (k2dint = 6) 514 !! 1092 515 !! 1093 516 !! ** Action : … … 1096 519 !! ! 07-03 (A. Weaver) 1097 520 !! ! 15-02 (M. Martin) Combined routine for surface types 521 !! ! 17-03 (M. Martin) Added horizontal averaging options 1098 522 !!----------------------------------------------------------------------- 1099 523 … … 1117 541 & psurfmask ! Land-sea mask 1118 542 LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data 543 REAL(KIND=wp), INTENT(IN) :: & 544 & plamscl, & ! Diameter in metres of obs footprint in E/W, N/S directions 545 & pphiscl ! This is the full width (rather than half-width) 546 LOGICAL, INTENT(IN) :: & 547 & lindegrees ! T=> plamscl and pphiscl are specified in degrees, F=> in metres 1119 548 1120 549 !! * Local declarations … … 1125 554 INTEGER :: isurf 1126 555 INTEGER :: iobs 556 INTEGER :: imaxifp, imaxjfp 557 INTEGER :: imodi, imodj 1127 558 INTEGER :: idayend 1128 559 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 1129 & igrdi, & 1130 & igrdj 560 & igrdi, & 561 & igrdj, & 562 & igrdip1, & 563 & igrdjp1 1131 564 INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 1132 565 & icount_night, & … … 1136 569 REAL(wp), DIMENSION(1) :: zext, zobsmask 1137 570 REAL(wp) :: zdaystp 1138 REAL(wp), DIMENSION(2,2,1) :: &1139 & zweig1140 571 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 572 & zweig, & 1141 573 & zmask, & 1142 574 & zsurf, & 1143 575 & zsurfm, & 576 & zsurftmp, & 1144 577 & zglam, & 1145 & zgphi 578 & zgphi, & 579 & zglamf, & 580 & zgphif 581 1146 582 REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 1147 583 & zintmp, & … … 1155 591 inrc = kt - kit000 + 2 1156 592 isurf = surfdataqc%nsstp(inrc) 593 594 ! Work out the maximum footprint size for the 595 ! interpolation/averaging in model grid-points - has to be even. 596 597 CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp ) 598 1157 599 1158 600 IF ( ldnightav ) THEN … … 1221 663 1222 664 ALLOCATE( & 1223 & igrdi(2,2,isurf), & 1224 & igrdj(2,2,isurf), & 1225 & zglam(2,2,isurf), & 1226 & zgphi(2,2,isurf), & 1227 & zmask(2,2,isurf), & 1228 & zsurf(2,2,isurf) & 665 & zweig(imaxifp,imaxjfp,1), & 666 & igrdi(imaxifp,imaxjfp,isurf), & 667 & igrdj(imaxifp,imaxjfp,isurf), & 668 & zglam(imaxifp,imaxjfp,isurf), & 669 & zgphi(imaxifp,imaxjfp,isurf), & 670 & zmask(imaxifp,imaxjfp,isurf), & 671 & zsurf(imaxifp,imaxjfp,isurf), & 672 & zsurftmp(imaxifp,imaxjfp,isurf), & 673 & zglamf(imaxifp+1,imaxjfp+1,isurf), & 674 & zgphif(imaxifp+1,imaxjfp+1,isurf), & 675 & igrdip1(imaxifp+1,imaxjfp+1,isurf), & 676 & igrdjp1(imaxifp+1,imaxjfp+1,isurf) & 1229 677 & ) 1230 678 1231 679 DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 1232 680 iobs = jobs - surfdataqc%nsurfup 1233 igrdi(1,1,iobs) = surfdataqc%mi(jobs)-1 1234 igrdj(1,1,iobs) = surfdataqc%mj(jobs)-1 1235 igrdi(1,2,iobs) = surfdataqc%mi(jobs)-1 1236 igrdj(1,2,iobs) = surfdataqc%mj(jobs) 1237 igrdi(2,1,iobs) = surfdataqc%mi(jobs) 1238 igrdj(2,1,iobs) = surfdataqc%mj(jobs)-1 1239 igrdi(2,2,iobs) = surfdataqc%mi(jobs) 1240 igrdj(2,2,iobs) = surfdataqc%mj(jobs) 681 DO ji = 0, imaxifp 682 imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1 683 684 !Deal with wrap around in longitude 685 IF ( imodi < 1 ) imodi = imodi + jpiglo 686 IF ( imodi > jpiglo ) imodi = imodi - jpiglo 687 688 DO jj = 0, imaxjfp 689 imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1 690 !If model values are out of the domain to the north/south then 691 !set them to be the edge of the domain 692 IF ( imodj < 1 ) imodj = 1 693 IF ( imodj > jpjglo ) imodj = jpjglo 694 695 igrdip1(ji+1,jj+1,iobs) = imodi 696 igrdjp1(ji+1,jj+1,iobs) = imodj 697 698 IF ( ji >= 1 .AND. jj >= 1 ) THEN 699 igrdi(ji,jj,iobs) = imodi 700 igrdj(ji,jj,iobs) = imodj 701 ENDIF 702 703 END DO 704 END DO 1241 705 END DO 1242 706 1243 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &707 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 1244 708 & igrdi, igrdj, glamt, zglam ) 1245 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &709 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 1246 710 & igrdi, igrdj, gphit, zgphi ) 1247 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &711 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 1248 712 & igrdi, igrdj, psurfmask, zmask ) 1249 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, &713 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 1250 714 & igrdi, igrdj, psurf, zsurf ) 715 CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 716 & igrdip1, igrdjp1, glamf, zglamf ) 717 CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 718 & igrdip1, igrdjp1, gphif, zgphif ) 1251 719 1252 720 ! At the end of the day get interpolated means 1253 IF (ldnightav ) THEN 1254 IF ( idayend == 0 ) THEN 1255 1256 ALLOCATE( & 1257 & zsurfm(2,2,isurf) & 1258 & ) 1259 1260 CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, igrdi, igrdj, & 1261 & surfdataqc%vdmean(:,:), zsurfm ) 1262 1263 ENDIF 721 IF ( idayend == 0 .AND. ldnightav ) THEN 722 723 ALLOCATE( & 724 & zsurfm(imaxifp,imaxjfp,isurf) & 725 & ) 726 727 CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, & 728 & surfdataqc%vdmean(:,:), zsurfm ) 729 1264 730 ENDIF 1265 731 … … 1290 756 zphi = surfdataqc%rphi(jobs) 1291 757 1292 ! Get weights to interpolate the model value to the observation point1293 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, &1294 & zglam(:,:,iobs), zgphi(:,:,iobs), &1295 & zmask(:,:,iobs), zweig, zobsmask )1296 1297 ! Interpolate the model field to the observation point1298 758 IF ( ldnightav .AND. idayend == 0 ) THEN 1299 759 ! Night-time averaged data 1300 CALL obs_int_h2d( 1, 1, zweig, zsurfm(:,:,iobs), zext)760 zsurftmp(:,:,iobs) = zsurfm(:,:,iobs) 1301 761 ELSE 1302 CALL obs_int_h2d( 1, 1, zweig, zsurf(:,:,iobs), zext ) 762 zsurftmp(:,:,iobs) = zsurf(:,:,iobs) 763 ENDIF 764 765 IF ( k2dint <= 4 ) THEN 766 767 ! Get weights to interpolate the model value to the observation point 768 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 769 & zglam(:,:,iobs), zgphi(:,:,iobs), & 770 & zmask(:,:,iobs), zweig, zobsmask ) 771 772 ! Interpolate the model value to the observation point 773 CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 774 775 ELSE 776 777 ! Get weights to average the model SLA to the observation footprint 778 CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam, zphi, & 779 & zglam(:,:,iobs), zgphi(:,:,iobs), & 780 & zglamf(:,:,iobs), zgphif(:,:,iobs), & 781 & zmask(:,:,iobs), plamscl, pphiscl, & 782 & lindegrees, zweig, zobsmask ) 783 784 ! Average the model SST to the observation footprint 785 CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 786 & zweig, zsurftmp(:,:,iobs), zext ) 787 1303 788 ENDIF 1304 789 … … 1310 795 surfdataqc%rmod(jobs,1) = zext(1) 1311 796 ENDIF 797 798 IF ( zext(1) == obfillflt ) THEN 799 ! If the observation value is a fill value, set QC flag to bad 800 surfdataqc%nqc(jobs) = 4 801 ENDIF 1312 802 1313 803 END DO … … 1315 805 ! Deallocate the data for interpolation 1316 806 DEALLOCATE( & 807 & zweig, & 1317 808 & igrdi, & 1318 809 & igrdj, & … … 1320 811 & zgphi, & 1321 812 & zmask, & 1322 & zsurf & 813 & zsurf, & 814 & zsurftmp, & 815 & zglamf, & 816 & zgphif, & 817 & igrdip1,& 818 & igrdjp1 & 1323 819 & ) 1324 820 1325 821 ! At the end of the day also deallocate night-time mean array 1326 IF ( ldnightav ) THEN 1327 IF ( idayend == 0 ) THEN 1328 DEALLOCATE( & 1329 & zsurfm & 1330 & ) 1331 ENDIF 822 IF ( idayend == 0 .AND. ldnightav ) THEN 823 DEALLOCATE( & 824 & zsurfm & 825 & ) 1332 826 ENDIF 1333 827 -
NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90
r11350 r11361 13 13 !! obs_sor : Sort the observation arrays 14 14 !!--------------------------------------------------------------------- 15 USE par_kind, ONLY : wp ! Precision variables 15 !! * Modules used 16 USE par_kind, ONLY : & ! Precision variables 17 & wp 16 18 USE in_out_manager ! I/O manager 17 19 USE obs_profiles_def ! Definitions for storage arrays for profiles … … 22 24 USE obs_inter_sup ! Interpolation support 23 25 USE obs_oper ! Observation operators 24 USE lib_mpp, ONLY : ctl_warn, ctl_stop 26 #if defined key_bdy 27 USE bdy_oce, ONLY : & ! Boundary information 28 idx_bdy, nb_bdy 29 #endif 30 USE lib_mpp, ONLY : & 31 & ctl_warn, ctl_stop 25 32 26 33 IMPLICIT NONE 34 35 !! * Routine accessibility 27 36 PRIVATE 28 37 29 PUBLIC obs_pre_prof ! First level check and screening of profile obs 30 PUBLIC obs_pre_surf ! First level check and screening of surface obs 31 PUBLIC calc_month_len ! Calculate the number of days in the months of a year 38 PUBLIC & 39 & obs_pre_prof, & ! First level check and screening of profile obs 40 & obs_pre_surf, & ! First level check and screening of surface obs 41 & calc_month_len ! Calculate the number of days in the months of a year 32 42 33 43 !!---------------------------------------------------------------------- … … 37 47 !!---------------------------------------------------------------------- 38 48 39 40 49 CONTAINS 41 50 42 SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea ) 51 SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject, & 52 kqc_cutoff ) 43 53 !!---------------------------------------------------------------------- 44 54 !! *** ROUTINE obs_pre_sla *** … … 57 67 !! ! 2015-02 (M. Martin) Combined routine for surface types. 58 68 !!---------------------------------------------------------------------- 69 !! * Modules used 59 70 USE par_oce ! Ocean parameters 60 USE dom_oce, ONLY : glamt, gphit, tmask, nproc ! Geographical information 71 USE dom_oce, ONLY : & ! Geographical information 72 & glamt, & 73 & gphit, & 74 & tmask, & 75 & nproc 61 76 !! * Arguments 62 77 TYPE(obs_surf), INTENT(INOUT) :: surfdata ! Full set of surface data 63 TYPE(obs_surf), INTENT(INOUT) :: surfdataqc ! Subset of surface data not failing screening 64 LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land 65 ! 78 TYPE(obs_surf), INTENT(INOUT) :: surfdataqc ! Subset of surface data not failing screening 79 LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land 80 LOGICAL, INTENT(IN) :: ld_bound_reject ! Switch for rejecting obs near the boundary 81 INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff ! cut off for QC value 82 !! * Local declarations 83 INTEGER :: iqc_cutoff = 255 ! cut off for QC value 66 84 INTEGER :: iyea0 ! Initial date 67 85 INTEGER :: imon0 ! - (year, month, day, hour, minute) … … 76 94 INTEGER :: inlasobs ! - close to land 77 95 INTEGER :: igrdobs ! - fail the grid search 96 INTEGER :: ibdysobs ! - close to open boundary 78 97 ! Global counters for observations that 79 98 INTEGER :: iotdobsmpp ! - outside time domain … … 82 101 INTEGER :: inlasobsmpp ! - close to land 83 102 INTEGER :: igrdobsmpp ! - fail the grid search 84 LOGICAL, DIMENSION(:), ALLOCATABLE :: llvalid ! SLA data selection 103 INTEGER :: ibdysobsmpp ! - close to open boundary 104 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 105 & llvalid ! SLA data selection 85 106 INTEGER :: jobs ! Obs. loop variable 86 107 INTEGER :: jstp ! Time loop variable 87 108 INTEGER :: inrc ! Time index variable 88 !!---------------------------------------------------------------------- 89 90 IF(lwp) WRITE(numout,*) 'obs_pre_surf : Preparing the surface observations...' 91 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 109 110 IF(lwp) WRITE(numout,*)'obs_pre_surf : Preparing the surface observations...' 111 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 92 112 93 113 ! Initial date initialization (year, month, day, hour, minute) … … 95 115 imon0 = ( ndate0 - iyea0 * 10000 ) / 100 96 116 iday0 = ndate0 - iyea0 * 10000 - imon0 * 100 97 ihou0 = nn_time0 / 10098 imin0 = ( nn_time0 - ihou0 * 100 )117 ihou0 = 0 118 imin0 = 0 99 119 100 120 icycle = no ! Assimilation cycle … … 107 127 ilansobs = 0 108 128 inlasobs = 0 129 ibdysobs = 0 130 131 ! Set QC cutoff to optional value if provided 132 IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff 109 133 110 134 ! ----------------------------------------------------------------------- … … 140 164 & tmask(:,:,1), surfdata%nqc, & 141 165 & iosdsobs, ilansobs, & 142 & inlasobs, ld_nea ) 166 & inlasobs, ld_nea, & 167 & ibdysobs, ld_bound_reject, & 168 & iqc_cutoff ) 143 169 144 170 CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 145 171 CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 146 172 CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 173 CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp ) 147 174 148 175 ! ----------------------------------------------------------------------- … … 155 182 ALLOCATE( llvalid(surfdata%nsurf) ) 156 183 157 ! We want all data which has qc flags <= 10158 159 llvalid(:) = ( surfdata%nqc(:) <= 10)184 ! We want all data which has qc flags <= iqc_cutoff 185 186 llvalid(:) = ( surfdata%nqc(:) <= iqc_cutoff ) 160 187 161 188 ! The actual copying … … 190 217 & inlasobsmpp 191 218 ENDIF 219 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near open boundary (removed) = ', & 220 & ibdysobsmpp 192 221 WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data accepted = ', & 193 222 & surfdataqc%nsurfmpp … … 222 251 223 252 224 SUBROUTINE obs_pre_prof( profdata, prodatqc, ld_var 1, ld_var2, &253 SUBROUTINE obs_pre_prof( profdata, prodatqc, ld_var, & 225 254 & kpi, kpj, kpk, & 226 & zmask 1, pglam1, pgphi1, zmask2, pglam2, pgphi2, &227 & ld_nea, kdailyavtypes)255 & zmask, pglam, pgphi, & 256 & ld_nea, ld_bound_reject, kdailyavtypes, kqc_cutoff ) 228 257 229 258 !!---------------------------------------------------------------------- … … 241 270 !! 242 271 !!---------------------------------------------------------------------- 243 USE par_oce ! Ocean parameters 244 USE dom_oce, ONLY : gdept_1d, nproc ! Geographical information 272 !! * Modules used 273 USE par_oce ! Ocean parameters 274 USE dom_oce, ONLY : & ! Geographical information 275 & gdept_1d, & 276 & nproc 245 277 246 278 !! * Arguments 247 279 TYPE(obs_prof), INTENT(INOUT) :: profdata ! Full set of profile data 248 280 TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening 249 LOGICAL, INTENT(IN) :: ld_var1 ! Observed variables switches250 LOGICAL, INTENT(IN) :: ld_var2281 LOGICAL, DIMENSION(profdata%nvar), INTENT(IN) :: & 282 & ld_var ! Observed variables switches 251 283 LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land 284 LOGICAL, INTENT(IN) :: ld_bound_reject ! Switch for rejecting observations near the boundary 252 285 INTEGER, INTENT(IN) :: kpi, kpj, kpk ! Local domain sizes 253 286 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 254 287 & kdailyavtypes ! Types for daily averages 255 REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 256 & zmask1, & 257 & zmask2 258 REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 259 & pglam1, & 260 & pglam2, & 261 & pgphi1, & 262 & pgphi2 288 REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,kpk,profdata%nvar) :: & 289 & zmask 290 REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,profdata%nvar) :: & 291 & pglam, & 292 & pgphi 293 INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff ! cut off for QC value 263 294 264 295 !! * Local declarations 296 INTEGER :: iqc_cutoff = 255 ! cut off for QC value 265 297 INTEGER :: iyea0 ! Initial date 266 298 INTEGER :: imon0 ! - (year, month, day, hour, minute) … … 269 301 INTEGER :: imin0 270 302 INTEGER :: icycle ! Current assimilation cycle 271 ! Counters for observations that are 272 INTEGER :: iotdobs ! - outside time domain 273 INTEGER :: iosdv1obs ! - outside space domain (variable 1) 274 INTEGER :: iosdv2obs ! - outside space domain (variable 2) 275 INTEGER :: ilanv1obs ! - within a model land cell (variable 1) 276 INTEGER :: ilanv2obs ! - within a model land cell (variable 2) 277 INTEGER :: inlav1obs ! - close to land (variable 1) 278 INTEGER :: inlav2obs ! - close to land (variable 2) 279 INTEGER :: igrdobs ! - fail the grid search 280 INTEGER :: iuvchku ! - reject u if v rejected and vice versa 281 INTEGER :: iuvchkv ! 282 ! Global counters for observations that are 283 INTEGER :: iotdobsmpp ! - outside time domain 284 INTEGER :: iosdv1obsmpp ! - outside space domain (variable 1) 285 INTEGER :: iosdv2obsmpp ! - outside space domain (variable 2) 286 INTEGER :: ilanv1obsmpp ! - within a model land cell (variable 1) 287 INTEGER :: ilanv2obsmpp ! - within a model land cell (variable 2) 288 INTEGER :: inlav1obsmpp ! - close to land (variable 1) 289 INTEGER :: inlav2obsmpp ! - close to land (variable 2) 290 INTEGER :: igrdobsmpp ! - fail the grid search 291 INTEGER :: iuvchkumpp ! - reject var1 if var2 rejected and vice versa 292 INTEGER :: iuvchkvmpp ! 303 ! Counters for observations that are 304 INTEGER :: iotdobs ! - outside time domain 305 INTEGER, DIMENSION(profdata%nvar) :: iosdvobs ! - outside space domain 306 INTEGER, DIMENSION(profdata%nvar) :: ilanvobs ! - within a model land cell 307 INTEGER, DIMENSION(profdata%nvar) :: inlavobs ! - close to land 308 INTEGER, DIMENSION(profdata%nvar) :: ibdyvobs ! - boundary 309 INTEGER :: igrdobs ! - fail the grid search 310 INTEGER :: iuvchku ! - reject UVEL if VVEL rejected 311 INTEGER :: iuvchkv ! - reject VVEL if UVEL rejected 312 ! Global counters for observations that are 313 INTEGER :: iotdobsmpp ! - outside time domain 314 INTEGER, DIMENSION(profdata%nvar) :: iosdvobsmpp ! - outside space domain 315 INTEGER, DIMENSION(profdata%nvar) :: ilanvobsmpp ! - within a model land cell 316 INTEGER, DIMENSION(profdata%nvar) :: inlavobsmpp ! - close to land 317 INTEGER, DIMENSION(profdata%nvar) :: ibdyvobsmpp ! - boundary 318 INTEGER :: igrdobsmpp ! - fail the grid search 319 INTEGER :: iuvchkumpp ! - reject UVEL if VVEL rejected 320 INTEGER :: iuvchkvmpp ! - reject VVEL if UVEL rejected 293 321 TYPE(obs_prof_valid) :: llvalid ! Profile selection 294 322 TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: & 295 & llvvalid ! var 1,var2selection323 & llvvalid ! vars selection 296 324 INTEGER :: jvar ! Variable loop variable 297 325 INTEGER :: jobs ! Obs. loop variable 298 326 INTEGER :: jstp ! Time loop variable 299 327 INTEGER :: inrc ! Time index variable 300 !!---------------------------------------------------------------------- 328 CHARACTER(LEN=256) :: cout1 ! Diagnostic output line 329 CHARACTER(LEN=256) :: cout2 ! Diagnostic output line 301 330 302 331 IF(lwp) WRITE(numout,*)'obs_pre_prof: Preparing the profile data...' … … 307 336 imon0 = ( ndate0 - iyea0 * 10000 ) / 100 308 337 iday0 = ndate0 - iyea0 * 10000 - imon0 * 100 309 ihou0 = nn_time0 / 100310 imin0 = ( nn_time0 - ihou0 * 100 )338 ihou0 = 0 339 imin0 = 0 311 340 312 341 icycle = no ! Assimilation cycle 313 342 314 ! Diagnotics counters for various failures. 315 316 iotdobs = 0 317 igrdobs = 0 318 iosdv1obs = 0 319 iosdv2obs = 0 320 ilanv1obs = 0 321 ilanv2obs = 0 322 inlav1obs = 0 323 inlav2obs = 0 324 iuvchku = 0 325 iuvchkv = 0 343 ! Diagnostics counters for various failures. 344 345 iotdobs = 0 346 igrdobs = 0 347 iosdvobs(:) = 0 348 ilanvobs(:) = 0 349 inlavobs(:) = 0 350 ibdyvobs(:) = 0 351 iuvchku = 0 352 iuvchkv = 0 353 354 355 ! Set QC cutoff to optional value if provided 356 IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff 326 357 327 358 ! ----------------------------------------------------------------------- … … 335 366 & profdata%nday, profdata%nhou, profdata%nmin, & 336 367 & profdata%ntyp, profdata%nqc, profdata%mstp, & 337 & iotdobs, kdailyavtypes = kdailyavtypes ) 368 & iotdobs, kdailyavtypes = kdailyavtypes, & 369 & kqc_cutoff = iqc_cutoff ) 338 370 ELSE 339 371 CALL obs_coo_tim_prof( icycle, & … … 342 374 & profdata%nday, profdata%nhou, profdata%nmin, & 343 375 & profdata%ntyp, profdata%nqc, profdata%mstp, & 344 & iotdobs )376 & iotdobs, kqc_cutoff = iqc_cutoff ) 345 377 ENDIF 346 378 … … 351 383 ! ----------------------------------------------------------------------- 352 384 353 CALL obs_coo_grd( profdata%nprof, profdata%mi(:,1), profdata%mj(:,1), &354 & profdata%nqc, igrdobs )355 CALL obs_coo_grd( profdata%nprof, profdata%mi(:,2), profdata%mj(:,2), &356 & profdata%nqc, igrdobs )385 DO jvar = 1, profdata%nvar 386 CALL obs_coo_grd( profdata%nprof, profdata%mi(:,jvar), profdata%mj(:,jvar), & 387 & profdata%nqc, igrdobs ) 388 END DO 357 389 358 390 CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) 359 391 360 392 ! ----------------------------------------------------------------------- 361 ! Reject all observations for profiles with nqc > 10362 ! ----------------------------------------------------------------------- 363 364 CALL obs_pro_rej( profdata )393 ! Reject all observations for profiles with nqc > iqc_cutoff 394 ! ----------------------------------------------------------------------- 395 396 CALL obs_pro_rej( profdata, kqc_cutoff = iqc_cutoff ) 365 397 366 398 ! ----------------------------------------------------------------------- … … 369 401 ! ----------------------------------------------------------------------- 370 402 371 ! Variable 1 372 CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(1), & 373 & profdata%npvsta(:,1), profdata%npvend(:,1), & 374 & jpi, jpj, & 375 & jpk, & 376 & profdata%mi, profdata%mj, & 377 & profdata%var(1)%mvk, & 378 & profdata%rlam, profdata%rphi, & 379 & profdata%var(1)%vdep, & 380 & pglam1, pgphi1, & 381 & gdept_1d, zmask1, & 382 & profdata%nqc, profdata%var(1)%nvqc, & 383 & iosdv1obs, ilanv1obs, & 384 & inlav1obs, ld_nea ) 385 386 CALL obs_mpp_sum_integer( iosdv1obs, iosdv1obsmpp ) 387 CALL obs_mpp_sum_integer( ilanv1obs, ilanv1obsmpp ) 388 CALL obs_mpp_sum_integer( inlav1obs, inlav1obsmpp ) 389 390 ! Variable 2 391 CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(2), & 392 & profdata%npvsta(:,2), profdata%npvend(:,2), & 393 & jpi, jpj, & 394 & jpk, & 395 & profdata%mi, profdata%mj, & 396 & profdata%var(2)%mvk, & 397 & profdata%rlam, profdata%rphi, & 398 & profdata%var(2)%vdep, & 399 & pglam2, pgphi2, & 400 & gdept_1d, zmask2, & 401 & profdata%nqc, profdata%var(2)%nvqc, & 402 & iosdv2obs, ilanv2obs, & 403 & inlav2obs, ld_nea ) 404 405 CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp ) 406 CALL obs_mpp_sum_integer( ilanv2obs, ilanv2obsmpp ) 407 CALL obs_mpp_sum_integer( inlav2obs, inlav2obsmpp ) 403 DO jvar = 1, profdata%nvar 404 CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(jvar), & 405 & profdata%npvsta(:,jvar), profdata%npvend(:,jvar), & 406 & jpi, jpj, & 407 & jpk, & 408 & profdata%mi, profdata%mj, & 409 & profdata%var(jvar)%mvk, & 410 & profdata%rlam, profdata%rphi, & 411 & profdata%var(jvar)%vdep, & 412 & pglam(:,:,jvar), pgphi(:,:,jvar), & 413 & gdept_1d, zmask(:,:,:,jvar), & 414 & profdata%nqc, profdata%var(jvar)%nvqc, & 415 & iosdvobs(jvar), ilanvobs(jvar), & 416 & inlavobs(jvar), ld_nea, & 417 & ibdyvobs(jvar), ld_bound_reject, & 418 & iqc_cutoff ) 419 420 CALL obs_mpp_sum_integer( iosdvobs(jvar), iosdvobsmpp(jvar) ) 421 CALL obs_mpp_sum_integer( ilanvobs(jvar), ilanvobsmpp(jvar) ) 422 CALL obs_mpp_sum_integer( inlavobs(jvar), inlavobsmpp(jvar) ) 423 CALL obs_mpp_sum_integer( ibdyvobs(jvar), ibdyvobsmpp(jvar) ) 424 END DO 408 425 409 426 ! ----------------------------------------------------------------------- … … 412 429 413 430 IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 414 CALL obs_uv_rej( profdata, iuvchku, iuvchkv )431 CALL obs_uv_rej( profdata, iuvchku, iuvchkv, iqc_cutoff ) 415 432 CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp ) 416 433 CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp ) … … 429 446 END DO 430 447 431 ! We want all data which has qc flags = 0432 433 llvalid%luse(:) = ( profdata%nqc(:) <= 10)448 ! We want all data which has qc flags <= iqc_cutoff 449 450 llvalid%luse(:) = ( profdata%nqc(:) <= iqc_cutoff ) 434 451 DO jvar = 1,profdata%nvar 435 llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10)452 llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= iqc_cutoff ) 436 453 END DO 437 454 … … 456 473 457 474 WRITE(numout,*) 458 WRITE(numout,*) ' Profiles outside time domain = ', &475 WRITE(numout,*) ' Profiles outside time domain = ', & 459 476 & iotdobsmpp 460 WRITE(numout,*) ' Remaining profiles that failed grid search = ', &477 WRITE(numout,*) ' Remaining profiles that failed grid search = ', & 461 478 & igrdobsmpp 462 WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data outside space domain = ', & 463 & iosdv1obsmpp 464 WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data at land points = ', & 465 & ilanv1obsmpp 466 IF (ld_nea) THEN 467 WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near land points (removed) = ',& 468 & inlav1obsmpp 469 ELSE 470 WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near land points (kept) = ',& 471 & inlav1obsmpp 472 ENDIF 473 IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 474 WRITE(numout,*) ' U observation rejected since V rejected = ', & 475 & iuvchku 476 ENDIF 477 WRITE(numout,*) ' '//prodatqc%cvars(1)//' data accepted = ', & 478 & prodatqc%nvprotmpp(1) 479 WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data outside space domain = ', & 480 & iosdv2obsmpp 481 WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data at land points = ', & 482 & ilanv2obsmpp 483 IF (ld_nea) THEN 484 WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (removed) = ',& 485 & inlav2obsmpp 486 ELSE 487 WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (kept) = ',& 488 & inlav2obsmpp 489 ENDIF 490 IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 491 WRITE(numout,*) ' V observation rejected since U rejected = ', & 492 & iuvchkv 493 ENDIF 494 WRITE(numout,*) ' '//prodatqc%cvars(2)//' data accepted = ', & 495 & prodatqc%nvprotmpp(2) 479 DO jvar = 1, profdata%nvar 480 WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data outside space domain = ', & 481 & iosdvobsmpp(jvar) 482 WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data at land points = ', & 483 & ilanvobsmpp(jvar) 484 IF (ld_nea) THEN 485 WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near land points (removed) = ',& 486 & inlavobsmpp(jvar) 487 ELSE 488 WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near land points (kept) = ',& 489 & inlavobsmpp(jvar) 490 ENDIF 491 IF ( TRIM(profdata%cvars(jvar)) == 'UVEL' ) THEN 492 WRITE(numout,*) ' U observation rejected since V rejected = ', & 493 & iuvchku 494 ELSE IF ( TRIM(profdata%cvars(jvar)) == 'VVEL' ) THEN 495 WRITE(numout,*) ' V observation rejected since U rejected = ', & 496 & iuvchkv 497 ENDIF 498 WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near open boundary (removed) = ',& 499 & ibdyvobsmpp(jvar) 500 WRITE(numout,*) ' '//prodatqc%cvars(jvar)//' data accepted = ', & 501 & prodatqc%nvprotmpp(jvar) 502 END DO 496 503 497 504 WRITE(numout,*) 498 505 WRITE(numout,*) ' Number of observations per time step :' 499 506 WRITE(numout,*) 500 WRITE(numout,'(10X,A,5X,A,5X,A,A)')'Time step','Profiles', & 501 & ' '//prodatqc%cvars(1)//' ', & 502 & ' '//prodatqc%cvars(2)//' ' 503 WRITE(numout,998) 507 WRITE(cout1,'(10X,A9,5X,A8)') 'Time step', 'Profiles' 508 WRITE(cout2,'(10X,A9,5X,A8)') '---------', '--------' 509 DO jvar = 1, prodatqc%nvar 510 WRITE(cout1,'(A,5X,A11)') TRIM(cout1), TRIM(prodatqc%cvars(jvar)) 511 WRITE(cout2,'(A,5X,A11)') TRIM(cout2), '-----------' 512 END DO 513 WRITE(numout,*) cout1 514 WRITE(numout,*) cout2 504 515 ENDIF 505 516 … … 528 539 DO jstp = nit000 - 1, nitend 529 540 inrc = jstp - nit000 + 2 530 WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), & 531 & prodatqc%nvstpmpp(inrc,1), & 532 & prodatqc%nvstpmpp(inrc,2) 541 WRITE(cout1,'(10X,I9,5X,I8)') jstp, prodatqc%npstpmpp(inrc) 542 DO jvar = 1, prodatqc%nvar 543 WRITE(cout1,'(A,5X,I11)') TRIM(cout1), prodatqc%nvstpmpp(inrc,jvar) 544 END DO 545 WRITE(numout,*) cout1 533 546 END DO 534 547 ENDIF 535 536 998 FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'----------------')537 999 FORMAT(10X,I9,5X,I8,5X,I11,5X,I8)538 548 539 549 END SUBROUTINE obs_pre_prof … … 644 654 & .AND. ( kobsmin(jobs) <= kmin0 ) ) ) THEN 645 655 kobsstp(jobs) = -1 646 kobsqc(jobs) = kobsqc(jobs) + 11656 kobsqc(jobs) = IBSET(kobsqc(jobs),13) 647 657 kotdobs = kotdobs + 1 648 658 CYCLE … … 695 705 IF ( ( kobsstp(jobs) < ( nit000 - 1 ) ) & 696 706 & .OR.( kobsstp(jobs) > nitend ) ) THEN 697 kobsqc(jobs) = kobsqc(jobs) + 12707 kobsqc(jobs) = IBSET(kobsqc(jobs),13) 698 708 kotdobs = kotdobs + 1 699 709 CYCLE … … 739 749 & kobsno, & 740 750 & kobsyea, kobsmon, kobsday, kobshou, kobsmin, & 741 & ktyp, kobsqc, kobsstp, kotdobs, kdailyavtypes ) 751 & ktyp, kobsqc, kobsstp, kotdobs, kdailyavtypes, & 752 & kqc_cutoff ) 742 753 !!---------------------------------------------------------------------- 743 754 !! *** ROUTINE obs_coo_tim *** … … 783 794 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 784 795 & kdailyavtypes ! Types for daily averages 796 INTEGER, OPTIONAL, INTENT(IN) :: kqc_cutoff ! QC cutoff value 797 785 798 !! * Local declarations 786 799 INTEGER :: jobs 800 INTEGER :: iqc_cutoff=255 787 801 788 802 !----------------------------------------------------------------------- … … 803 817 DO jobs = 1, kobsno 804 818 805 IF ( kobsqc(jobs) <= 10) THEN819 IF ( kobsqc(jobs) <= iqc_cutoff ) THEN 806 820 807 821 IF ( ( kobsstp(jobs) == (nit000 - 1) ).AND.& 808 822 & ( ANY (kdailyavtypes(:) == ktyp(jobs)) ) ) THEN 809 kobsqc(jobs) = kobsqc(jobs) + 14823 kobsqc(jobs) = IBSET(kobsqc(jobs),13) 810 824 kotdobs = kotdobs + 1 811 825 CYCLE … … 850 864 DO jobs = 1, kobsno 851 865 IF ( ( kobsi(jobs) <= 0 ) .AND. ( kobsj(jobs) <= 0 ) ) THEN 852 kobsqc(jobs) = kobsqc(jobs) + 18866 kobsqc(jobs) = IBSET(kobsqc(jobs),12) 853 867 kgrdobs = kgrdobs + 1 854 868 ENDIF … … 861 875 & plam, pphi, pmask, & 862 876 & kobsqc, kosdobs, klanobs, & 863 & knlaobs,ld_nea ) 877 & knlaobs,ld_nea, & 878 & kbdyobs,ld_bound_reject, & 879 & kqc_cutoff ) 864 880 !!---------------------------------------------------------------------- 865 881 !! *** ROUTINE obs_coo_spc_2d *** … … 894 910 INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: & 895 911 & kobsqc ! Observation quality control 896 INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain 897 INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell 898 INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land 899 LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land 912 INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain 913 INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell 914 INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land 915 INTEGER, INTENT(INOUT) :: kbdyobs ! Observations near boundary 916 LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land 917 LOGICAL, INTENT(IN) :: ld_bound_reject ! Flag observations near open boundary 918 INTEGER, INTENT(IN) :: kqc_cutoff ! Cutoff QC value 919 900 920 !! * Local declarations 901 921 REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 902 922 & zgmsk ! Grid mask 923 #if defined key_bdy 924 REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 925 & zbmsk ! Boundary mask 926 REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 927 #endif 903 928 REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 904 929 & zglam, & ! Model longitude at grid points … … 917 942 ! For invalid points use 2,2 918 943 919 IF ( kobsqc(jobs) >= 10) THEN944 IF ( kobsqc(jobs) >= kqc_cutoff ) THEN 920 945 921 946 igrdi(1,1,jobs) = 1 … … 942 967 943 968 END DO 969 970 #if defined key_bdy 971 ! Create a mask grid points in boundary rim 972 IF (ld_bound_reject) THEN 973 zbdymask(:,:) = 1.0_wp 974 DO ji = 1, nb_bdy 975 DO jj = 1, idx_bdy(ji)%nblen(1) 976 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 977 ENDDO 978 ENDDO 979 980 CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk ) 981 ENDIF 982 #endif 944 983 945 984 CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pmask, zgmsk ) … … 950 989 951 990 ! Skip bad observations 952 IF ( kobsqc(jobs) >= 10) CYCLE991 IF ( kobsqc(jobs) >= kqc_cutoff ) CYCLE 953 992 954 993 ! Flag if the observation falls outside the model spatial domain … … 957 996 & .OR. ( pobsphi(jobs) < -90. ) & 958 997 & .OR. ( pobsphi(jobs) > 90. ) ) THEN 959 kobsqc(jobs) = kobsqc(jobs) + 11998 kobsqc(jobs) = IBSET(kobsqc(jobs),11) 960 999 kosdobs = kosdobs + 1 961 1000 CYCLE … … 964 1003 ! Flag if the observation falls with a model land cell 965 1004 IF ( SUM( zgmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 966 kobsqc(jobs) = kobsqc(jobs) + 121005 kobsqc(jobs) = IBSET(kobsqc(jobs),10) 967 1006 klanobs = klanobs + 1 968 1007 CYCLE … … 978 1017 IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) & 979 1018 & .AND. & 980 & ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp )&981 & ) THEN1019 & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) & 1020 & < 1.0e-6_wp ) ) THEN 982 1021 lgridobs = .TRUE. 983 1022 iig = ji … … 986 1025 END DO 987 1026 END DO 988 989 ! For observations on the grid reject them if their are at 990 ! a masked point 991 1027 992 1028 IF (lgridobs) THEN 993 1029 IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN 994 kobsqc(jobs) = kobsqc(jobs) + 121030 kobsqc(jobs) = IBSET(kobsqc(jobs),10) 995 1031 klanobs = klanobs + 1 996 1032 CYCLE 997 1033 ENDIF 998 1034 ENDIF 999 1035 1036 1000 1037 ! Flag if the observation falls is close to land 1001 1038 IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN 1002 IF (ld_nea) kobsqc(jobs) = kobsqc(jobs) + 141003 1039 knlaobs = knlaobs + 1 1004 CYCLE 1040 IF (ld_nea) THEN 1041 kobsqc(jobs) = IBSET(kobsqc(jobs),9) 1042 CYCLE 1043 ENDIF 1005 1044 ENDIF 1045 1046 #if defined key_bdy 1047 ! Flag if the observation falls close to the boundary rim 1048 IF (ld_bound_reject) THEN 1049 IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 1050 kobsqc(jobs) = IBSET(kobsqc(jobs),8) 1051 kbdyobs = kbdyobs + 1 1052 CYCLE 1053 ENDIF 1054 ! for observations on the grid... 1055 IF (lgridobs) THEN 1056 IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN 1057 kobsqc(jobs) = IBSET(kobsqc(jobs),8) 1058 kbdyobs = kbdyobs + 1 1059 CYCLE 1060 ENDIF 1061 ENDIF 1062 ENDIF 1063 #endif 1006 1064 1007 1065 END DO … … 1015 1073 & plam, pphi, pdep, pmask, & 1016 1074 & kpobsqc, kobsqc, kosdobs, & 1017 & klanobs, knlaobs, ld_nea ) 1075 & klanobs, knlaobs, ld_nea, & 1076 & kbdyobs, ld_bound_reject, & 1077 & kqc_cutoff ) 1018 1078 !!---------------------------------------------------------------------- 1019 1079 !! *** ROUTINE obs_coo_spc_3d *** … … 1040 1100 & gdepw_1d, & 1041 1101 & gdepw_0, & 1042 & gdepw_n, & 1102 & gdepw_n, & 1103 #if defined key_vvl 1043 1104 & gdept_n, & 1105 #endif 1044 1106 & ln_zco, & 1045 & ln_zps 1107 & ln_zps, & 1108 & ln_linssh 1046 1109 1047 1110 !! * Arguments … … 1077 1140 INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell 1078 1141 INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land 1142 INTEGER, INTENT(INOUT) :: kbdyobs ! Observations near boundary 1079 1143 LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land 1144 LOGICAL, INTENT(IN) :: ld_bound_reject ! Flag observations near open boundary 1145 INTEGER, INTENT(IN) :: kqc_cutoff ! Cutoff QC value 1146 1080 1147 !! * Local declarations 1081 1148 REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 1082 1149 & zgmsk ! Grid mask 1150 #if defined key_bdy 1151 REAL(KIND=wp), DIMENSION(2,2,kprofno) :: & 1152 & zbmsk ! Boundary mask 1153 REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 1154 #endif 1083 1155 REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 1084 & zgdepw 1156 & zgdepw 1085 1157 REAL(KIND=wp), DIMENSION(2,2,kprofno) :: & 1086 1158 & zglam, & ! Model longitude at grid points … … 1100 1172 ! For invalid points use 2,2 1101 1173 1102 IF ( kpobsqc(jobs) >= 10) THEN1174 IF ( kpobsqc(jobs) >= kqc_cutoff ) THEN 1103 1175 1104 1176 igrdi(1,1,jobs) = 1 … … 1125 1197 1126 1198 END DO 1199 1200 #if defined key_bdy 1201 ! Create a mask grid points in boundary rim 1202 IF (ld_bound_reject) THEN 1203 zbdymask(:,:) = 1.0_wp 1204 DO ji = 1, nb_bdy 1205 DO jj = 1, idx_bdy(ji)%nblen(1) 1206 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 1207 ENDDO 1208 ENDDO 1209 ENDIF 1210 1211 CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk ) 1212 #endif 1127 1213 1128 1214 CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk ) 1129 1215 CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam ) 1130 1216 CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi ) 1131 IF ( .NOT.( ln_zps .OR. ln_zco ) ) THEN 1132 ! Need to know the bathy depth for each observation for sco 1133 CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), & 1134 & zgdepw ) 1135 ENDIF 1217 ! Need to know the bathy depth for each observation for sco 1218 CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), & 1219 & zgdepw ) 1136 1220 1137 1221 DO jobs = 1, kprofno 1138 1222 1139 1223 ! Skip bad profiles 1140 IF ( kpobsqc(jobs) >= 10) CYCLE1224 IF ( kpobsqc(jobs) >= kqc_cutoff ) CYCLE 1141 1225 1142 1226 ! Check if this observation is on a grid point … … 1149 1233 IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) & 1150 1234 & .AND. & 1151 & ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) &1235 & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) < 1.0e-6_wp ) & 1152 1236 & ) THEN 1153 1237 lgridobs = .TRUE. … … 1158 1242 END DO 1159 1243 1160 ! Check if next to land 1161 IF ( ANY( zgmsk(1:2,1:2,1,jobs) == 0.0_wp ) ) THEN 1162 ll_next_to_land=.TRUE. 1163 ELSE 1164 ll_next_to_land=.FALSE. 1244 ! Check if next to land 1245 IF ( ANY( zgmsk(1:2,1:2,1,jobs) == 0.0_wp ) ) THEN 1246 ll_next_to_land=.TRUE. 1247 ELSE 1248 ll_next_to_land=.FALSE. 1165 1249 ENDIF 1166 1250 1167 1251 ! Reject observations 1168 1252 … … 1176 1260 & .OR. ( pobsdep(jobsp) < 0.0 ) & 1177 1261 & .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN 1178 kobsqc(jobsp) = kobsqc(jobsp) + 111262 kobsqc(jobsp) = IBSET(kobsqc(jobsp),11) 1179 1263 kosdobs = kosdobs + 1 1180 1264 CYCLE 1181 1265 ENDIF 1182 1266 1183 ! To check if an observations falls within land there are two cases: 1184 ! 1: z-coordibnates, where the check uses the mask 1185 ! 2: terrain following (eg s-coordinates), 1186 ! where we use the depth of the bottom cell to mask observations 1267 ! To check if an observations falls within land there are two cases: 1268 ! 1: z-coordibnates, where the check uses the mask 1269 ! 2: terrain following (eg s-coordinates), 1270 ! where we use the depth of the bottom cell to mask observations 1187 1271 1188 IF( ln_zps .OR. ln_zco ) THEN !(CASE 1) 1189 1190 ! Flag if the observation falls with a model land cell 1191 IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 1192 & == 0.0_wp ) THEN 1193 kobsqc(jobsp) = kobsqc(jobsp) + 12 1272 IF( (ln_linssh) .AND. ( ln_zps .OR. ln_zco ) ) THEN !(CASE 1) 1273 1274 ! Flag if the observation falls with a model land cell 1275 IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 1276 & == 0.0_wp ) THEN 1277 kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 1278 klanobs = klanobs + 1 1279 CYCLE 1280 ENDIF 1281 1282 ! Flag if the observation is close to land 1283 IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == & 1284 & 0.0_wp) THEN 1285 knlaobs = knlaobs + 1 1286 IF (ld_nea) THEN 1287 kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 1288 ENDIF 1289 ENDIF 1290 1291 ELSE ! Case 2 1292 ! Flag if the observation is deeper than the bathymetry 1293 ! Or if it is within the mask 1294 IF ( ANY( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) & 1295 & .OR. & 1296 & ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 1297 & == 0.0_wp) ) THEN 1298 kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 1299 klanobs = klanobs + 1 1300 CYCLE 1301 ENDIF 1302 1303 ! Flag if the observation is close to land 1304 IF ( ll_next_to_land ) THEN 1305 knlaobs = knlaobs + 1 1306 IF (ld_nea) THEN 1307 kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 1308 ENDIF 1309 ENDIF 1310 1311 ENDIF 1312 1313 ! For observations on the grid reject them if their are at 1314 ! a masked point 1315 1316 IF (lgridobs) THEN 1317 IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN 1318 kobsqc(jobsp) = IBSET(kobsqc(jobs),10) 1194 1319 klanobs = klanobs + 1 1195 1320 CYCLE 1196 1321 ENDIF 1197 1198 ! Flag if the observation is close to land 1199 IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == & 1200 & 0.0_wp) THEN 1201 knlaobs = knlaobs + 1 1202 IF (ld_nea) THEN 1203 kobsqc(jobsp) = kobsqc(jobsp) + 14 1204 ENDIF 1205 ENDIF 1206 1207 ELSE ! Case 2 1208 1209 ! Flag if the observation is deeper than the bathymetry 1210 ! Or if it is within the mask 1211 IF ( ALL( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) & 1212 & .OR. & 1213 & ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 1214 & == 0.0_wp) ) THEN 1215 kobsqc(jobsp) = kobsqc(jobsp) + 12 1216 klanobs = klanobs + 1 1217 CYCLE 1218 ENDIF 1219 1220 ! Flag if the observation is close to land 1221 IF ( ll_next_to_land ) THEN 1222 knlaobs = knlaobs + 1 1223 IF (ld_nea) THEN 1224 kobsqc(jobsp) = kobsqc(jobsp) + 14 1225 ENDIF 1226 ENDIF 1227 ENDIF 1228 1229 ! For observations on the grid reject them if their are at 1230 ! a masked point 1231 1232 IF (lgridobs) THEN 1233 IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN 1234 kobsqc(jobsp) = kobsqc(jobsp) + 12 1235 klanobs = klanobs + 1 1236 CYCLE 1237 ENDIF 1238 ENDIF 1239 1240 ! Flag if the observation falls is close to land 1241 IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == & 1242 & 0.0_wp) THEN 1243 IF (ld_nea) kobsqc(jobsp) = kobsqc(jobsp) + 14 1244 knlaobs = knlaobs + 1 1245 ENDIF 1246 1322 ENDIF 1323 1247 1324 ! Set observation depth equal to that of the first model depth 1248 1325 IF ( pobsdep(jobsp) <= pdep(1) ) THEN … … 1250 1327 ENDIF 1251 1328 1329 #if defined key_bdy 1330 ! Flag if the observation falls close to the boundary rim 1331 IF (ld_bound_reject) THEN 1332 IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 1333 kobsqc(jobsp) = IBSET(kobsqc(jobs),8) 1334 kbdyobs = kbdyobs + 1 1335 CYCLE 1336 ENDIF 1337 ! for observations on the grid... 1338 IF (lgridobs) THEN 1339 IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN 1340 kobsqc(jobsp) = IBSET(kobsqc(jobs),8) 1341 kbdyobs = kbdyobs + 1 1342 CYCLE 1343 ENDIF 1344 ENDIF 1345 ENDIF 1346 #endif 1347 1252 1348 END DO 1253 1349 END DO … … 1255 1351 END SUBROUTINE obs_coo_spc_3d 1256 1352 1257 SUBROUTINE obs_pro_rej( profdata )1353 SUBROUTINE obs_pro_rej( profdata, kqc_cutoff ) 1258 1354 !!---------------------------------------------------------------------- 1259 1355 !! *** ROUTINE obs_pro_rej *** … … 1273 1369 !! * Arguments 1274 1370 TYPE(obs_prof), INTENT(INOUT) :: profdata ! Profile data 1371 INTEGER, INTENT(IN) :: kqc_cutoff ! QC cutoff value 1372 1275 1373 !! * Local declarations 1276 1374 INTEGER :: jprof … … 1282 1380 DO jprof = 1, profdata%nprof 1283 1381 1284 IF ( profdata%nqc(jprof) > 10) THEN1382 IF ( profdata%nqc(jprof) > kqc_cutoff ) THEN 1285 1383 1286 1384 DO jvar = 1, profdata%nvar … … 1290 1388 1291 1389 profdata%var(jvar)%nvqc(jobs) = & 1292 & profdata%var(jvar)%nvqc(jobs) + 261390 & IBSET(profdata%var(jvar)%nvqc(jobs),14) 1293 1391 1294 1392 END DO … … 1302 1400 END SUBROUTINE obs_pro_rej 1303 1401 1304 SUBROUTINE obs_uv_rej( profdata, knumu, knumv )1402 SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff ) 1305 1403 !!---------------------------------------------------------------------- 1306 1404 !! *** ROUTINE obs_uv_rej *** … … 1322 1420 INTEGER, INTENT(INOUT) :: knumu ! Number of u rejected 1323 1421 INTEGER, INTENT(INOUT) :: knumv ! Number of v rejected 1422 INTEGER, INTENT(IN) :: kqc_cutoff ! QC cutoff value 1423 1324 1424 !! * Local declarations 1325 1425 INTEGER :: jprof … … 1341 1441 DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1) 1342 1442 1343 IF ( ( profdata%var(1)%nvqc(jobs) > 10) .AND. &1344 & ( profdata%var(2)%nvqc(jobs) <= 10) ) THEN1345 profdata%var(2)%nvqc(jobs) = profdata%var(2)%nvqc(jobs) + 421443 IF ( ( profdata%var(1)%nvqc(jobs) > kqc_cutoff ) .AND. & 1444 & ( profdata%var(2)%nvqc(jobs) <= kqc_cutoff) ) THEN 1445 profdata%var(2)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15) 1346 1446 knumv = knumv + 1 1347 1447 ENDIF 1348 IF ( ( profdata%var(2)%nvqc(jobs) > 10) .AND. &1349 & ( profdata%var(1)%nvqc(jobs) <= 10) ) THEN1350 profdata%var(1)%nvqc(jobs) = profdata%var(1)%nvqc(jobs) + 421448 IF ( ( profdata%var(2)%nvqc(jobs) > kqc_cutoff ) .AND. & 1449 & ( profdata%var(1)%nvqc(jobs) <= kqc_cutoff) ) THEN 1450 profdata%var(1)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15) 1351 1451 knumu = knumu + 1 1352 1452 ENDIF -
NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_profiles_def.F90
r11350 r11361 104 104 ! Bookkeeping arrays with sizes equal to number of variables 105 105 106 CHARACTER(len= 6), POINTER, DIMENSION(:) :: &106 CHARACTER(len=8), POINTER, DIMENSION(:) :: & 107 107 & cvars !: Variable names 108 108 -
NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_altbias.F90
r11350 r11361 128 128 ! Get the Alt bias data 129 129 130 CALL iom_get( numaltbias, jpdom_ data, 'altbias', z_altbias(:,:), 1 )130 CALL iom_get( numaltbias, jpdom_autoglo, 'altbias', z_altbias(:,:), 1 ) 131 131 132 132 ! Close the file -
NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90
r11350 r11361 45 45 SUBROUTINE obs_rea_prof( profdata, knumfiles, cdfilenames, & 46 46 & kvars, kextr, kstp, ddobsini, ddobsend, & 47 & ldvar 1, ldvar2, ldignmis, ldsatt, &47 & ldvar, ldignmis, ldsatt, & 48 48 & ldmod, kdailyavtypes ) 49 49 !!--------------------------------------------------------------------- … … 74 74 INTEGER, INTENT(IN) :: kextr ! Number of extra fields for each var 75 75 INTEGER, INTENT(IN) :: kstp ! Ocean time-step index 76 LOGICAL, INTENT(IN) :: ldvar1 ! Observed variables switches 77 LOGICAL, INTENT(IN) :: ldvar2 76 LOGICAL, DIMENSION(kvars), INTENT(IN) :: ldvar ! Observed variables switches 78 77 LOGICAL, INTENT(IN) :: ldignmis ! Ignore missing files 79 78 LOGICAL, INTENT(IN) :: ldsatt ! Compute salinity at all temperature points … … 87 86 CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_prof' 88 87 CHARACTER(len=8) :: clrefdate 89 CHARACTER(len= 6), DIMENSION(:), ALLOCATABLE :: clvars88 CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars 90 89 INTEGER :: jvar 91 90 INTEGER :: ji … … 105 104 INTEGER :: iprof 106 105 INTEGER :: iproftot 107 INTEGER :: ivar1t0 108 INTEGER :: ivar2t0 109 INTEGER :: ivar1t 110 INTEGER :: ivar2t 106 INTEGER, DIMENSION(kvars) :: ivart0 107 INTEGER, DIMENSION(kvars) :: ivart 111 108 INTEGER :: ip3dt 112 109 INTEGER :: ios 113 110 INTEGER :: ioserrcount 114 INTEGER :: ivar1tmpp 115 INTEGER :: ivar2tmpp 111 INTEGER, DIMENSION(kvars) :: ivartmpp 116 112 INTEGER :: ip3dtmpp 117 113 INTEGER :: itype 118 114 INTEGER, DIMENSION(knumfiles) :: & 119 115 & irefdate 120 INTEGER, DIMENSION(ntyp1770+1) :: & 121 & itypvar1, & 122 & itypvar1mpp, & 123 & itypvar2, & 124 & itypvar2mpp 116 INTEGER, DIMENSION(ntyp1770+1,kvars) :: & 117 & itypvar, & 118 & itypvarmpp 119 INTEGER, DIMENSION(:,:), ALLOCATABLE :: & 120 & iobsi, & 121 & iobsj, & 122 & iproc 125 123 INTEGER, DIMENSION(:), ALLOCATABLE :: & 126 & iobsi1, &127 & iobsj1, &128 & iproc1, &129 & iobsi2, &130 & iobsj2, &131 & iproc2, &132 124 & iindx, & 133 125 & ifileidx, & … … 147 139 LOGICAL :: llvalprof 148 140 LOGICAL :: lldavtimset 141 LOGICAL :: llcycle 149 142 TYPE(obfbdata), POINTER, DIMENSION(:) :: & 150 143 & inpfiles … … 152 145 ! Local initialization 153 146 iprof = 0 154 ivar1t0 = 0 155 ivar2t0 = 0 147 ivart0(:) = 0 156 148 ip3dt = 0 157 149 … … 219 211 & ldgrid = .TRUE. ) 220 212 221 IF ( inpfiles(jj)%nvar < 2) THEN213 IF ( inpfiles(jj)%nvar /= kvars ) THEN 222 214 CALL ctl_stop( 'Feedback format error: ', & 223 & ' less than 2vars in profile file' )215 & ' unexpected number of vars in profile file' ) 224 216 ENDIF 225 217 … … 307 299 inowin = 0 308 300 DO ji = 1, inpfiles(jj)%nobs 309 IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 310 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 311 & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE 301 IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 302 llcycle = .TRUE. 303 DO jvar = 1, kvars 304 IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 305 llcycle = .FALSE. 306 EXIT 307 ENDIF 308 END DO 309 IF ( llcycle ) CYCLE 312 310 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 313 311 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN … … 317 315 ALLOCATE( zlam(inowin) ) 318 316 ALLOCATE( zphi(inowin) ) 319 ALLOCATE( iobsi1(inowin) ) 320 ALLOCATE( iobsj1(inowin) ) 321 ALLOCATE( iproc1(inowin) ) 322 ALLOCATE( iobsi2(inowin) ) 323 ALLOCATE( iobsj2(inowin) ) 324 ALLOCATE( iproc2(inowin) ) 317 ALLOCATE( iobsi(inowin,kvars) ) 318 ALLOCATE( iobsj(inowin,kvars) ) 319 ALLOCATE( iproc(inowin,kvars) ) 325 320 inowin = 0 326 321 DO ji = 1, inpfiles(jj)%nobs 327 IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 328 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 329 & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE 322 IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 323 llcycle = .TRUE. 324 DO jvar = 1, kvars 325 IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 326 llcycle = .FALSE. 327 EXIT 328 ENDIF 329 END DO 330 IF ( llcycle ) CYCLE 330 331 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 331 332 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN … … 336 337 END DO 337 338 338 IF ( TRIM( inpfiles(jj)%cname(1) ) == 'POTM' ) THEN 339 CALL obs_grid_search( inowin, zlam, zphi, iobsi1, iobsj1, & 340 & iproc1, 'T' ) 341 iobsi2(:) = iobsi1(:) 342 iobsj2(:) = iobsj1(:) 343 iproc2(:) = iproc1(:) 344 ELSEIF ( TRIM( inpfiles(jj)%cname(1) ) == 'UVEL' ) THEN 345 CALL obs_grid_search( inowin, zlam, zphi, iobsi1, iobsj1, & 346 & iproc1, 'U' ) 347 CALL obs_grid_search( inowin, zlam, zphi, iobsi2, iobsj2, & 348 & iproc2, 'V' ) 339 ! Assume anything other than velocity is on T grid 340 IF ( TRIM( inpfiles(jj)%cname(1) ) == 'UVEL' ) THEN 341 CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,1), iobsj(:,1), & 342 & iproc(:,1), 'U' ) 343 CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,2), iobsj(:,2), & 344 & iproc(:,2), 'V' ) 345 ELSE 346 CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,1), iobsj(:,1), & 347 & iproc(:,1), 'T' ) 348 IF ( kvars > 1 ) THEN 349 DO jvar = 2, kvars 350 iobsi(:,jvar) = iobsi(:,1) 351 iobsj(:,jvar) = iobsj(:,1) 352 iproc(:,jvar) = iproc(:,1) 353 END DO 354 ENDIF 349 355 ENDIF 350 356 351 357 inowin = 0 352 358 DO ji = 1, inpfiles(jj)%nobs 353 IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 354 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 355 & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE 359 IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 360 llcycle = .TRUE. 361 DO jvar = 1, kvars 362 IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 363 llcycle = .FALSE. 364 EXIT 365 ENDIF 366 END DO 367 IF ( llcycle ) CYCLE 356 368 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 357 369 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN 358 370 inowin = inowin + 1 359 inpfiles(jj)%iproc(ji,1) = iproc1(inowin) 360 inpfiles(jj)%iobsi(ji,1) = iobsi1(inowin) 361 inpfiles(jj)%iobsj(ji,1) = iobsj1(inowin) 362 inpfiles(jj)%iproc(ji,2) = iproc2(inowin) 363 inpfiles(jj)%iobsi(ji,2) = iobsi2(inowin) 364 inpfiles(jj)%iobsj(ji,2) = iobsj2(inowin) 365 IF ( inpfiles(jj)%iproc(ji,1) /= & 366 & inpfiles(jj)%iproc(ji,2) ) THEN 367 CALL ctl_stop( 'Error in obs_read_prof:', & 368 & 'var1 and var2 observation on different processors') 371 DO jvar = 1, kvars 372 inpfiles(jj)%iproc(ji,jvar) = iproc(inowin,jvar) 373 inpfiles(jj)%iobsi(ji,jvar) = iobsi(inowin,jvar) 374 inpfiles(jj)%iobsj(ji,jvar) = iobsj(inowin,jvar) 375 END DO 376 IF ( kvars > 1 ) THEN 377 DO jvar = 2, kvars 378 IF ( inpfiles(jj)%iproc(ji,jvar) /= & 379 & inpfiles(jj)%iproc(ji,1) ) THEN 380 CALL ctl_stop( 'Error in obs_read_prof:', & 381 & 'observation on different processors for different vars') 382 ENDIF 383 END DO 369 384 ENDIF 370 385 ENDIF 371 386 END DO 372 DEALLOCATE( zlam, zphi, iobsi 1, iobsj1, iproc1, iobsi2, iobsj2, iproc2)387 DEALLOCATE( zlam, zphi, iobsi, iobsj, iproc ) 373 388 374 389 DO ji = 1, inpfiles(jj)%nobs 375 IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 376 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 377 & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE 390 IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 391 llcycle = .TRUE. 392 DO jvar = 1, kvars 393 IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 394 llcycle = .FALSE. 395 EXIT 396 ENDIF 397 END DO 398 IF ( llcycle ) CYCLE 378 399 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 379 400 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN … … 384 405 ENDIF 385 406 llvalprof = .FALSE. 386 IF ( ldvar1 ) THEN 387 loop_t_count : DO ij = 1,inpfiles(jj)%nlev 388 IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 389 & CYCLE 390 IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 391 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 392 ivar1t0 = ivar1t0 + 1 393 ENDIF 394 END DO loop_t_count 395 ENDIF 396 IF ( ldvar2 ) THEN 397 loop_s_count : DO ij = 1,inpfiles(jj)%nlev 398 IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 399 & CYCLE 400 IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 401 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 402 ivar2t0 = ivar2t0 + 1 403 ENDIF 404 END DO loop_s_count 405 ENDIF 406 loop_p_count : DO ij = 1,inpfiles(jj)%nlev 407 DO jvar = 1, kvars 408 IF ( ldvar(jvar) ) THEN 409 DO ij = 1,inpfiles(jj)%nlev 410 IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 411 & CYCLE 412 IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. & 413 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 414 ivart0(jvar) = ivart0(jvar) + 1 415 ENDIF 416 END DO 417 ENDIF 418 END DO 419 DO ij = 1,inpfiles(jj)%nlev 407 420 IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 408 421 & CYCLE 409 IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. &410 & ( inpfiles(jj)%idqc(ij,ji) <= 2 )) .AND. &411 & ldvar1 ) .OR. &412 & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. &413 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. &414 & ldvar2 ) ) THEN415 ip3dt = ip3dt + 1416 llvalprof = .TRUE.417 END IF418 END DO loop_p_count422 DO jvar = 1, kvars 423 IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. & 424 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 425 & ldvar(jvar) ) ) THEN 426 ip3dt = ip3dt + 1 427 llvalprof = .TRUE. 428 EXIT 429 ENDIF 430 END DO 431 END DO 419 432 420 433 IF ( llvalprof ) iprof = iprof + 1 … … 437 450 DO jj = 1, inobf 438 451 DO ji = 1, inpfiles(jj)%nobs 439 IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 440 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 441 & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE 452 IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 453 llcycle = .TRUE. 454 DO jvar = 1, kvars 455 IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 456 llcycle = .FALSE. 457 EXIT 458 ENDIF 459 END DO 460 IF ( llcycle ) CYCLE 442 461 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 443 462 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN … … 452 471 DO jj = 1, inobf 453 472 DO ji = 1, inpfiles(jj)%nobs 454 IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 455 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 456 & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE 473 IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 474 llcycle = .TRUE. 475 DO jvar = 1, kvars 476 IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 477 llcycle = .FALSE. 478 EXIT 479 ENDIF 480 END DO 481 IF ( llcycle ) CYCLE 457 482 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 458 483 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN … … 470 495 iv3dt(:) = -1 471 496 IF (ldsatt) THEN 472 iv3dt(1) = ip3dt 473 iv3dt(2) = ip3dt 497 iv3dt(:) = ip3dt 474 498 ELSE 475 iv3dt(1) = ivar1t0 476 iv3dt(2) = ivar2t0 499 iv3dt(:) = ivart0(:) 477 500 ENDIF 478 501 CALL obs_prof_alloc( profdata, kvars, kextr, iprof, iv3dt, & … … 487 510 488 511 ip3dt = 0 489 ivar1t = 0 490 ivar2t = 0 491 itypvar1 (:) = 0 492 itypvar1mpp(:) = 0 493 494 itypvar2 (:) = 0 495 itypvar2mpp(:) = 0 512 ivart(:) = 0 513 itypvar (:,:) = 0 514 itypvarmpp(:,:) = 0 496 515 497 516 ioserrcount = 0 … … 501 520 ji = iprofidx(iindx(jk)) 502 521 503 IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 504 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 505 & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE 522 IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 523 llcycle = .TRUE. 524 DO jvar = 1, kvars 525 IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 526 llcycle = .FALSE. 527 EXIT 528 ENDIF 529 END DO 530 IF ( llcycle ) CYCLE 506 531 507 532 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & … … 518 543 IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 519 544 520 IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 521 & ( inpfiles(jj)%ivqc(ji,2) > 2 ) ) CYCLE 545 IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 546 llcycle = .TRUE. 547 DO jvar = 1, kvars 548 IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 549 llcycle = .FALSE. 550 EXIT 551 ENDIF 552 END DO 553 IF ( llcycle ) CYCLE 522 554 523 555 loop_prof : DO ij = 1, inpfiles(jj)%nlev … … 526 558 & CYCLE 527 559 528 IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 529 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 530 531 llvalprof = .TRUE. 532 EXIT loop_prof 533 534 ENDIF 535 536 IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 537 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 538 539 llvalprof = .TRUE. 540 EXIT loop_prof 541 542 ENDIF 560 DO jvar = 1, kvars 561 IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. & 562 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 563 564 llvalprof = .TRUE. 565 EXIT loop_prof 566 567 ENDIF 568 END DO 543 569 544 570 END DO loop_prof … … 572 598 573 599 ! Coordinate search parameters 574 profdata%mi (iprof,1) = inpfiles(jj)%iobsi(ji,1)575 profdata%mj (iprof,1) = inpfiles(jj)%iobsj(ji,1)576 profdata%mi (iprof,2) = inpfiles(jj)%iobsi(ji,2)577 profdata%mj (iprof,2) = inpfiles(jj)%iobsj(ji,2)600 DO jvar = 1, kvars 601 profdata%mi (iprof,jvar) = inpfiles(jj)%iobsi(ji,jvar) 602 profdata%mj (iprof,jvar) = inpfiles(jj)%iobsj(ji,jvar) 603 END DO 578 604 579 605 ! Profile WMO number … … 615 641 IF (ldsatt) THEN 616 642 617 IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 618 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 619 & ldvar1 ) .OR. & 620 & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 621 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 622 & ldvar2 ) ) THEN 623 ip3dt = ip3dt + 1 624 ELSE 625 CYCLE 643 DO jvar = 1, kvars 644 IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. & 645 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 646 & ldvar(jvar) ) ) THEN 647 ip3dt = ip3dt + 1 648 EXIT 649 ELSE IF ( jvar == kvars ) THEN 650 CYCLE loop_p 651 ENDIF 652 END DO 653 654 ENDIF 655 656 DO jvar = 1, kvars 657 658 IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. & 659 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 660 & ldvar(jvar) ) .OR. ldsatt ) THEN 661 662 IF (ldsatt) THEN 663 664 ivart(jvar) = ip3dt 665 666 ELSE 667 668 ivart(jvar) = ivart(jvar) + 1 669 670 ENDIF 671 672 ! Depth of jvar observation 673 profdata%var(jvar)%vdep(ivart(jvar)) = & 674 & inpfiles(jj)%pdep(ij,ji) 675 676 ! Depth of jvar observation QC 677 profdata%var(jvar)%idqc(ivart(jvar)) = & 678 & inpfiles(jj)%idqc(ij,ji) 679 680 ! Depth of jvar observation QC flags 681 profdata%var(jvar)%idqcf(:,ivart(jvar)) = & 682 & inpfiles(jj)%idqcf(:,ij,ji) 683 684 ! Profile index 685 profdata%var(jvar)%nvpidx(ivart(jvar)) = iprof 686 687 ! Vertical index in original profile 688 profdata%var(jvar)%nvlidx(ivart(jvar)) = ij 689 690 ! Profile jvar value 691 IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. & 692 & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 693 profdata%var(jvar)%vobs(ivart(jvar)) = & 694 & inpfiles(jj)%pob(ij,ji,jvar) 695 IF ( ldmod ) THEN 696 profdata%var(jvar)%vmod(ivart(jvar)) = & 697 & inpfiles(jj)%padd(ij,ji,1,jvar) 698 ENDIF 699 ! Count number of profile var1 data as function of type 700 itypvar( profdata%ntyp(iprof) + 1, jvar ) = & 701 & itypvar( profdata%ntyp(iprof) + 1, jvar ) + 1 702 ELSE 703 profdata%var(jvar)%vobs(ivart(jvar)) = fbrmdi 704 ENDIF 705 706 ! Profile jvar qc 707 profdata%var(jvar)%nvqc(ivart(jvar)) = & 708 & inpfiles(jj)%ivlqc(ij,ji,jvar) 709 710 ! Profile jvar qc flags 711 profdata%var(jvar)%nvqcf(:,ivart(jvar)) = & 712 & inpfiles(jj)%ivlqcf(:,ij,ji,jvar) 713 714 ! Profile insitu T value 715 IF ( TRIM( inpfiles(jj)%cname(jvar) ) == 'POTM' ) THEN 716 profdata%var(jvar)%vext(ivart(jvar),1) = & 717 & inpfiles(jj)%pext(ij,ji,1) 718 ENDIF 719 626 720 ENDIF 627 628 ENDIF 629 630 IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 631 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 632 & ldvar1 ) .OR. ldsatt ) THEN 633 634 IF (ldsatt) THEN 635 636 ivar1t = ip3dt 637 638 ELSE 639 640 ivar1t = ivar1t + 1 641 642 ENDIF 643 644 ! Depth of var1 observation 645 profdata%var(1)%vdep(ivar1t) = & 646 & inpfiles(jj)%pdep(ij,ji) 647 648 ! Depth of var1 observation QC 649 profdata%var(1)%idqc(ivar1t) = & 650 & inpfiles(jj)%idqc(ij,ji) 651 652 ! Depth of var1 observation QC flags 653 profdata%var(1)%idqcf(:,ivar1t) = & 654 & inpfiles(jj)%idqcf(:,ij,ji) 655 656 ! Profile index 657 profdata%var(1)%nvpidx(ivar1t) = iprof 658 659 ! Vertical index in original profile 660 profdata%var(1)%nvlidx(ivar1t) = ij 661 662 ! Profile var1 value 663 IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 664 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 665 profdata%var(1)%vobs(ivar1t) = & 666 & inpfiles(jj)%pob(ij,ji,1) 667 IF ( ldmod ) THEN 668 profdata%var(1)%vmod(ivar1t) = & 669 & inpfiles(jj)%padd(ij,ji,1,1) 670 ENDIF 671 ! Count number of profile var1 data as function of type 672 itypvar1( profdata%ntyp(iprof) + 1 ) = & 673 & itypvar1( profdata%ntyp(iprof) + 1 ) + 1 674 ELSE 675 profdata%var(1)%vobs(ivar1t) = fbrmdi 676 ENDIF 677 678 ! Profile var1 qc 679 profdata%var(1)%nvqc(ivar1t) = & 680 & inpfiles(jj)%ivlqc(ij,ji,1) 681 682 ! Profile var1 qc flags 683 profdata%var(1)%nvqcf(:,ivar1t) = & 684 & inpfiles(jj)%ivlqcf(:,ij,ji,1) 685 686 ! Profile insitu T value 687 IF ( TRIM( inpfiles(jj)%cname(1) ) == 'POTM' ) THEN 688 profdata%var(1)%vext(ivar1t,1) = & 689 & inpfiles(jj)%pext(ij,ji,1) 690 ENDIF 691 692 ENDIF 693 694 IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 695 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 696 & ldvar2 ) .OR. ldsatt ) THEN 697 698 IF (ldsatt) THEN 699 700 ivar2t = ip3dt 701 702 ELSE 703 704 ivar2t = ivar2t + 1 705 706 ENDIF 707 708 ! Depth of var2 observation 709 profdata%var(2)%vdep(ivar2t) = & 710 & inpfiles(jj)%pdep(ij,ji) 711 712 ! Depth of var2 observation QC 713 profdata%var(2)%idqc(ivar2t) = & 714 & inpfiles(jj)%idqc(ij,ji) 715 716 ! Depth of var2 observation QC flags 717 profdata%var(2)%idqcf(:,ivar2t) = & 718 & inpfiles(jj)%idqcf(:,ij,ji) 719 720 ! Profile index 721 profdata%var(2)%nvpidx(ivar2t) = iprof 722 723 ! Vertical index in original profile 724 profdata%var(2)%nvlidx(ivar2t) = ij 725 726 ! Profile var2 value 727 IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 728 & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 729 profdata%var(2)%vobs(ivar2t) = & 730 & inpfiles(jj)%pob(ij,ji,2) 731 IF ( ldmod ) THEN 732 profdata%var(2)%vmod(ivar2t) = & 733 & inpfiles(jj)%padd(ij,ji,1,2) 734 ENDIF 735 ! Count number of profile var2 data as function of type 736 itypvar2( profdata%ntyp(iprof) + 1 ) = & 737 & itypvar2( profdata%ntyp(iprof) + 1 ) + 1 738 ELSE 739 profdata%var(2)%vobs(ivar2t) = fbrmdi 740 ENDIF 741 742 ! Profile var2 qc 743 profdata%var(2)%nvqc(ivar2t) = & 744 & inpfiles(jj)%ivlqc(ij,ji,2) 745 746 ! Profile var2 qc flags 747 profdata%var(2)%nvqcf(:,ivar2t) = & 748 & inpfiles(jj)%ivlqcf(:,ij,ji,2) 749 750 ENDIF 721 722 END DO 751 723 752 724 END DO loop_p … … 762 734 !----------------------------------------------------------------------- 763 735 764 CALL obs_mpp_sum_integer ( ivar1t0, ivar1tmpp ) 765 CALL obs_mpp_sum_integer ( ivar2t0, ivar2tmpp ) 736 DO jvar = 1, kvars 737 CALL obs_mpp_sum_integer ( ivart0(jvar), ivartmpp(jvar) ) 738 END DO 766 739 CALL obs_mpp_sum_integer ( ip3dt, ip3dtmpp ) 767 740 768 CALL obs_mpp_sum_integers( itypvar1, itypvar1mpp, ntyp1770 + 1 ) 769 CALL obs_mpp_sum_integers( itypvar2, itypvar2mpp, ntyp1770 + 1 ) 741 DO jvar = 1, kvars 742 CALL obs_mpp_sum_integers( itypvar(:,jvar), itypvarmpp(:,jvar), ntyp1770 + 1 ) 743 END DO 770 744 771 745 !----------------------------------------------------------------------- … … 777 751 WRITE(numout,'(1X,A)') '------------' 778 752 WRITE(numout,*) 779 WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(1) ) 780 WRITE(numout,'(1X,A)') '------------------------' 781 DO ji = 0, ntyp1770 782 IF ( itypvar1mpp(ji+1) > 0 ) THEN 783 WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), & 784 & cwmonam1770(ji)(1:52),' = ', & 785 & itypvar1mpp(ji+1) 786 ENDIF 753 DO jvar = 1, kvars 754 WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(jvar) ) 755 WRITE(numout,'(1X,A)') '------------------------' 756 DO ji = 0, ntyp1770 757 IF ( itypvarmpp(ji+1,jvar) > 0 ) THEN 758 WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), & 759 & cwmonam1770(ji)(1:52),' = ', & 760 & itypvarmpp(ji+1,jvar) 761 ENDIF 762 END DO 763 WRITE(numout,'(1X,A)') & 764 & '---------------------------------------------------------------' 765 WRITE(numout,'(1X,A55,I8)') & 766 & 'Total profile data for variable '//TRIM( profdata%cvars(jvar) )// & 767 & ' = ', ivartmpp(jvar) 768 WRITE(numout,'(1X,A)') & 769 & '---------------------------------------------------------------' 770 WRITE(numout,*) 787 771 END DO 788 WRITE(numout,'(1X,A)') & 789 & '---------------------------------------------------------------' 790 WRITE(numout,'(1X,A55,I8)') & 791 & 'Total profile data for variable '//TRIM( profdata%cvars(1) )// & 792 & ' = ', ivar1tmpp 793 WRITE(numout,'(1X,A)') & 794 & '---------------------------------------------------------------' 795 WRITE(numout,*) 796 WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(2) ) 797 WRITE(numout,'(1X,A)') '------------------------' 798 DO ji = 0, ntyp1770 799 IF ( itypvar2mpp(ji+1) > 0 ) THEN 800 WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), & 801 & cwmonam1770(ji)(1:52),' = ', & 802 & itypvar2mpp(ji+1) 803 ENDIF 772 ENDIF 773 774 IF (ldsatt) THEN 775 profdata%nvprot(:) = ip3dt 776 profdata%nvprotmpp(:) = ip3dtmpp 777 ELSE 778 DO jvar = 1, kvars 779 profdata%nvprot(jvar) = ivart(jvar) 780 profdata%nvprotmpp(jvar) = ivartmpp(jvar) 804 781 END DO 805 WRITE(numout,'(1X,A)') &806 & '---------------------------------------------------------------'807 WRITE(numout,'(1X,A55,I8)') &808 & 'Total profile data for variable '//TRIM( profdata%cvars(2) )// &809 & ' = ', ivar2tmpp810 WRITE(numout,'(1X,A)') &811 & '---------------------------------------------------------------'812 WRITE(numout,*)813 ENDIF814 815 IF (ldsatt) THEN816 profdata%nvprot(1) = ip3dt817 profdata%nvprot(2) = ip3dt818 profdata%nvprotmpp(1) = ip3dtmpp819 profdata%nvprotmpp(2) = ip3dtmpp820 ELSE821 profdata%nvprot(1) = ivar1t822 profdata%nvprot(2) = ivar2t823 profdata%nvprotmpp(1) = ivar1tmpp824 profdata%nvprotmpp(2) = ivar2tmpp825 782 ENDIF 826 783 profdata%nprof = iprof … … 829 786 ! Model level search 830 787 !----------------------------------------------------------------------- 831 IF ( ldvar1 ) THEN 832 CALL obs_level_search( jpk, gdept_1d, & 833 & profdata%nvprot(1), profdata%var(1)%vdep, & 834 & profdata%var(1)%mvk ) 835 ENDIF 836 IF ( ldvar2 ) THEN 837 CALL obs_level_search( jpk, gdept_1d, & 838 & profdata%nvprot(2), profdata%var(2)%vdep, & 839 & profdata%var(2)%mvk ) 840 ENDIF 788 DO jvar = 1, kvars 789 IF ( ldvar(jvar) ) THEN 790 CALL obs_level_search( jpk, gdept_1d, & 791 & profdata%nvprot(jvar), profdata%var(jvar)%vdep, & 792 & profdata%var(jvar)%mvk ) 793 ENDIF 794 END DO 841 795 842 796 !----------------------------------------------------------------------- -
NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90
r11350 r11361 77 77 CHARACTER(LEN=11), PARAMETER :: cpname='obs_rea_surf' 78 78 CHARACTER(len=8) :: clrefdate 79 CHARACTER(len= 6), DIMENSION(:), ALLOCATABLE :: clvars79 CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars 80 80 INTEGER :: ji 81 81 INTEGER :: jj … … 97 97 INTEGER :: ios 98 98 INTEGER :: ioserrcount 99 INTEGER :: iextr 99 100 INTEGER, PARAMETER :: jpsurfmaxtype = 1024 100 101 INTEGER, DIMENSION(knumfiles) :: irefdate … … 108 109 & iindx, & 109 110 & ifileidx, & 110 & isurfidx 111 & isurfidx, & 112 & iadd_std 111 113 REAL(wp), DIMENSION(:), ALLOCATABLE :: & 112 114 & zphi, & … … 122 124 123 125 ! Local initialization 124 iobs = 0125 126 iobs = 0 127 iextr = kextr 126 128 !----------------------------------------------------------------------- 127 129 ! Count the number of files needed and allocate the obfbdata type … … 131 133 132 134 ALLOCATE( inpfiles(inobf) ) 135 ALLOCATE( iadd_std(inobf) ) 133 136 134 137 surf_files : DO jj = 1, inobf … … 172 175 173 176 !------------------------------------------------------------------ 174 ! Read the profile file into inpfiles177 ! Read the surface file into inpfiles 175 178 !------------------------------------------------------------------ 176 179 CALL init_obfbdata( inpfiles(jj) ) … … 197 200 ENDIF 198 201 199 IF (lwp) WRITE(numout,*)'Observation file contains ',inpfiles(jj)%nobs,' observations' 202 iadd_std(jj) = -1 203 IF ( inpfiles(jj)%nadd > 0 ) THEN 204 DO ji = 1, inpfiles(jj)%nadd 205 IF ( TRIM( inpfiles(jj)%caddname(ji) ) == 'STD' ) THEN 206 iextr = kextr + 1 207 iadd_std(jj) = ji 208 ENDIF 209 END DO 210 ENDIF 211 212 IF(lwp) THEN 213 IF ( iadd_std(jj) /= -1 ) THEN 214 WRITE(numout,*) ' STD variable available in input file so passing it through the obs oper' 215 WRITE(numout,*) 216 ENDIF 217 ENDIF 200 218 201 219 !------------------------------------------------------------------ … … 296 314 ENDIF 297 315 llvalprof = .FALSE. 298 IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. & 299 & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN 316 IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(1,ji,1),2) ) THEN 300 317 iobs = iobs + 1 301 318 ENDIF … … 342 359 & iindx ) 343 360 344 CALL obs_surf_alloc( surfdata, iobs, kvars, kextr, kstp, jpi, jpj )361 CALL obs_surf_alloc( surfdata, iobs, kvars, iextr, kstp, jpi, jpj ) 345 362 346 363 ! Read obs/positions, QC, all variable and assign to surfdata 347 364 348 365 iobs = 0 349 350 366 surfdata%cvars(:) = clvars(:) 367 IF ( ldmod .AND. ( TRIM( surfdata%cvars(1) ) == 'SLA' ) ) THEN 368 surfdata%cext(1) = 'SSH' 369 surfdata%cext(2) = 'MDT' 370 ENDIF 371 IF ( iextr > kextr ) surfdata%cext(iextr) = 'STD' 351 372 352 373 ityp (:) = 0 … … 370 391 ! Set observation information 371 392 372 IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. & 373 & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN 393 IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(1,ji,1),2) ) THEN 374 394 375 395 iobs = iobs + 1 … … 431 451 surfdata%robs(iobs,1) = inpfiles(jj)%pob(1,ji,1) 432 452 433 434 453 ! Model and MDT is set to fbrmdi unless read from file 435 454 IF ( ldmod ) THEN … … 442 461 surfdata%rmod(iobs,1) = fbrmdi 443 462 IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) surfdata%rext(iobs,:) = fbrmdi 463 ENDIF 464 465 ! STD (obs error standard deviation) read from file and passed through obs operator 466 IF ( iadd_std(jj) /= -1 ) THEN 467 surfdata%rext(iobs,iextr) = inpfiles(jj)%padd(1,ji,iadd_std(jj),1) 444 468 ENDIF 445 469 ENDIF … … 493 517 END DO 494 518 DEALLOCATE( inpfiles ) 519 DEALLOCATE( iadd_std ) 495 520 496 521 END SUBROUTINE obs_rea_surf -
NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_readmdt.F90
r11350 r11361 24 24 USE lib_mpp ! MPP library 25 25 USE dom_oce, ONLY : & ! Domain variables 26 & tmask, tmask_i, e1 e2t, gphit, glamt26 & tmask, tmask_i, e1t, e2t, gphit, glamt 27 27 USE obs_const, ONLY : obfillflt ! Fillvalue 28 28 USE oce , ONLY : sshn ! Model variables … … 154 154 155 155 ! mark any masked data with a QC flag 156 IF( zobsmask(1) == 0 ) sladata%nqc(jobs) = 11156 IF( zobsmask(1) == 0 ) sladata%nqc(jobs) = IBSET(sladata%nqc(jobs),15) 157 157 158 158 END DO … … 219 219 DO jj = 1, jpj 220 220 DO ji = 1, jpi 221 zdxdy = e1 e2t(ji,jj) * zpromsk(ji,jj)221 zdxdy = e1t(ji,jj) * e2t(ji,jj) * zpromsk(ji,jj) 222 222 zarea = zarea + zdxdy 223 223 zeta1 = zeta1 + mdt(ji,jj) * zdxdy -
NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90
r11350 r11361 50 50 INTEGER :: npj 51 51 INTEGER :: nsurfup !: Observation counter used in obs_oper 52 INTEGER :: nrec !: Number of surface observation records in window 52 53 53 54 ! Arrays with size equal to the number of surface observations … … 56 57 & mi, & !: i-th grid coord. for interpolating to surface observation 57 58 & mj, & !: j-th grid coord. for interpolating to surface observation 59 & mt, & !: time record number for gridded data 58 60 & nsidx,& !: Surface observation number 59 61 & nsfil,& !: Surface observation number in file … … 67 69 & ntyp !: Type of surface observation product 68 70 69 CHARACTER(len= 6), POINTER, DIMENSION(:) :: &71 CHARACTER(len=8), POINTER, DIMENSION(:) :: & 70 72 & cvars !: Variable names 73 74 CHARACTER(len=8), POINTER, DIMENSION(:) :: & 75 & cext !: Extra field names 71 76 72 77 CHARACTER(LEN=8), POINTER, DIMENSION(:) :: & … … 93 98 & nsstpmpp !: Global number of surface observations per time step 94 99 100 ! Arrays with size equal to the number of observation records in the window 101 INTEGER, POINTER, DIMENSION(:) :: & 102 & mrecstp ! Time step of the records 103 95 104 ! Arrays used to store source indices when 96 105 ! compressing obs_surf derived types … … 100 109 INTEGER, POINTER, DIMENSION(:) :: & 101 110 & nsind !: Source indices of surface data in compressed data 111 112 ! Is this a gridded product? 113 114 LOGICAL :: lgrid 102 115 103 116 END TYPE obs_surf … … 160 173 & surf%mi(ksurf), & 161 174 & surf%mj(ksurf), & 175 & surf%mt(ksurf), & 162 176 & surf%nsidx(ksurf), & 163 177 & surf%nsfil(ksurf), & … … 176 190 & ) 177 191 192 surf%mt(:) = -1 193 178 194 179 195 ! Allocate arrays of number of surface data size * number of variables … … 187 203 188 204 ALLOCATE( & 189 & surf%rext(ksurf,kextra) & 190 & ) 205 & surf%rext(ksurf,kextra), & 206 & surf%cext(kextra) & 207 & ) 208 209 surf%rext(:,:) = 0.0_wp 210 211 DO ji = 1, kextra 212 surf%cext(ji) = "NotSet" 213 END DO 191 214 192 215 ! Allocate arrays of number of time step size … … 217 240 218 241 surf%nsurfup = 0 242 243 ! Not gridded by default 244 245 surf%lgrid = .FALSE. 219 246 220 247 END SUBROUTINE obs_surf_alloc … … 242 269 & surf%mi, & 243 270 & surf%mj, & 271 & surf%mt, & 244 272 & surf%nsidx, & 245 273 & surf%nsfil, & … … 268 296 269 297 DEALLOCATE( & 270 & surf%rext & 298 & surf%rext, & 299 & surf%cext & 271 300 & ) 272 301 … … 370 399 newsurf%mi(insurf) = surf%mi(ji) 371 400 newsurf%mj(insurf) = surf%mj(ji) 401 newsurf%mt(insurf) = surf%mt(ji) 372 402 newsurf%nsidx(insurf) = surf%nsidx(ji) 373 403 newsurf%nsfil(insurf) = surf%nsfil(ji) … … 414 444 newsurf%nstp = surf%nstp 415 445 newsurf%cvars(:) = surf%cvars(:) 446 newsurf%cext(:) = surf%cext(:) 447 448 ! Set gridded stuff 449 450 newsurf%mt(insurf) = surf%mt(ji) 416 451 417 452 ! Deallocate temporary data … … 454 489 oldsurf%mi(jj) = surf%mi(ji) 455 490 oldsurf%mj(jj) = surf%mj(ji) 491 oldsurf%mt(jj) = surf%mt(ji) 456 492 oldsurf%nsidx(jj) = surf%nsidx(ji) 457 493 oldsurf%nsfil(jj) = surf%nsfil(ji) -
NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90
r11350 r11361 8 8 !! obs_wri_prof : Write profile observations in feedback format 9 9 !! obs_wri_surf : Write surface observations in feedback format 10 !! obs_wri_stats : Print basic statistics on the data being written out10 !! obs_wri_stats : Print basic statistics on the data being written out 11 11 !!---------------------------------------------------------------------- 12 12 … … 83 83 TYPE(obfbdata) :: fbdata 84 84 CHARACTER(LEN=40) :: clfname 85 CHARACTER(LEN=6) :: clfiletype 85 CHARACTER(LEN=10) :: clfiletype 86 CHARACTER(LEN=ilenlong) :: cllongname ! Long name of variable 87 CHARACTER(LEN=ilenunit) :: clunits ! Units of variable 88 CHARACTER(LEN=ilengrid) :: clgrid ! Grid of variable 86 89 INTEGER :: ilevel 87 90 INTEGER :: jvar … … 111 114 ! Find maximum level 112 115 ilevel = 0 113 DO jvar = 1, 2116 DO jvar = 1, profdata%nvar 114 117 ilevel = MAX( ilevel, MAXVAL( profdata%var(jvar)%nvlidx(:) ) ) 115 118 END DO … … 175 178 END DO 176 179 180 CASE('PLCHLTOT') 181 182 clfiletype = 'plchltotfb' 183 cllongname = 'log10(chlorophyll concentration)' 184 clunits = 'log10(mg/m3)' 185 clgrid = 'T' 186 187 CASE('PCHLTOT') 188 189 clfiletype = 'pchltotfb' 190 cllongname = 'chlorophyll concentration' 191 clunits = 'mg/m3' 192 clgrid = 'T' 193 194 CASE('PNO3') 195 196 clfiletype = 'pno3fb' 197 cllongname = 'nitrate' 198 clunits = 'mmol/m3' 199 clgrid = 'T' 200 201 CASE('PSI4') 202 203 clfiletype = 'psi4fb' 204 cllongname = 'silicate' 205 clunits = 'mmol/m3' 206 clgrid = 'T' 207 208 CASE('PPO4') 209 210 clfiletype = 'ppo4fb' 211 cllongname = 'phosphate' 212 clunits = 'mmol/m3' 213 clgrid = 'T' 214 215 CASE('PDIC') 216 217 clfiletype = 'pdicfb' 218 cllongname = 'dissolved inorganic carbon' 219 clunits = 'mmol/m3' 220 clgrid = 'T' 221 222 CASE('PALK') 223 224 clfiletype = 'palkfb' 225 cllongname = 'alkalinity' 226 clunits = 'meq/m3' 227 clgrid = 'T' 228 229 CASE('PPH') 230 231 clfiletype = 'pphfb' 232 cllongname = 'pH' 233 clunits = '-' 234 clgrid = 'T' 235 236 CASE('PO2') 237 238 clfiletype = 'po2fb' 239 cllongname = 'dissolved oxygen' 240 clunits = 'mmol/m3' 241 clgrid = 'T' 242 177 243 END SELECT 244 245 IF ( ( TRIM(profdata%cvars(1)) /= 'POTM' ) .AND. & 246 & ( TRIM(profdata%cvars(1)) /= 'UVEL' ) ) THEN 247 CALL alloc_obfbdata( fbdata, 1, profdata%nprof, ilevel, & 248 & 1 + iadd, iext, .TRUE. ) 249 fbdata%cname(1) = profdata%cvars(1) 250 fbdata%coblong(1) = cllongname 251 fbdata%cobunit(1) = clunits 252 fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(cllongname) 253 fbdata%caddunit(1,1) = clunits 254 fbdata%cgrid(:) = clgrid 255 DO je = 1, iext 256 fbdata%cextname(je) = pext%cdname(je) 257 fbdata%cextlong(je) = pext%cdlong(je,1) 258 fbdata%cextunit(je) = pext%cdunit(je,1) 259 END DO 260 DO ja = 1, iadd 261 fbdata%caddname(1+ja) = padd%cdname(ja) 262 fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 263 fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 264 END DO 265 ENDIF 178 266 179 267 fbdata%caddname(1) = 'Hx' … … 196 284 fbdata%ivqc(jo,:) = profdata%ivqc(jo,:) 197 285 fbdata%ivqcf(:,jo,:) = profdata%ivqcf(:,jo,:) 198 IF ( profdata%nqc(jo) > 10) THEN199 fbdata%ioqc(jo) = 4286 IF ( profdata%nqc(jo) > 255 ) THEN 287 fbdata%ioqc(jo) = IBSET(profdata%nqc(jo),2) 200 288 fbdata%ioqcf(1,jo) = profdata%nqcf(1,jo) 201 fbdata%ioqcf(2,jo) = profdata%nqc(jo) - 10289 fbdata%ioqcf(2,jo) = profdata%nqc(jo) 202 290 ELSE 203 291 fbdata%ioqc(jo) = profdata%nqc(jo) … … 228 316 & krefdate = 19500101 ) 229 317 ! Reform the profiles arrays for output 230 DO jvar = 1, 2318 DO jvar = 1, profdata%nvar 231 319 DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar) 232 320 ik = profdata%var(jvar)%nvlidx(jk) … … 236 324 fbdata%idqc(ik,jo) = profdata%var(jvar)%idqc(jk) 237 325 fbdata%idqcf(:,ik,jo) = profdata%var(jvar)%idqcf(:,jk) 238 IF ( profdata%var(jvar)%nvqc(jk) > 10) THEN239 fbdata%ivlqc(ik,jo,jvar) = 4326 IF ( profdata%var(jvar)%nvqc(jk) > 255 ) THEN 327 fbdata%ivlqc(ik,jo,jvar) = IBSET(profdata%var(jvar)%nvqc(jk),2) 240 328 fbdata%ivlqcf(1,ik,jo,jvar) = profdata%var(jvar)%nvqcf(1,jk) 241 fbdata%ivlqcf(2,ik,jo,jvar) = profdata%var(jvar)%nvqc(jk) - 10329 fbdata%ivlqcf(2,ik,jo,jvar) = IAND(profdata%var(jvar)%nvqc(jk),b'0000 0000 1111 1111') 242 330 ELSE 243 331 fbdata%ivlqc(ik,jo,jvar) = profdata%var(jvar)%nvqc(jk) … … 320 408 TYPE(obfbdata) :: fbdata 321 409 CHARACTER(LEN=40) :: clfname ! netCDF filename 322 CHARACTER(LEN= 6):: clfiletype410 CHARACTER(LEN=10) :: clfiletype 323 411 CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf' 412 CHARACTER(LEN=ilenlong) :: cllongname ! Long name of variable 413 CHARACTER(LEN=ilenunit) :: clunits ! Units of variable 414 CHARACTER(LEN=ilengrid) :: clgrid ! Grid of variable 324 415 INTEGER :: jo 325 416 INTEGER :: ja … … 327 418 INTEGER :: iadd 328 419 INTEGER :: iext 420 INTEGER :: indx_std 421 INTEGER :: iadd_std 329 422 330 423 IF ( PRESENT( padd ) ) THEN … … 340 433 ENDIF 341 434 435 iadd_std = 0 436 indx_std = -1 437 IF ( surfdata%nextra > 0 ) THEN 438 DO je = 1, surfdata%nextra 439 IF ( TRIM( surfdata%cext(je) ) == 'STD' ) THEN 440 iadd_std = 1 441 indx_std = je 442 ENDIF 443 END DO 444 ENDIF 445 342 446 CALL init_obfbdata( fbdata ) 343 447 … … 345 449 CASE('SLA') 346 450 451 ! SLA needs special treatment because of MDT, so is all done here 452 ! Other variables are done more generically 453 347 454 CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 348 & 2 + iadd , 1 + iext, .TRUE. )455 & 2 + iadd + iadd_std, 1 + iext, .TRUE. ) 349 456 350 457 clfiletype = 'slafb' … … 367 474 fbdata%cgrid(1) = 'T' 368 475 DO ja = 1, iadd 369 fbdata%caddname(2+ ja) = padd%cdname(ja)370 fbdata%caddlong(2+ ja,1) = padd%cdlong(ja,1)371 fbdata%caddunit(2+ ja,1) = padd%cdunit(ja,1)476 fbdata%caddname(2+iadd_std+ja) = padd%cdname(ja) 477 fbdata%caddlong(2+iadd_std+ja,1) = padd%cdlong(ja,1) 478 fbdata%caddunit(2+iadd_std+ja,1) = padd%cdunit(ja,1) 372 479 END DO 373 480 374 481 CASE('SST') 375 482 483 clfiletype = 'sstfb' 484 cllongname = 'Sea surface temperature' 485 clunits = 'Degree centigrade' 486 clgrid = 'T' 487 488 CASE('ICECONC') 489 490 clfiletype = 'sicfb' 491 cllongname = 'Sea ice' 492 clunits = 'Fraction' 493 clgrid = 'T' 494 495 CASE('SSS') 496 497 clfiletype = 'sssfb' 498 cllongname = 'Sea surface salinity' 499 clunits = 'psu' 500 clgrid = 'T' 501 502 CASE('SLCHLTOT','LOGCHL','LogChl','logchl') 503 504 clfiletype = 'slchltotfb' 505 cllongname = 'Surface total log10(chlorophyll)' 506 clunits = 'log10(mg/m3)' 507 clgrid = 'T' 508 509 CASE('SLCHLDIA') 510 511 clfiletype = 'slchldiafb' 512 cllongname = 'Surface diatom log10(chlorophyll)' 513 clunits = 'log10(mg/m3)' 514 clgrid = 'T' 515 516 CASE('SLCHLNON') 517 518 clfiletype = 'slchlnonfb' 519 cllongname = 'Surface non-diatom log10(chlorophyll)' 520 clunits = 'log10(mg/m3)' 521 clgrid = 'T' 522 523 CASE('SLCHLDIN') 524 525 clfiletype = 'slchldinfb' 526 cllongname = 'Surface dinoflagellate log10(chlorophyll)' 527 clunits = 'log10(mg/m3)' 528 clgrid = 'T' 529 530 CASE('SLCHLMIC') 531 532 clfiletype = 'slchlmicfb' 533 cllongname = 'Surface microphytoplankton log10(chlorophyll)' 534 clunits = 'log10(mg/m3)' 535 clgrid = 'T' 536 537 CASE('SLCHLNAN') 538 539 clfiletype = 'slchlnanfb' 540 cllongname = 'Surface nanophytoplankton log10(chlorophyll)' 541 clunits = 'log10(mg/m3)' 542 clgrid = 'T' 543 544 CASE('SLCHLPIC') 545 546 clfiletype = 'slchlpicfb' 547 cllongname = 'Surface picophytoplankton log10(chlorophyll)' 548 clunits = 'log10(mg/m3)' 549 clgrid = 'T' 550 551 CASE('SCHLTOT') 552 553 clfiletype = 'schltotfb' 554 cllongname = 'Surface total chlorophyll' 555 clunits = 'mg/m3' 556 clgrid = 'T' 557 558 CASE('SLPHYTOT') 559 560 clfiletype = 'slphytotfb' 561 cllongname = 'Surface total log10(phytoplankton carbon)' 562 clunits = 'log10(mmolC/m3)' 563 clgrid = 'T' 564 565 CASE('SLPHYDIA') 566 567 clfiletype = 'slphydiafb' 568 cllongname = 'Surface diatom log10(phytoplankton carbon)' 569 clunits = 'log10(mmolC/m3)' 570 clgrid = 'T' 571 572 CASE('SLPHYNON') 573 574 clfiletype = 'slphynonfb' 575 cllongname = 'Surface non-diatom log10(phytoplankton carbon)' 576 clunits = 'log10(mmolC/m3)' 577 clgrid = 'T' 578 579 CASE('SSPM') 580 581 clfiletype = 'sspmfb' 582 cllongname = 'Surface suspended particulate matter' 583 clunits = 'g/m3' 584 clgrid = 'T' 585 586 CASE('SFCO2','FCO2','fCO2','fco2') 587 588 clfiletype = 'sfco2fb' 589 cllongname = 'Surface fugacity of carbon dioxide' 590 clunits = 'uatm' 591 clgrid = 'T' 592 593 CASE('SPCO2') 594 595 clfiletype = 'spco2fb' 596 cllongname = 'Surface partial pressure of carbon dioxide' 597 clunits = 'uatm' 598 clgrid = 'T' 599 600 CASE DEFAULT 601 602 CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' ) 603 604 END SELECT 605 606 ! SLA needs special treatment because of MDT, so is done above 607 ! Remaining variables treated more generically 608 609 IF ( TRIM(surfdata%cvars(1)) /= 'SLA' ) THEN 610 376 611 CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 377 & 1 + iadd, iext, .TRUE. ) 378 379 clfiletype = 'sstfb' 612 & 1 + iadd + iadd_std, iext, .TRUE. ) 613 380 614 fbdata%cname(1) = surfdata%cvars(1) 381 fbdata%coblong(1) = 'Sea surface temperature'382 fbdata%cobunit(1) = 'Degree centigrade'615 fbdata%coblong(1) = cllongname 616 fbdata%cobunit(1) = clunits 383 617 DO je = 1, iext 384 618 fbdata%cextname(je) = pext%cdname(je) … … 386 620 fbdata%cextunit(je) = pext%cdunit(je,1) 387 621 END DO 388 fbdata%caddlong(1,1) = 'Model interpolated SST' 389 fbdata%caddunit(1,1) = 'Degree centigrade' 390 fbdata%cgrid(1) = 'T' 622 IF ( TRIM(surfdata%cvars(1)) == 'ICECONC' ) THEN 623 fbdata%caddlong(1,1) = 'Model interpolated ICE' 624 ELSE 625 fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(surfdata%cvars(1)) 626 ENDIF 627 fbdata%caddunit(1,1) = clunits 628 fbdata%cgrid(1) = clgrid 391 629 DO ja = 1, iadd 392 fbdata%caddname(1+ja) = padd%cdname(ja) 393 fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 394 fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 395 END DO 396 397 CASE('ICECON') 398 399 CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 400 & 1 + iadd, iext, .TRUE. ) 401 402 clfiletype = 'sicfb' 403 fbdata%cname(1) = surfdata%cvars(1) 404 fbdata%coblong(1) = 'Sea ice' 405 fbdata%cobunit(1) = 'Fraction' 406 DO je = 1, iext 407 fbdata%cextname(je) = pext%cdname(je) 408 fbdata%cextlong(je) = pext%cdlong(je,1) 409 fbdata%cextunit(je) = pext%cdunit(je,1) 410 END DO 411 fbdata%caddlong(1,1) = 'Model interpolated ICE' 412 fbdata%caddunit(1,1) = 'Fraction' 413 fbdata%cgrid(1) = 'T' 414 DO ja = 1, iadd 415 fbdata%caddname(1+ja) = padd%cdname(ja) 416 fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 417 fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 418 END DO 419 420 END SELECT 421 630 fbdata%caddname(1+iadd_std+ja) = padd%cdname(ja) 631 fbdata%caddlong(1+iadd_std+ja,1) = padd%cdlong(ja,1) 632 fbdata%caddunit(1+iadd_std+ja,1) = padd%cdunit(ja,1) 633 END DO 634 635 ENDIF 636 422 637 fbdata%caddname(1) = 'Hx' 423 638 IF ( indx_std /= -1 ) THEN 639 IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) iadd_std = iadd_std + 1 640 fbdata%caddname(1+iadd_std) = surfdata%cext(indx_std) 641 fbdata%caddlong(1+iadd_std,1) = 'Obs error standard deviation' 642 fbdata%caddunit(1+iadd_std,1) = fbdata%cobunit(1) 643 ENDIF 644 424 645 WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc 425 646 … … 439 660 fbdata%ivqc(jo,:) = 0 440 661 fbdata%ivqcf(:,jo,:) = 0 441 IF ( surfdata%nqc(jo) > 10) THEN662 IF ( surfdata%nqc(jo) > 255 ) THEN 442 663 fbdata%ioqc(jo) = 4 443 664 fbdata%ioqcf(1,jo) = 0 444 fbdata%ioqcf(2,jo) = surfdata%nqc(jo) - 10665 fbdata%ioqcf(2,jo) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111') 445 666 ELSE 446 667 fbdata%ioqc(jo) = surfdata%nqc(jo) … … 474 695 fbdata%idqc(1,jo) = 0 475 696 fbdata%idqcf(:,1,jo) = 0 476 IF ( surfdata%nqc(jo) > 10) THEN697 IF ( surfdata%nqc(jo) > 255 ) THEN 477 698 fbdata%ivqc(jo,1) = 4 478 699 fbdata%ivlqc(1,jo,1) = 4 479 700 fbdata%ivlqcf(1,1,jo,1) = 0 480 fbdata%ivlqcf(2,1,jo,1) = surfdata%nqc(jo) - 10701 fbdata%ivlqcf(2,1,jo,1) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111') 481 702 ELSE 482 703 fbdata%ivqc(jo,1) = surfdata%nqc(jo) … … 486 707 fbdata%iobsk(1,jo,1) = 0 487 708 IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2) 709 IF ( indx_std /= -1 ) THEN 710 fbdata%padd(1,jo,1+iadd_std,1) = surfdata%rext(jo,indx_std) 711 ENDIF 712 488 713 DO ja = 1, iadd 489 fbdata%padd(1,jo,2+ ja,1) = &714 fbdata%padd(1,jo,2+iadd_std+ja,1) = & 490 715 & surfdata%rext(jo,padd%ipoint(ja)) 491 716 END DO -
NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_h2d.h90
r11350 r11361 1240 1240 & zdum, & 1241 1241 & zaamax 1242 1242 1243 imax = -1 1243 1244 ! Main computation 1244 1245 pflt = 1.0_wp
Note: See TracChangeset
for help on using the changeset viewer.
