- Timestamp:
- 2015-12-04T11:56:46+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
r5997 r5998 6 6 !!====================================================================== 7 7 8 !!----------------------------------------------------------------------9 !! 'key_diaobs' : Switch on the observation diagnostic computation10 8 !!---------------------------------------------------------------------- 11 9 !! dia_obs_init : Reading and prepare observations … … 16 14 !! fin_date : Compute the final date YYYYMMDD.HHMMSS 17 15 !!---------------------------------------------------------------------- 18 !! * Modules used 16 !! * Modules used 19 17 USE wrk_nemo ! Memory Allocation 20 18 USE par_kind ! Precision variables … … 22 20 USE par_oce 23 21 USE dom_oce ! Ocean space and time domain variables 24 USE obs_fbm, ONLY: ln_cl4 ! Class 4 diagnostic switch 25 USE obs_read_prof ! Reading and allocation of observations (Coriolis) 26 USE obs_read_sla ! Reading and allocation of SLA observations 27 USE obs_read_sst ! Reading and allocation of SST observations 22 USE obs_read_prof ! Reading and allocation of profile obs 23 USE obs_read_surf ! Reading and allocation of surface obs 28 24 USE obs_sstbias ! Bias correction routine for SST 29 25 USE obs_readmdt ! Reading and allocation of MDT for SLA. 30 USE obs_read_seaice ! Reading and allocation of Sea Ice observations31 USE obs_read_vel ! Reading and allocation of velocity component observations32 26 USE obs_prep ! Preparation of obs. (grid search etc). 33 27 USE obs_oper ! Observation operators … … 36 30 USE obs_read_altbias ! Bias treatment for altimeter 37 31 USE obs_profiles_def ! Profile data definitions 38 USE obs_profiles ! Profile data storage39 32 USE obs_surf_def ! Surface data definitions 40 USE obs_sla ! SLA data storage41 USE obs_sst ! SST data storage42 USE obs_seaice ! Sea Ice data storage43 33 USE obs_types ! Definitions for observation types 44 34 USE mpp_map ! MPP mapping … … 55 45 & calc_date ! Compute the date of a timestep 56 46 57 !! * Shared Module variables58 LOGICAL, PUBLIC, PARAMETER :: &59 #if defined key_diaobs60 & lk_diaobs = .TRUE. !: Logical switch for observation diangostics61 #else62 & lk_diaobs = .FALSE. !: Logical switch for observation diangostics63 #endif64 65 47 !! * Module variables 66 LOGICAL, PUBLIC :: ln_t3d !: Logical switch for temperature profiles 67 LOGICAL, PUBLIC :: ln_s3d !: Logical switch for salinity profiles 68 LOGICAL, PUBLIC :: ln_ena !: Logical switch for the ENACT data set 69 LOGICAL, PUBLIC :: ln_cor !: Logical switch for the Coriolis data set 70 LOGICAL, PUBLIC :: ln_profb !: Logical switch for profile feedback datafiles 71 LOGICAL, PUBLIC :: ln_sla !: Logical switch for sea level anomalies 72 LOGICAL, PUBLIC :: ln_sladt !: Logical switch for SLA from AVISO files 73 LOGICAL, PUBLIC :: ln_slafb !: Logical switch for SLA from feedback files 74 LOGICAL, PUBLIC :: ln_sst !: Logical switch for sea surface temperature 75 LOGICAL, PUBLIC :: ln_reysst !: Logical switch for Reynolds sea surface temperature 76 LOGICAL, PUBLIC :: ln_ghrsst !: Logical switch for GHRSST data 77 LOGICAL, PUBLIC :: ln_sstfb !: Logical switch for SST from feedback files 78 LOGICAL, PUBLIC :: ln_seaice !: Logical switch for sea ice concentration 79 LOGICAL, PUBLIC :: ln_vel3d !: Logical switch for velocity component (u,v) observations 80 LOGICAL, PUBLIC :: ln_velavcur !: Logical switch for raw daily averaged netCDF current meter vel. data 81 LOGICAL, PUBLIC :: ln_velhrcur !: Logical switch for raw high freq netCDF current meter vel. data 82 LOGICAL, PUBLIC :: ln_velavadcp !: Logical switch for raw daily averaged netCDF ADCP vel. data 83 LOGICAL, PUBLIC :: ln_velhradcp !: Logical switch for raw high freq netCDF ADCP vel. data 84 LOGICAL, PUBLIC :: ln_velfb !: Logical switch for velocities from feedback files 85 LOGICAL, PUBLIC :: ln_ssh !: Logical switch for sea surface height 86 LOGICAL, PUBLIC :: ln_sss !: Logical switch for sea surface salinity 87 LOGICAL, PUBLIC :: ln_sstnight !: Logical switch for night mean SST observations 88 LOGICAL, PUBLIC :: ln_nea !: Remove observations near land 89 LOGICAL, PUBLIC :: ln_altbias !: Logical switch for altimeter bias 90 LOGICAL, PUBLIC :: ln_ignmis !: Logical switch for ignoring missing files 91 LOGICAL, PUBLIC :: ln_s_at_t !: Logical switch to compute model S at T observations 92 LOGICAL, PUBLIC :: ln_sstbias !: Logical switch for bias corection of SST 48 LOGICAL, PUBLIC :: ln_diaobs !: Logical switch for the obs operator 49 LOGICAL :: ln_sstnight !: Logical switch for night mean SST obs 93 50 94 REAL(KIND=dp), PUBLIC :: dobsini !: Observation window start date YYYYMMDD.HHMMSS 95 REAL(KIND=dp), PUBLIC :: dobsend !: Observation window end date YYYYMMDD.HHMMSS 96 97 INTEGER, PUBLIC :: n1dint !: Vertical interpolation method 98 INTEGER, PUBLIC :: n2dint !: Horizontal interpolation method 99 51 INTEGER :: nn_1dint !: Vertical interpolation method 52 INTEGER :: nn_2dint !: Horizontal interpolation method 100 53 INTEGER, DIMENSION(imaxavtypes) :: & 101 & endailyavtypes !: ENACT data types which are daily average 102 103 INTEGER, PARAMETER :: MaxNumFiles = 1000 104 LOGICAL, DIMENSION(MaxNumFiles) :: & 105 & ln_profb_ena, & !: Is the feedback files from ENACT data ? 106 ! !: If so use endailyavtypes 107 & ln_profb_enatim !: Change tim for 820 enact data set. 108 109 INTEGER, DIMENSION(MaxNumFiles), PUBLIC :: sstbias_type !SST bias type 110 111 LOGICAL, DIMENSION(MaxNumFiles) :: & 112 & ln_velfb_av !: Is the velocity feedback files daily average? 113 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 114 & ld_enact !: Profile data is ENACT so use endailyavtypes 115 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 116 & ld_velav !: Velocity data is daily averaged 117 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 118 & ld_sstnight !: SST observation corresponds to night mean 54 & nn_profdavtypes !: Profile data types representing a daily average 55 INTEGER :: nproftypes !: Number of profile obs types 56 INTEGER :: nsurftypes !: Number of surface obs types 57 INTEGER, DIMENSION(:), ALLOCATABLE :: & 58 & nvarsprof, & !: Number of profile variables 59 & nvarssurf !: Number of surface variables 60 INTEGER, DIMENSION(:), ALLOCATABLE :: & 61 & nextrprof, & !: Number of profile extra variables 62 & nextrsurf !: Number of surface extra variables 63 INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:),:: sstbias_type !SST bias type 64 TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: & 65 & surfdata, & !: Initial surface data 66 & surfdataqc !: Surface data after quality control 67 TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: & 68 & profdata, & !: Initial profile data 69 & profdataqc !: Profile data after quality control 70 71 CHARACTER(len=6), PUBLIC, DIMENSION(:), ALLOCATABLE :: & 72 & cobstypesprof, & !: Profile obs types 73 & cobstypessurf !: Surface obs types 119 74 120 75 !!---------------------------------------------------------------------- … … 144 99 !! ! 07-03 (K. Mogensen) General handling of profiles 145 100 !! ! 14-08 (J.While) Incorporated SST bias correction 101 !! ! 15-02 (M. Martin) Simplification of namelist and code 146 102 !!---------------------------------------------------------------------- 147 103 … … 149 105 150 106 !! * Local declarations 151 CHARACTER(len=128) :: enactfiles(MaxNumFiles) 152 CHARACTER(len=128) :: coriofiles(MaxNumFiles) 153 CHARACTER(len=128) :: profbfiles(MaxNumFiles) 154 CHARACTER(len=128) :: sstfiles(MaxNumFiles) 155 CHARACTER(len=128) :: sstfbfiles(MaxNumFiles) 156 CHARACTER(len=128) :: sstbias_files(MaxNumFiles) 157 CHARACTER(len=128) :: slafilesact(MaxNumFiles) 158 CHARACTER(len=128) :: slafilespas(MaxNumFiles) 159 CHARACTER(len=128) :: slafbfiles(MaxNumFiles) 160 CHARACTER(len=128) :: seaicefiles(MaxNumFiles) 161 CHARACTER(len=128) :: velcurfiles(MaxNumFiles) 162 CHARACTER(len=128) :: veladcpfiles(MaxNumFiles) 163 CHARACTER(len=128) :: velavcurfiles(MaxNumFiles) 164 CHARACTER(len=128) :: velhrcurfiles(MaxNumFiles) 165 CHARACTER(len=128) :: velavadcpfiles(MaxNumFiles) 166 CHARACTER(len=128) :: velhradcpfiles(MaxNumFiles) 167 CHARACTER(len=128) :: velfbfiles(MaxNumFiles) 168 CHARACTER(LEN=128) :: reysstname 169 CHARACTER(LEN=12) :: reysstfmt 170 CHARACTER(LEN=128) :: bias_file 171 CHARACTER(LEN=20) :: datestr=" ", timestr=" " 172 NAMELIST/namobs/ln_ena, ln_cor, ln_profb, ln_t3d, ln_s3d, & 173 & ln_sla, ln_sladt, ln_slafb, & 174 & ln_ssh, ln_sst, ln_sstfb, ln_sss, ln_nea, & 175 & enactfiles, coriofiles, profbfiles, & 176 & slafilesact, slafilespas, slafbfiles, & 177 & sstfiles, sstfbfiles, & 178 & ln_seaice, seaicefiles, & 179 & dobsini, dobsend, n1dint, n2dint, & 180 & nmsshc, mdtcorr, mdtcutoff, & 181 & ln_reysst, ln_ghrsst, reysstname, reysstfmt, & 182 & ln_sstnight, & 107 INTEGER, PARAMETER :: & 108 & jpmaxnfiles = 1000 ! Maximum number of files for each obs type 109 INTEGER, DIMENSION(:), ALLOCATABLE :: & 110 & ifilesprof, & ! Number of profile files 111 & ifilessurf ! Number of surface files 112 INTEGER :: ios ! Local integer output status for namelist read 113 INTEGER :: jtype ! Counter for obs types 114 INTEGER :: jvar ! Counter for variables 115 INTEGER :: jfile ! Counter for files 116 117 CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: & 118 & cn_profbfiles, & ! T/S profile input filenames 119 & cn_sstfbfiles, & ! Sea surface temperature input filenames 120 & cn_slafbfiles, & ! Sea level anomaly input filenames 121 & cn_sicfbfiles, & ! Seaice concentration input filenames 122 & cn_velfbfiles ! Velocity profile input filenames 123 CHARACTER(LEN=128) :: & 124 & cn_altbiasfile ! Altimeter bias input filename 125 & cn_sstbias_files ! Altimeter bias input filenames 126 CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: & 127 & clproffiles, & ! Profile filenames 128 & clsurffiles ! Surface filenames 129 130 LOGICAL :: ln_t3d ! Logical switch for temperature profiles 131 LOGICAL :: ln_s3d ! Logical switch for salinity profiles 132 LOGICAL :: ln_sla ! Logical switch for sea level anomalies 133 LOGICAL :: ln_sst ! Logical switch for sea surface temperature 134 LOGICAL :: ln_sic ! Logical switch for sea ice concentration 135 LOGICAL :: ln_vel3d ! Logical switch for velocity (u,v) obs 136 LOGICAL :: ln_nea ! Logical switch to remove obs near land 137 LOGICAL :: ln_altbias ! Logical switch for altimeter bias 138 LOGICAL :: ln_sstbias !: Logical switch for bias corection of SST 139 LOGICAL :: ln_ignmis ! Logical switch for ignoring missing files 140 LOGICAL :: ln_s_at_t ! Logical switch to compute model S at T obs 141 LOGICAL :: llvar1 ! Logical for profile variable 1 142 LOGICAL :: llvar2 ! Logical for profile variable 1 143 LOGICAL :: llnightav ! Logical for calculating night-time averages 144 145 REAL(dp) :: rn_dobsini ! Obs window start date YYYYMMDD.HHMMSS 146 REAL(dp) :: rn_dobsend ! Obs window end date YYYYMMDD.HHMMSS 147 REAL(wp), POINTER, DIMENSION(:,:) :: & 148 & zglam1, & ! Model longitudes for profile variable 1 149 & zglam2 ! Model longitudes for profile variable 2 150 REAL(wp), POINTER, DIMENSION(:,:) :: & 151 & zgphi1, & ! Model latitudes for profile variable 1 152 & zgphi2 ! Model latitudes for profile variable 2 153 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 154 & zmask1, & ! Model land/sea mask associated with variable 1 155 & zmask2 ! Model land/sea mask associated with variable 2 156 157 NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla, & 158 & ln_sst, ln_sic, ln_vel3d, & 159 & ln_altbias, ln_nea, ln_grid_global, & 183 160 & ln_grid_search_lookup, & 184 & grid_search_file, grid_search_res, & 185 & ln_grid_global, bias_file, ln_altbias, & 186 & endailyavtypes, ln_s_at_t, ln_profb_ena, & 187 & ln_vel3d, ln_velavcur, velavcurfiles, & 188 & ln_velhrcur, velhrcurfiles, & 189 & ln_velavadcp, velavadcpfiles, & 190 & ln_velhradcp, velhradcpfiles, & 191 & ln_velfb, velfbfiles, ln_velfb_av, & 192 & ln_profb_enatim, ln_ignmis, ln_cl4, & 193 & ln_sstbias, sstbias_files 194 195 INTEGER :: jprofset 196 INTEGER :: jveloset 197 INTEGER :: jvar 198 INTEGER :: jnumenact 199 INTEGER :: jnumcorio 200 INTEGER :: jnumprofb 201 INTEGER :: jnumslaact 202 INTEGER :: jnumslapas 203 INTEGER :: jnumslafb 204 INTEGER :: jnumsst 205 INTEGER :: jnumsstfb 206 INTEGER :: jnumsstbias 207 INTEGER :: jnumseaice 208 INTEGER :: jnumvelavcur 209 INTEGER :: jnumvelhrcur 210 INTEGER :: jnumvelavadcp 211 INTEGER :: jnumvelhradcp 212 INTEGER :: jnumvelfb 213 INTEGER :: ji 214 INTEGER :: jset 215 INTEGER :: ios ! Local integer output status for namelist read 216 LOGICAL :: lmask(MaxNumFiles), ll_u3d, ll_v3d 161 & ln_ignmis, ln_s_at_t, ln_sstnight, & 162 & cn_profbfiles, cn_slafbfiles, & 163 & cn_sstfbfiles, cn_sicfbfiles, & 164 & cn_velfbfiles, cn_altbiasfile, & 165 & cn_gridsearchfile, rn_gridsearchres, & 166 & rn_dobsini, rn_dobsend, nn_1dint, nn_2dint, & 167 & nn_msshc, rn_mdtcorr, rn_mdtcutoff, & 168 & nn_profdavtypes, ln_sstbias, sstbias_files 169 170 INTEGER :: jnumsstbias !TG - Is this still needed 171 CALL wrk_alloc( jpi, jpj, zglam1 ) 172 CALL wrk_alloc( jpi, jpj, zglam2 ) 173 CALL wrk_alloc( jpi, jpj, zgphi1 ) 174 CALL wrk_alloc( jpi, jpj, zgphi2 ) 175 CALL wrk_alloc( jpi, jpj, jpk, zmask1 ) 176 CALL wrk_alloc( jpi, jpj, jpk, zmask2 ) 217 177 218 178 !----------------------------------------------------------------------- … … 221 181 222 182 !Initalise all values in namelist arrays 223 enactfiles(:) = '' 224 coriofiles(:) = '' 225 profbfiles(:) = '' 226 slafilesact(:) = '' 227 slafilespas(:) = '' 228 slafbfiles(:) = '' 229 sstfiles(:) = '' 230 sstfbfiles(:) = '' 231 seaicefiles(:) = '' 232 velcurfiles(:) = '' 233 veladcpfiles(:) = '' 234 velavcurfiles(:) = '' 235 velhrcurfiles(:) = '' 236 velavadcpfiles(:) = '' 237 velhradcpfiles(:) = '' 238 velfbfiles(:) = '' 239 velcurfiles(:) = '' 240 veladcpfiles(:) = '' 241 endailyavtypes(:) = -1 242 endailyavtypes(1) = 820 243 ln_profb_ena(:) = .FALSE. 244 ln_profb_enatim(:) = .TRUE. 245 ln_velfb_av(:) = .FALSE. 246 ln_ignmis = .FALSE. 247 248 CALL ini_date( dobsini ) 249 CALL fin_date( dobsend ) 250 251 ! Read Namelist namobs : control observation diagnostics 252 REWIND( numnam_ref ) ! Namelist namobs in reference namelist : Diagnostic: control observation 183 ALLOCATE(sstbias_type(jpmaxnumfiles)) 184 ! Some namelist arrays need initialising 185 cn_profbfiles(:) = '' 186 cn_slafbfiles(:) = '' 187 cn_sstfbfiles(:) = '' 188 cn_sicfbfiles(:) = '' 189 cn_velfbfiles(:) = '' 190 cn_sstbias_files(:) = '' 191 nn_profdavtypes(:) = -1 192 193 CALL ini_date( rn_dobsini ) 194 CALL fin_date( rn_dobsend ) 195 196 ! Read namelist namobs : control observation diagnostics 197 REWIND( numnam_ref ) ! Namelist namobs in reference namelist 253 198 READ ( numnam_ref, namobs, IOSTAT = ios, ERR = 901) 254 199 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp ) 255 200 256 REWIND( numnam_cfg ) ! Namelist namobs in configuration namelist : Diagnostic: control observation201 REWIND( numnam_cfg ) ! Namelist namobs in configuration namelist 257 202 READ ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 ) 258 203 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp ) 259 204 IF(lwm) WRITE ( numond, namobs ) 260 205 261 ! Count number of files for each type 262 IF (ln_ena) THEN 263 lmask(:) = .FALSE. 264 WHERE (enactfiles(:) /= '') lmask(:) = .TRUE. 265 jnumenact = COUNT(lmask) 266 ENDIF 267 IF (ln_cor) THEN 268 lmask(:) = .FALSE. 269 WHERE (coriofiles(:) /= '') lmask(:) = .TRUE. 270 jnumcorio = COUNT(lmask) 271 ENDIF 272 IF (ln_profb) THEN 273 lmask(:) = .FALSE. 274 WHERE (profbfiles(:) /= '') lmask(:) = .TRUE. 275 jnumprofb = COUNT(lmask) 276 ENDIF 277 IF (ln_sladt) THEN 278 lmask(:) = .FALSE. 279 WHERE (slafilesact(:) /= '') lmask(:) = .TRUE. 280 jnumslaact = COUNT(lmask) 281 lmask(:) = .FALSE. 282 WHERE (slafilespas(:) /= '') lmask(:) = .TRUE. 283 jnumslapas = COUNT(lmask) 284 ENDIF 285 IF (ln_slafb) THEN 286 lmask(:) = .FALSE. 287 WHERE (slafbfiles(:) /= '') lmask(:) = .TRUE. 288 jnumslafb = COUNT(lmask) 289 lmask(:) = .FALSE. 290 ENDIF 291 IF (ln_ghrsst) THEN 292 lmask(:) = .FALSE. 293 WHERE (sstfiles(:) /= '') lmask(:) = .TRUE. 294 jnumsst = COUNT(lmask) 295 ENDIF 296 IF (ln_sstfb) THEN 297 lmask(:) = .FALSE. 298 WHERE (sstfbfiles(:) /= '') lmask(:) = .TRUE. 299 jnumsstfb = COUNT(lmask) 300 lmask(:) = .FALSE. 301 ENDIF 206 IF ( .NOT. ln_diaobs ) THEN 207 IF(lwp) WRITE(numout,cform_war) 208 IF(lwp) WRITE(numout,*)' ln_diaobs is set to false so not calling dia_obs' 209 RETURN 210 ENDIF 211 212 !----------------------------------------------------------------------- 213 ! Set up list of observation types to be used 214 ! and the files associated with each type 215 !----------------------------------------------------------------------- 216 217 nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 218 nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic /) ) 219 302 220 IF (ln_sstbias) THEN 303 221 lmask(:) = .FALSE. … … 306 224 lmask(:) = .FALSE. 307 225 ENDIF 308 IF (ln_seaice) THEN 309 lmask(:) = .FALSE. 310 WHERE (seaicefiles(:) /= '') lmask(:) = .TRUE. 311 jnumseaice = COUNT(lmask) 312 ENDIF 313 IF (ln_velavcur) THEN 314 lmask(:) = .FALSE. 315 WHERE (velavcurfiles(:) /= '') lmask(:) = .TRUE. 316 jnumvelavcur = COUNT(lmask) 317 ENDIF 318 IF (ln_velhrcur) THEN 319 lmask(:) = .FALSE. 320 WHERE (velhrcurfiles(:) /= '') lmask(:) = .TRUE. 321 jnumvelhrcur = COUNT(lmask) 322 ENDIF 323 IF (ln_velavadcp) THEN 324 lmask(:) = .FALSE. 325 WHERE (velavadcpfiles(:) /= '') lmask(:) = .TRUE. 326 jnumvelavadcp = COUNT(lmask) 327 ENDIF 328 IF (ln_velhradcp) THEN 329 lmask(:) = .FALSE. 330 WHERE (velhradcpfiles(:) /= '') lmask(:) = .TRUE. 331 jnumvelhradcp = COUNT(lmask) 332 ENDIF 333 IF (ln_velfb) THEN 334 lmask(:) = .FALSE. 335 WHERE (velfbfiles(:) /= '') lmask(:) = .TRUE. 336 jnumvelfb = COUNT(lmask) 337 lmask(:) = .FALSE. 338 ENDIF 339 340 ! Control print 226 227 IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 228 IF(lwp) WRITE(numout,cform_war) 229 IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 230 & ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', & 231 & ' are set to .FALSE. so turning off calls to dia_obs' 232 nwarn = nwarn + 1 233 ln_diaobs = .FALSE. 234 RETURN 235 ENDIF 236 237 IF ( nproftypes > 0 ) THEN 238 239 ALLOCATE( cobstypesprof(nproftypes) ) 240 ALLOCATE( ifilesprof(nproftypes) ) 241 ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 242 243 jtype = 0 244 IF (ln_t3d .OR. ln_s3d) THEN 245 jtype = jtype + 1 246 clproffiles(jtype,:) = cn_profbfiles(:) 247 cobstypesprof(jtype) = 'prof ' 248 ifilesprof(jtype) = 0 249 DO jfile = 1, jpmaxnfiles 250 IF ( trim(clproffiles(jtype,jfile)) /= '' ) & 251 ifilesprof(jtype) = ifilesprof(jtype) + 1 252 END DO 253 ENDIF 254 IF (ln_vel3d) THEN 255 jtype = jtype + 1 256 clproffiles(jtype,:) = cn_velfbfiles(:) 257 cobstypesprof(jtype) = 'vel ' 258 ifilesprof(jtype) = 0 259 DO jfile = 1, jpmaxnfiles 260 IF ( trim(clproffiles(jtype,jfile)) /= '' ) & 261 ifilesprof(jtype) = ifilesprof(jtype) + 1 262 END DO 263 ENDIF 264 265 ENDIF 266 267 IF ( nsurftypes > 0 ) THEN 268 269 ALLOCATE( cobstypessurf(nsurftypes) ) 270 ALLOCATE( ifilessurf(nsurftypes) ) 271 ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 272 273 jtype = 0 274 IF (ln_sla) THEN 275 jtype = jtype + 1 276 clsurffiles(jtype,:) = cn_slafbfiles(:) 277 cobstypessurf(jtype) = 'sla ' 278 ifilessurf(jtype) = 0 279 DO jfile = 1, jpmaxnfiles 280 IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 281 ifilessurf(jtype) = ifilessurf(jtype) + 1 282 END DO 283 ENDIF 284 IF (ln_sst) THEN 285 jtype = jtype + 1 286 clsurffiles(jtype,:) = cn_sstfbfiles(:) 287 cobstypessurf(jtype) = 'sst ' 288 ifilessurf(jtype) = 0 289 DO jfile = 1, jpmaxnfiles 290 IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 291 ifilessurf(jtype) = ifilessurf(jtype) + 1 292 END DO 293 ENDIF 294 #if defined key_lim2 || defined key_lim3 295 IF (ln_sic) THEN 296 jtype = jtype + 1 297 clsurffiles(jtype,:) = cn_sicfbfiles(:) 298 cobstypessurf(jtype) = 'sic ' 299 ifilessurf(jtype) = 0 300 DO jfile = 1, jpmaxnfiles 301 IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 302 ifilessurf(jtype) = ifilessurf(jtype) + 1 303 END DO 304 ENDIF 305 #endif 306 307 ENDIF 308 309 !Write namelist settings to stdout 341 310 IF(lwp) THEN 342 311 WRITE(numout,*) … … 344 313 WRITE(numout,*) '~~~~~~~~~~~~' 345 314 WRITE(numout,*) ' Namelist namobs : set observation diagnostic parameters' 346 WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d 347 WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d 348 WRITE(numout,*) ' Logical switch for ENACT insitu data set ln_ena = ', ln_ena 349 WRITE(numout,*) ' Logical switch for Coriolis insitu data set ln_cor = ', ln_cor 350 WRITE(numout,*) ' Logical switch for feedback insitu data set ln_profb = ', ln_profb 351 WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla 352 WRITE(numout,*) ' Logical switch for AVISO SLA data ln_sladt = ', ln_sladt 353 WRITE(numout,*) ' Logical switch for feedback SLA data ln_slafb = ', ln_slafb 354 WRITE(numout,*) ' Logical switch for SSH observations ln_ssh = ', ln_ssh 355 WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst 356 WRITE(numout,*) ' Logical switch for Reynolds observations ln_reysst = ', ln_reysst 357 WRITE(numout,*) ' Logical switch for GHRSST observations ln_ghrsst = ', ln_ghrsst 358 WRITE(numout,*) ' Logical switch for feedback SST data ln_sstfb = ', ln_sstfb 315 WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d 316 WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d 317 WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla 318 WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst 319 WRITE(numout,*) ' Logical switch for Sea Ice observations ln_sic = ', ln_sic 320 WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d 321 WRITE(numout,*) ' Global distribution of observations ln_grid_global = ',ln_grid_global 359 322 WRITE(numout,*) ' Logical switch for SST bias correction ln_sstbias = ', ln_sstbias 360 WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight 361 WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss 362 WRITE(numout,*) ' Logical switch for Sea Ice observations ln_seaice = ', ln_seaice 363 WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d 364 WRITE(numout,*) ' Logical switch for velocity daily av. cur. ln_velavcur = ', ln_velavcur 365 WRITE(numout,*) ' Logical switch for velocity high freq. cur. ln_velhrcur = ', ln_velhrcur 366 WRITE(numout,*) ' Logical switch for velocity daily av. ADCP ln_velavadcp = ', ln_velavadcp 367 WRITE(numout,*) ' Logical switch for velocity high freq. ADCP ln_velhradcp = ', ln_velhradcp 368 WRITE(numout,*) ' Logical switch for feedback velocity data ln_velfb = ', ln_velfb 369 WRITE(numout,*) ' Global distribtion of observations ln_grid_global = ',ln_grid_global 370 WRITE(numout,*) & 371 ' Logical switch for obs grid search w/lookup table ln_grid_search_lookup = ',ln_grid_search_lookup 323 WRITE(numout,*) ' Logical switch for obs grid search lookup ln_grid_search_lookup = ',ln_grid_search_lookup 372 324 IF (ln_grid_search_lookup) & 373 WRITE(numout,*) ' Grid search lookup file header grid_search_file = ', grid_search_file 374 IF (ln_ena) THEN 375 DO ji = 1, jnumenact 376 WRITE(numout,'(1X,2A)') ' ENACT input observation file name enactfiles = ', & 377 TRIM(enactfiles(ji)) 325 WRITE(numout,*) ' Grid search lookup file header cn_gridsearchfile = ', cn_gridsearchfile 326 WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS rn_dobsini = ', rn_dobsini 327 WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS rn_dobsend = ', rn_dobsend 328 WRITE(numout,*) ' Type of vertical interpolation method nn_1dint = ', nn_1dint 329 WRITE(numout,*) ' Type of horizontal interpolation method nn_2dint = ', nn_2dint 330 WRITE(numout,*) ' Rejection of observations near land switch ln_nea = ', ln_nea 331 WRITE(numout,*) ' MSSH correction scheme nn_msshc = ', nn_msshc 332 WRITE(numout,*) ' MDT correction rn_mdtcorr = ', rn_mdtcorr 333 WRITE(numout,*) ' MDT cutoff for computed correction rn_mdtcutoff = ', rn_mdtcutoff 334 WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias 335 WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis 336 WRITE(numout,*) ' Daily average types nn_profdavtypes = ', nn_profdavtypes 337 WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight 338 WRITE(numout,*) ' Number of profile obs types: ',nproftypes 339 340 IF ( nproftypes > 0 ) THEN 341 DO jtype = 1, nproftypes 342 DO jfile = 1, ifilesprof(jtype) 343 WRITE(numout,'(1X,2A)') ' '//cobstypesprof(jtype)//' input observation file names = ', & 344 TRIM(clproffiles(jtype,jfile)) 345 END DO 378 346 END DO 379 347 ENDIF 380 IF (ln_cor) THEN 381 DO ji = 1, jnumcorio 382 WRITE(numout,'(1X,2A)') ' Coriolis input observation file name coriofiles = ', & 383 TRIM(coriofiles(ji)) 348 349 WRITE(numout,*)' Number of surface obs types: ',nsurftypes 350 IF ( nsurftypes > 0 ) THEN 351 DO jtype = 1, nsurftypes 352 DO jfile = 1, ifilessurf(jtype) 353 WRITE(numout,'(1X,2A)') ' '//cobstypessurf(jtype)//' input observation file names = ', & 354 TRIM(clsurffiles(jtype,jfile)) 355 END DO 384 356 END DO 385 357 ENDIF 386 IF (ln_profb) THEN 387 DO ji = 1, jnumprofb 388 IF (ln_profb_ena(ji)) THEN 389 WRITE(numout,'(1X,2A)') ' Enact feedback input observation file name profbfiles = ', & 390 TRIM(profbfiles(ji)) 391 ELSE 392 WRITE(numout,'(1X,2A)') ' Feedback input observation file name profbfiles = ', & 393 TRIM(profbfiles(ji)) 394 ENDIF 395 WRITE(numout,'(1X,2A)') ' Enact feedback input time setting switch ln_profb_enatim = ', ln_profb_enatim(ji) 396 END DO 397 ENDIF 398 IF (ln_sladt) THEN 399 DO ji = 1, jnumslaact 400 WRITE(numout,'(1X,2A)') ' Active SLA input observation file name slafilesact = ', & 401 TRIM(slafilesact(ji)) 402 END DO 403 DO ji = 1, jnumslapas 404 WRITE(numout,'(1X,2A)') ' Passive SLA input observation file name slafilespas = ', & 405 TRIM(slafilespas(ji)) 406 END DO 407 ENDIF 408 IF (ln_slafb) THEN 409 DO ji = 1, jnumslafb 410 WRITE(numout,'(1X,2A)') ' Feedback SLA input observation file name slafbfiles = ', & 411 TRIM(slafbfiles(ji)) 412 END DO 413 ENDIF 414 IF (ln_ghrsst) THEN 415 DO ji = 1, jnumsst 416 WRITE(numout,'(1X,2A)') ' GHRSST input observation file name sstfiles = ', & 417 TRIM(sstfiles(ji)) 418 END DO 419 ENDIF 420 IF (ln_sstfb) THEN 421 DO ji = 1, jnumsstfb 422 WRITE(numout,'(1X,2A)') ' Feedback SST input observation file name sstfbfiles = ', & 423 TRIM(sstfbfiles(ji)) 424 END DO 425 ENDIF 426 IF (ln_seaice) THEN 427 DO ji = 1, jnumseaice 428 WRITE(numout,'(1X,2A)') ' Sea Ice input observation file name seaicefiles = ', & 429 TRIM(seaicefiles(ji)) 430 END DO 431 ENDIF 432 IF (ln_velavcur) THEN 433 DO ji = 1, jnumvelavcur 434 WRITE(numout,'(1X,2A)') ' Vel. cur. daily av. input file name velavcurfiles = ', & 435 TRIM(velavcurfiles(ji)) 436 END DO 437 ENDIF 438 IF (ln_velhrcur) THEN 439 DO ji = 1, jnumvelhrcur 440 WRITE(numout,'(1X,2A)') ' Vel. cur. high freq. input file name velhvcurfiles = ', & 441 TRIM(velhrcurfiles(ji)) 442 END DO 443 ENDIF 444 IF (ln_velavadcp) THEN 445 DO ji = 1, jnumvelavadcp 446 WRITE(numout,'(1X,2A)') ' Vel. ADCP daily av. input file name velavadcpfiles = ', & 447 TRIM(velavadcpfiles(ji)) 448 END DO 449 ENDIF 450 IF (ln_velhradcp) THEN 451 DO ji = 1, jnumvelhradcp 452 WRITE(numout,'(1X,2A)') ' Vel. ADCP high freq. input file name velhvadcpfiles = ', & 453 TRIM(velhradcpfiles(ji)) 454 END DO 455 ENDIF 456 IF (ln_velfb) THEN 457 DO ji = 1, jnumvelfb 458 IF (ln_velfb_av(ji)) THEN 459 WRITE(numout,'(1X,2A)') ' Vel. feedback daily av. input file name velfbfiles = ', & 460 TRIM(velfbfiles(ji)) 461 ELSE 462 WRITE(numout,'(1X,2A)') ' Vel. feedback input observation file name velfbfiles = ', & 463 TRIM(velfbfiles(ji)) 464 ENDIF 465 END DO 466 ENDIF 467 WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS dobsini = ', dobsini 468 WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS dobsend = ', dobsend 469 WRITE(numout,*) ' Type of vertical interpolation method n1dint = ', n1dint 470 WRITE(numout,*) ' Type of horizontal interpolation method n2dint = ', n2dint 471 WRITE(numout,*) ' Rejection of observations near land swithch ln_nea = ', ln_nea 472 WRITE(numout,*) ' MSSH correction scheme nmsshc = ', nmsshc 473 WRITE(numout,*) ' MDT correction mdtcorr = ', mdtcorr 474 WRITE(numout,*) ' MDT cutoff for computed correction mdtcutoff = ', mdtcutoff 475 WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias 476 WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis 477 WRITE(numout,*) ' ENACT daily average types = ',endailyavtypes 478 479 ENDIF 480 358 WRITE(numout,*) '~~~~~~~~~~~~' 359 360 ENDIF 361 362 !----------------------------------------------------------------------- 363 ! Obs operator parameter checking and initialisations 364 !----------------------------------------------------------------------- 365 481 366 IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN 482 367 CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' ) … … 484 369 ENDIF 485 370 486 CALL obs_typ_init 487 488 CALL mppmap_init 489 490 ! Parameter control 491 #if defined key_diaobs 492 IF ( ( .NOT. ln_t3d ).AND.( .NOT. ln_s3d ).AND.( .NOT. ln_sla ).AND. & 493 & ( .NOT. ln_vel3d ).AND. & 494 & ( .NOT. ln_ssh ).AND.( .NOT. ln_sst ).AND.( .NOT. ln_sss ).AND. & 495 & ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ) ) THEN 496 IF(lwp) WRITE(numout,cform_war) 497 IF(lwp) WRITE(numout,*) ' key_diaobs is activated but logical flags', & 498 & ' ln_t3d, ln_s3d, ln_sla, ln_ssh, ln_sst, ln_sss, ln_seaice, ln_vel3d are all set to .FALSE.' 499 nwarn = nwarn + 1 500 ENDIF 501 #endif 502 503 CALL obs_grid_setup( ) 504 IF ( ( n1dint < 0 ).OR.( n1dint > 1 ) ) THEN 371 IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN 505 372 CALL ctl_stop(' Choice of vertical (1D) interpolation method', & 506 373 & ' is not available') 507 374 ENDIF 508 IF ( ( n2dint < 0 ).OR.( n2dint > 4 ) ) THEN 375 376 IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 4 ) ) THEN 509 377 CALL ctl_stop(' Choice of horizontal (2D) interpolation method', & 510 378 & ' is not available') 511 379 ENDIF 512 380 381 CALL obs_typ_init 382 383 CALL mppmap_init 384 385 CALL obs_grid_setup( ) 386 513 387 !----------------------------------------------------------------------- 514 388 ! Depending on switches read the various observation types 515 389 !----------------------------------------------------------------------- 516 ! - Temperature/salinity profiles 517 518 IF ( ln_t3d .OR. ln_s3d ) THEN 519 520 ! Set the number of variables for profiles to 2 (T and S) 521 nprofvars = 2 522 ! Set the number of extra variables for profiles to 1 (insitu temp). 523 nprofextr = 1 524 525 ! Count how may insitu data sets we have and allocate data. 526 jprofset = 0 527 IF ( ln_ena ) jprofset = jprofset + 1 528 IF ( ln_cor ) jprofset = jprofset + 1 529 IF ( ln_profb ) jprofset = jprofset + jnumprofb 530 nprofsets = jprofset 531 IF ( nprofsets > 0 ) THEN 532 ALLOCATE(ld_enact(nprofsets)) 533 ALLOCATE(profdata(nprofsets)) 534 ALLOCATE(prodatqc(nprofsets)) 535 ENDIF 536 537 jprofset = 0 538 539 ! ENACT insitu data 540 541 IF ( ln_ena ) THEN 542 543 jprofset = jprofset + 1 544 545 ld_enact(jprofset) = .TRUE. 546 547 CALL obs_rea_pro_dri( 1, profdata(jprofset), & 548 & jnumenact, enactfiles(1:jnumenact), & 549 & nprofvars, nprofextr, & 550 & nitend-nit000+2, & 551 & dobsini, dobsend, ln_t3d, ln_s3d, & 552 & ln_ignmis, ln_s_at_t, .TRUE., .FALSE., & 553 & kdailyavtypes = endailyavtypes ) 554 555 DO jvar = 1, 2 556 557 CALL obs_prof_staend( profdata(jprofset), jvar ) 558 390 391 IF ( nproftypes > 0 ) THEN 392 393 ALLOCATE(profdata(nproftypes)) 394 ALLOCATE(profdataqc(nproftypes)) 395 ALLOCATE(nvarsprof(nproftypes)) 396 ALLOCATE(nextrprof(nproftypes)) 397 398 DO jtype = 1, nproftypes 399 400 nvarsprof(jtype) = 2 401 IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN 402 nextrprof(jtype) = 1 403 llvar1 = ln_t3d 404 llvar2 = ln_s3d 405 zglam1 = glamt 406 zgphi1 = gphit 407 zmask1 = tmask 408 zglam2 = glamt 409 zgphi2 = gphit 410 zmask2 = tmask 411 ENDIF 412 IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN 413 nextrprof(jtype) = 2 414 llvar1 = ln_vel3d 415 llvar2 = ln_vel3d 416 zglam1 = glamu 417 zgphi1 = gphiu 418 zmask1 = umask 419 zglam2 = glamv 420 zgphi2 = gphiv 421 zmask2 = vmask 422 ENDIF 423 424 !Read in profile or profile obs types 425 CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype), & 426 & clproffiles(jtype,1:ifilesprof(jtype)), & 427 & nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, & 428 & rn_dobsini, rn_dobsend, llvar1, llvar2, & 429 & ln_ignmis, ln_s_at_t, .FALSE., & 430 & kdailyavtypes = nn_profdavtypes ) 431 432 DO jvar = 1, nvarsprof(jtype) 433 CALL obs_prof_staend( profdata(jtype), jvar ) 559 434 END DO 560 435 561 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 562 & ln_t3d, ln_s3d, ln_nea, & 563 & kdailyavtypes=endailyavtypes ) 564 565 ENDIF 566 567 ! Coriolis insitu data 568 569 IF ( ln_cor ) THEN 570 571 jprofset = jprofset + 1 572 573 ld_enact(jprofset) = .FALSE. 574 575 CALL obs_rea_pro_dri( 2, profdata(jprofset), & 576 & jnumcorio, coriofiles(1:jnumcorio), & 577 & nprofvars, nprofextr, & 578 & nitend-nit000+2, & 579 & dobsini, dobsend, ln_t3d, ln_s3d, & 580 & ln_ignmis, ln_s_at_t, .FALSE., .FALSE. ) 581 582 DO jvar = 1, 2 583 584 CALL obs_prof_staend( profdata(jprofset), jvar ) 585 586 END DO 587 588 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 589 & ln_t3d, ln_s3d, ln_nea ) 590 591 ENDIF 592 593 ! Feedback insitu data 594 595 IF ( ln_profb ) THEN 596 597 DO jset = 1, jnumprofb 598 599 jprofset = jprofset + 1 600 ld_enact (jprofset) = ln_profb_ena(jset) 601 602 CALL obs_rea_pro_dri( 0, profdata(jprofset), & 603 & 1, profbfiles(jset:jset), & 604 & nprofvars, nprofextr, & 605 & nitend-nit000+2, & 606 & dobsini, dobsend, ln_t3d, ln_s3d, & 607 & ln_ignmis, ln_s_at_t, & 608 & ld_enact(jprofset).AND.& 609 & ln_profb_enatim(jset), & 610 & .FALSE., kdailyavtypes = endailyavtypes ) 611 612 DO jvar = 1, 2 613 614 CALL obs_prof_staend( profdata(jprofset), jvar ) 615 616 END DO 617 618 IF ( ld_enact(jprofset) ) THEN 619 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 620 & ln_t3d, ln_s3d, ln_nea, & 621 & kdailyavtypes = endailyavtypes ) 622 ELSE 623 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 624 & ln_t3d, ln_s3d, ln_nea ) 625 ENDIF 626 627 END DO 628 629 ENDIF 630 631 ENDIF 632 633 ! - Sea level anomalies 634 IF ( ln_sla ) THEN 635 ! Set the number of variables for sla to 1 636 nslavars = 1 637 638 ! Set the number of extra variables for sla to 2 639 nslaextr = 2 640 641 ! Set the number of sla data sets to 2 642 nslasets = 0 643 IF ( ln_sladt ) THEN 644 nslasets = nslasets + 2 645 ENDIF 646 IF ( ln_slafb ) THEN 647 nslasets = nslasets + jnumslafb 648 ENDIF 649 650 ALLOCATE(sladata(nslasets)) 651 ALLOCATE(sladatqc(nslasets)) 652 sladata(:)%nsurf=0 653 sladatqc(:)%nsurf=0 654 655 nslasets = 0 656 657 ! AVISO SLA data 658 659 IF ( ln_sladt ) THEN 660 661 ! Active SLA observations 662 663 nslasets = nslasets + 1 664 665 CALL obs_rea_sla( 1, sladata(nslasets), jnumslaact, & 666 & slafilesact(1:jnumslaact), & 667 & nslavars, nslaextr, nitend-nit000+2, & 668 & dobsini, dobsend, ln_ignmis, .FALSE. ) 669 CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), & 670 & ln_sla, ln_nea ) 671 672 ! Passive SLA observations 673 674 nslasets = nslasets + 1 675 676 CALL obs_rea_sla( 1, sladata(nslasets), jnumslapas, & 677 & slafilespas(1:jnumslapas), & 678 & nslavars, nslaextr, nitend-nit000+2, & 679 & dobsini, dobsend, ln_ignmis, .FALSE. ) 680 681 CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), & 682 & ln_sla, ln_nea ) 683 684 ENDIF 685 686 ! Feedback SLA data 687 688 IF ( ln_slafb ) THEN 689 690 DO jset = 1, jnumslafb 691 692 nslasets = nslasets + 1 693 694 CALL obs_rea_sla( 0, sladata(nslasets), 1, & 695 & slafbfiles(jset:jset), & 696 & nslavars, nslaextr, nitend-nit000+2, & 697 & dobsini, dobsend, ln_ignmis, .FALSE. ) 698 CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), & 699 & ln_sla, ln_nea ) 700 701 END DO 702 703 ENDIF 704 705 CALL obs_rea_mdt( nslasets, sladatqc, n2dint ) 706 707 ! read in altimeter bias 708 709 IF ( ln_altbias ) THEN 710 CALL obs_rea_altbias ( nslasets, sladatqc, n2dint, bias_file ) 711 ENDIF 712 713 ENDIF 714 715 ! - Sea surface height 716 IF ( ln_ssh ) THEN 717 IF(lwp) WRITE(numout,*) ' SSH currently not available' 718 ENDIF 719 720 ! - Sea surface temperature 721 IF ( ln_sst ) THEN 722 723 ! Set the number of variables for sst to 1 724 nsstvars = 1 725 726 ! Set the number of extra variables for sst to 0 727 nsstextr = 0 728 729 nsstsets = 0 730 731 IF (ln_reysst) nsstsets = nsstsets + 1 732 IF (ln_ghrsst) nsstsets = nsstsets + 1 733 IF ( ln_sstfb ) THEN 734 nsstsets = nsstsets + jnumsstfb 735 ENDIF 736 737 ALLOCATE(sstdata(nsstsets)) 738 ALLOCATE(sstdatqc(nsstsets)) 739 ALLOCATE(ld_sstnight(nsstsets)) 740 sstdata(:)%nsurf=0 741 sstdatqc(:)%nsurf=0 742 ld_sstnight(:)=.false. 743 744 nsstsets = 0 745 746 IF (ln_reysst) THEN 747 748 nsstsets = nsstsets + 1 749 750 ld_sstnight(nsstsets) = ln_sstnight 751 752 CALL obs_rea_sst_rey( reysstname, reysstfmt, sstdata(nsstsets), & 753 & nsstvars, nsstextr, & 754 & nitend-nit000+2, dobsini, dobsend ) 755 CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, & 756 & ln_nea ) 757 758 ENDIF 759 760 IF (ln_ghrsst) THEN 761 762 nsstsets = nsstsets + 1 763 764 ld_sstnight(nsstsets) = ln_sstnight 765 766 CALL obs_rea_sst( 1, sstdata(nsstsets), jnumsst, & 767 & sstfiles(1:jnumsst), & 768 & nsstvars, nsstextr, nitend-nit000+2, & 769 & dobsini, dobsend, ln_ignmis, .FALSE. ) 770 CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, & 771 & ln_nea ) 772 773 ENDIF 774 775 ! Feedback SST data 776 777 IF ( ln_sstfb ) THEN 778 779 DO jset = 1, jnumsstfb 780 781 nsstsets = nsstsets + 1 782 783 ld_sstnight(nsstsets) = ln_sstnight 784 785 CALL obs_rea_sst( 0, sstdata(nsstsets), 1, & 786 & sstfbfiles(jset:jset), & 787 & nsstvars, nsstextr, nitend-nit000+2, & 788 & dobsini, dobsend, ln_ignmis, .FALSE. ) 789 CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), & 790 & ln_sst, ln_nea ) 791 792 END DO 793 794 ENDIF 436 CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), & 437 & llvar1, llvar2, & 438 & jpi, jpj, jpk, & 439 & zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2, & 440 & ln_nea, kdailyavtypes = nn_profdavtypes ) 441 442 END DO 443 444 DEALLOCATE( ifilesprof, clproffiles ) 445 446 ENDIF 447 448 IF ( nsurftypes > 0 ) THEN 449 450 ALLOCATE(surfdata(nsurftypes)) 451 ALLOCATE(surfdataqc(nsurftypes)) 452 ALLOCATE(nvarssurf(nsurftypes)) 453 ALLOCATE(nextrsurf(nsurftypes)) 454 455 DO jtype = 1, nsurftypes 456 457 nvarssurf(jtype) = 1 458 nextrsurf(jtype) = 0 459 llnightav = .FALSE. 460 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2 461 IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav = ln_sstnight 462 463 !Read in surface obs types 464 CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), & 465 & clsurffiles(jtype,1:ifilessurf(jtype)), & 466 & nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 467 & rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav ) 795 468 796 469 !Read in bias field and correct SST. … … 804 477 ENDIF 805 478 806 ENDIF 807 808 ! - Sea surface salinity 809 IF ( ln_sss ) THEN 810 IF(lwp) WRITE(numout,*) ' SSS currently not available' 811 ENDIF 812 813 ! - Sea Ice Concentration 814 815 IF ( ln_seaice ) THEN 816 817 ! Set the number of variables for seaice to 1 818 nseaicevars = 1 819 820 ! Set the number of extra variables for seaice to 0 821 nseaiceextr = 0 822 823 ! Set the number of data sets to 1 824 nseaicesets = 1 825 826 ALLOCATE(seaicedata(nseaicesets)) 827 ALLOCATE(seaicedatqc(nseaicesets)) 828 seaicedata(:)%nsurf=0 829 seaicedatqc(:)%nsurf=0 830 831 CALL obs_rea_seaice( 1, seaicedata(nseaicesets), jnumseaice, & 832 & seaicefiles(1:jnumseaice), & 833 & nseaicevars, nseaiceextr, nitend-nit000+2, & 834 & dobsini, dobsend, ln_ignmis, .FALSE. ) 835 836 CALL obs_pre_seaice( seaicedata(nseaicesets), seaicedatqc(nseaicesets), & 837 & ln_seaice, ln_nea ) 838 839 ENDIF 840 841 IF (ln_vel3d) THEN 842 843 ! Set the number of variables for profiles to 2 (U and V) 844 nvelovars = 2 845 846 ! Set the number of extra variables for profiles to 2 to store 847 ! rotation parameters 848 nveloextr = 2 849 850 jveloset = 0 851 852 IF ( ln_velavcur ) jveloset = jveloset + 1 853 IF ( ln_velhrcur ) jveloset = jveloset + 1 854 IF ( ln_velavadcp ) jveloset = jveloset + 1 855 IF ( ln_velhradcp ) jveloset = jveloset + 1 856 IF (ln_velfb) jveloset = jveloset + jnumvelfb 857 858 nvelosets = jveloset 859 IF ( nvelosets > 0 ) THEN 860 ALLOCATE( velodata(nvelosets) ) 861 ALLOCATE( veldatqc(nvelosets) ) 862 ALLOCATE( ld_velav(nvelosets) ) 863 ENDIF 864 865 jveloset = 0 866 867 ! Daily averaged data 868 869 IF ( ln_velavcur ) THEN 870 871 jveloset = jveloset + 1 872 873 ld_velav(jveloset) = .TRUE. 874 875 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavcur, & 876 & velavcurfiles(1:jnumvelavcur), & 877 & nvelovars, nveloextr, & 878 & nitend-nit000+2, & 879 & dobsini, dobsend, ln_ignmis, & 880 & ld_velav(jveloset), & 881 & .FALSE. ) 882 883 DO jvar = 1, 2 884 CALL obs_prof_staend( velodata(jveloset), jvar ) 885 END DO 886 887 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 888 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 889 890 ENDIF 891 892 ! High frequency data 893 894 IF ( ln_velhrcur ) THEN 895 896 jveloset = jveloset + 1 897 898 ld_velav(jveloset) = .FALSE. 899 900 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhrcur, & 901 & velhrcurfiles(1:jnumvelhrcur), & 902 & nvelovars, nveloextr, & 903 & nitend-nit000+2, & 904 & dobsini, dobsend, ln_ignmis, & 905 & ld_velav(jveloset), & 906 & .FALSE. ) 907 908 DO jvar = 1, 2 909 CALL obs_prof_staend( velodata(jveloset), jvar ) 910 END DO 911 912 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 913 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 914 915 ENDIF 916 917 ! Daily averaged data 918 919 IF ( ln_velavadcp ) THEN 920 921 jveloset = jveloset + 1 922 923 ld_velav(jveloset) = .TRUE. 924 925 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavadcp, & 926 & velavadcpfiles(1:jnumvelavadcp), & 927 & nvelovars, nveloextr, & 928 & nitend-nit000+2, & 929 & dobsini, dobsend, ln_ignmis, & 930 & ld_velav(jveloset), & 931 & .FALSE. ) 932 933 DO jvar = 1, 2 934 CALL obs_prof_staend( velodata(jveloset), jvar ) 935 END DO 936 937 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 938 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 939 940 ENDIF 941 942 ! High frequency data 943 944 IF ( ln_velhradcp ) THEN 945 946 jveloset = jveloset + 1 947 948 ld_velav(jveloset) = .FALSE. 949 950 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhradcp, & 951 & velhradcpfiles(1:jnumvelhradcp), & 952 & nvelovars, nveloextr, & 953 & nitend-nit000+2, & 954 & dobsini, dobsend, ln_ignmis, & 955 & ld_velav(jveloset), & 956 & .FALSE. ) 957 958 DO jvar = 1, 2 959 CALL obs_prof_staend( velodata(jveloset), jvar ) 960 END DO 961 962 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 963 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 964 965 ENDIF 966 967 IF ( ln_velfb ) THEN 968 969 DO jset = 1, jnumvelfb 970 971 jveloset = jveloset + 1 972 973 ld_velav(jveloset) = ln_velfb_av(jset) 974 975 CALL obs_rea_vel_dri( 0, velodata(jveloset), 1, & 976 & velfbfiles(jset:jset), & 977 & nvelovars, nveloextr, & 978 & nitend-nit000+2, & 979 & dobsini, dobsend, ln_ignmis, & 980 & ld_velav(jveloset), & 981 & .FALSE. ) 982 983 DO jvar = 1, 2 984 CALL obs_prof_staend( velodata(jveloset), jvar ) 985 END DO 986 987 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 988 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 989 990 991 END DO 992 993 ENDIF 994 995 ENDIF 996 479 CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea ) 480 481 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 482 CALL obs_rea_mdt( surfdataqc(jtype), nn_2dint ) 483 IF ( ln_altbias ) CALL obs_rea_altbias ( surfdataqc(jtype), nn_2dint, cn_altbiasfile ) 484 ENDIF 485 486 END DO 487 488 DEALLOCATE( ifilessurf, clsurffiles ) 489 490 ENDIF 491 492 CALL wrk_dealloc( jpi, jpj, zglam1 ) 493 CALL wrk_dealloc( jpi, jpj, zglam2 ) 494 CALL wrk_dealloc( jpi, jpj, zgphi1 ) 495 CALL wrk_dealloc( jpi, jpj, zgphi2 ) 496 CALL wrk_dealloc( jpi, jpj, jpk, zmask1 ) 497 CALL wrk_dealloc( jpi, jpj, jpk, zmask2 ) 498 997 499 END SUBROUTINE dia_obs_init 998 500 … … 1004 506 !! 1005 507 !! ** Method : Call the observation operators on each time step to 1006 !! compute the model equivalent of the following date: 1007 !! - T profiles 1008 !! - S profiles 1009 !! - Sea surface height (referenced to a mean) 1010 !! - Sea surface temperature 1011 !! - Sea surface salinity 1012 !! - Velocity component (U,V) profiles 1013 !! 1014 !! ** Action : 508 !! compute the model equivalent of the following data: 509 !! - Profile data, currently T/S or U/V 510 !! - Surface data, currently SST, SLA or sea-ice concentration. 511 !! 512 !! ** Action : 1015 513 !! 1016 514 !! History : … … 1023 521 !! ! 14-08 (J. While) observation operator for profiles in 1024 522 !! generalised vertical coordinates 523 !! ! 15-08 (M. Martin) Combined surface/profile routines. 1025 524 !!---------------------------------------------------------------------- 1026 525 !! * Modules used 1027 526 USE dom_oce, ONLY : & ! Ocean space and time domain variables 1028 & rdt, &1029 & gdept_1d, &1030 527 #if defined key_vvl 1031 528 & gdept_n, & … … 1033 530 & gdept_1d, & 1034 531 #endif 1035 & tmask, umask, vmask1036 532 USE phycst, ONLY : & ! Physical constants 1037 533 & rday … … 1039 535 & tsn, & 1040 536 & un, vn, & 537 USE phycst, ONLY : & ! Physical constants 538 & rday 1041 539 & sshn 1042 540 #if defined key_lim3 1043 USE ice, ONLY : & ! LIMIce model variables541 USE ice, ONLY : & ! LIM3 Ice model variables 1044 542 & frld 1045 543 #endif 1046 544 #if defined key_lim2 1047 USE ice_2, ONLY : & ! LIMIce model variables545 USE ice_2, ONLY : & ! LIM2 Ice model variables 1048 546 & frld 1049 547 #endif … … 1051 549 1052 550 !! * Arguments 1053 INTEGER, INTENT(IN) :: kstp 551 INTEGER, INTENT(IN) :: kstp ! Current timestep 1054 552 !! * Local declarations 1055 INTEGER :: idaystp ! Number of timesteps per day 1056 INTEGER :: jprofset ! Profile data set loop variable 1057 INTEGER :: jslaset ! SLA data set loop variable 1058 INTEGER :: jsstset ! SST data set loop variable 1059 INTEGER :: jseaiceset ! sea ice data set loop variable 1060 INTEGER :: jveloset ! velocity profile data loop variable 1061 INTEGER :: jvar ! Variable number 553 INTEGER :: idaystp ! Number of timesteps per day 554 INTEGER :: jtype ! Data loop variable 555 INTEGER :: jvar ! Variable number 556 INTEGER :: ji, jj ! Loop counters 557 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 558 & zprofvar1, & ! Model values for 1st variable in a prof ob 559 & zprofvar2 ! Model values for 2nd variable in a prof ob 560 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 561 & zprofmask1, & ! Mask associated with zprofvar1 562 & zprofmask2 ! Mask associated with zprofvar2 563 REAL(wp), POINTER, DIMENSION(:,:) :: & 564 & zsurfvar ! Model values equivalent to surface ob. 565 REAL(wp), POINTER, DIMENSION(:,:) :: & 566 & zglam1, & ! Model longitudes for prof variable 1 567 & zglam2, & ! Model longitudes for prof variable 2 568 & zgphi1, & ! Model latitudes for prof variable 1 569 & zgphi2 ! Model latitudes for prof variable 2 1062 570 #if ! defined key_lim2 && ! defined key_lim3 1063 REAL(wp), POINTER, DIMENSION(:,:) :: frld 571 REAL(wp), POINTER, DIMENSION(:,:) :: frld 1064 572 #endif 1065 CHARACTER(LEN=20) :: datestr=" ",timestr=" " 1066 573 LOGICAL :: llnightav ! Logical for calculating night-time average 574 575 !Allocate local work arrays 576 CALL wrk_alloc( jpi, jpj, jpk, zprofvar1 ) 577 CALL wrk_alloc( jpi, jpj, jpk, zprofvar2 ) 578 CALL wrk_alloc( jpi, jpj, jpk, zprofmask1 ) 579 CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 ) 580 CALL wrk_alloc( jpi, jpj, zsurfvar ) 581 CALL wrk_alloc( jpi, jpj, zglam1 ) 582 CALL wrk_alloc( jpi, jpj, zglam2 ) 583 CALL wrk_alloc( jpi, jpj, zgphi1 ) 584 CALL wrk_alloc( jpi, jpj, zgphi2 ) 1067 585 #if ! defined key_lim2 && ! defined key_lim3 1068 586 CALL wrk_alloc(jpi,jpj,frld) … … 1084 602 #endif 1085 603 !----------------------------------------------------------------------- 1086 ! Depending on switches call various observation operators 1087 !----------------------------------------------------------------------- 1088 604 ! Call the profile and surface observation operators 605 !----------------------------------------------------------------------- 606 607 <<<<<<< .working 1089 608 ! - Temperature/salinity profiles 1090 609 IF ( ln_t3d .OR. ln_s3d ) THEN … … 1121 640 END DO 1122 641 ENDIF 1123 1124 ! - Sea surface anomaly 1125 IF ( ln_sla ) THEN 1126 DO jslaset = 1, nslasets 1127 CALL obs_sla_opt( sladatqc(jslaset), & 1128 & kstp, jpi, jpj, nit000, sshn, & 1129 & tmask(:,:,1), n2dint ) 1130 END DO 1131 ENDIF 1132 1133 ! - Sea surface temperature 1134 IF ( ln_sst ) THEN 1135 DO jsstset = 1, nsstsets 1136 CALL obs_sst_opt( sstdatqc(jsstset), & 1137 & kstp, jpi, jpj, nit000, idaystp, & 1138 & tsn(:,:,1,jp_tem), tmask(:,:,1), & 1139 & n2dint, ld_sstnight(jsstset) ) 642 ======= 643 IF ( nproftypes > 0 ) THEN 644 >>>>>>> .merge-right.r5997 645 646 DO jtype = 1, nproftypes 647 648 SELECT CASE ( TRIM(cobstypesprof(jtype)) ) 649 CASE('prof') 650 zprofvar1(:,:,:) = tsn(:,:,:,jp_tem) 651 zprofvar2(:,:,:) = tsn(:,:,:,jp_sal) 652 zprofmask1(:,:,:) = tmask(:,:,:) 653 zprofmask2(:,:,:) = tmask(:,:,:) 654 zglam1(:,:) = glamt(:,:) 655 zglam2(:,:) = glamt(:,:) 656 zgphi1(:,:) = gphit(:,:) 657 zgphi2(:,:) = gphit(:,:) 658 CASE('vel') 659 zprofvar1(:,:,:) = un(:,:,:) 660 zprofvar2(:,:,:) = vn(:,:,:) 661 zprofmask1(:,:,:) = umask(:,:,:) 662 zprofmask2(:,:,:) = vmask(:,:,:) 663 zglam1(:,:) = glamu(:,:) 664 zglam2(:,:) = glamv(:,:) 665 zgphi1(:,:) = gphiu(:,:) 666 zgphi2(:,:) = gphiv(:,:) 667 END SELECT 668 669 IF( ln_zco .OR. ln_zps ) THEN 670 CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk, & 671 & nit000, idaystp, & 672 & zprofvar1, zprofvar2, & 673 & gdept_1d, zprofmask1, zprofmask2, & 674 & zglam1, zglam2, zgphi1, zgphi2, & 675 & nn_1dint, nn_2dint, & 676 & kdailyavtypes = nn_profdavtypes ) 677 ELSE IF(TRIM(cobstypesprof(jtype)) == 'prof') 678 !TG - THIS NEEDS MODIFICATION TO MATCH SIMPLIFICATION 679 CALL obs_pro_sco_opt( prodatqc(jtype), & 680 & kstp, jpi, jpj, jpk, nit000, idaystp, & 681 & zprofvar1, zprofvar2, & 682 & fsdept(:,:,:), fsdepw(:,:,:), & 683 & tmask, nn_1dint, nn_2dint, & 684 & kdailyavtypes = nn_profdavtypes ) 685 ELSE 686 ctl_stop('DIA_OBS: Generalised vertical interpolation not'// & 687 'yet working for velocity date (turn off velocity observations') 688 ENDIF 689 1140 690 END DO 1141 ENDIF 1142 1143 ! - Sea surface salinity 1144 IF ( ln_sss ) THEN 1145 IF(lwp) WRITE(numout,*) ' SSS currently not available' 1146 ENDIF 1147 691 692 ENDIF 693 694 IF ( nsurftypes > 0 ) THEN 695 696 DO jtype = 1, nsurftypes 697 698 SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 699 CASE('sst') 700 zsurfvar(:,:) = tsn(:,:,1,jp_tem) 701 llnightav = ln_sstnight 702 CASE('sla') 703 zsurfvar(:,:) = sshn(:,:) 704 llnightav = .FALSE. 1148 705 #if defined key_lim2 || defined key_lim3 1149 IF ( ln_seaice ) THEN 1150 DO jseaiceset = 1, nseaicesets 1151 CALL obs_seaice_opt( seaicedatqc(jseaiceset), & 1152 & kstp, jpi, jpj, nit000, 1.-frld, & 1153 & tmask(:,:,1), n2dint ) 706 CASE('sic') 707 IF ( kstp == 0 ) THEN 708 IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN 709 CALL ctl_warn( 'Sea-ice not initialised on zeroth '// & 710 & 'time-step but some obs are valid then.' ) 711 WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & 712 & ' sea-ice obs will be missed' 713 ENDIF 714 surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + & 715 & surfdataqc(jtype)%nsstp(1) 716 CYCLE 717 ELSE 718 zsurfvar(:,:) = 1._wp - frld(:,:) 719 ENDIF 720 721 llnightav = .FALSE. 722 #endif 723 END SELECT 724 725 CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj, & 726 & nit000, idaystp, zsurfvar, tmask(:,:,1), & 727 & nn_2dint, llnightav ) 728 1154 729 END DO 1155 ENDIF 730 731 ENDIF 732 733 CALL wrk_dealloc( jpi, jpj, jpk, zprofvar1 ) 734 CALL wrk_dealloc( jpi, jpj, jpk, zprofvar2 ) 735 CALL wrk_dealloc( jpi, jpj, jpk, zprofmask1 ) 736 CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 ) 737 CALL wrk_dealloc( jpi, jpj, zsurfvar ) 738 CALL wrk_dealloc( jpi, jpj, zglam1 ) 739 CALL wrk_dealloc( jpi, jpj, zglam2 ) 740 CALL wrk_dealloc( jpi, jpj, zgphi1 ) 741 CALL wrk_dealloc( jpi, jpj, zgphi2 ) 742 #if ! defined key_lim2 && ! defined key_lim3 743 CALL wrk_dealloc(jpi,jpj,frld) 1156 744 #endif 1157 745 1158 ! - Velocity profiles1159 IF ( ln_vel3d ) THEN1160 DO jveloset = 1, nvelosets1161 ! zonal component of velocity1162 CALL obs_vel_opt( veldatqc(jveloset), kstp, jpi, jpj, jpk, &1163 & nit000, idaystp, un, vn, gdept_1d, umask, vmask, &1164 n1dint, n2dint, ld_velav(jveloset) )1165 END DO1166 ENDIF1167 1168 #if ! defined key_lim2 && ! defined key_lim31169 CALL wrk_dealloc(jpi,jpj,frld)1170 #endif1171 1172 746 END SUBROUTINE dia_obs 1173 1174 SUBROUTINE dia_obs_wri 747 748 SUBROUTINE dia_obs_wri 1175 749 !!---------------------------------------------------------------------- 1176 750 !! *** ROUTINE dia_obs_wri *** … … 1180 754 !! ** Method : Call observation diagnostic output routines 1181 755 !! 1182 !! ** Action : 756 !! ** Action : 1183 757 !! 1184 758 !! History : … … 1188 762 !! ! 07-03 (K. Mogensen) General handling of profiles 1189 763 !! ! 08-09 (M. Valdivieso) Velocity component (U,V) profiles 1190 !!---------------------------------------------------------------------- 764 !! ! 15-08 (M. Martin) Combined writing for prof and surf types 765 !!---------------------------------------------------------------------- 766 !! * Modules used 767 USE obs_rot_vel ! Rotation of velocities 768 1191 769 IMPLICIT NONE 1192 770 1193 771 !! * Local declarations 1194 1195 INTEGER :: jprofset ! Profile data set loop variable 1196 INTEGER :: jveloset ! Velocity data set loop variable 1197 INTEGER :: jslaset ! SLA data set loop variable 1198 INTEGER :: jsstset ! SST data set loop variable 1199 INTEGER :: jseaiceset ! Sea Ice data set loop variable 1200 INTEGER :: jset 1201 INTEGER :: jfbini 1202 CHARACTER(LEN=20) :: datestr=" ",timestr=" " 1203 CHARACTER(LEN=10) :: cdtmp 772 INTEGER :: jtype ! Data set loop variable 773 INTEGER :: jo, jvar, jk 774 REAL(wp), DIMENSION(:), ALLOCATABLE :: & 775 & zu, & 776 & zv 777 1204 778 !----------------------------------------------------------------------- 1205 779 ! Depending on switches call various observation output routines 1206 780 !----------------------------------------------------------------------- 1207 781 1208 ! - Temperature/salinity profiles 1209 1210 IF( ln_t3d .OR. ln_s3d ) THEN 1211 1212 ! Copy data from prodatqc to profdata structures 1213 DO jprofset = 1, nprofsets 1214 1215 CALL obs_prof_decompress( prodatqc(jprofset), & 1216 & profdata(jprofset), .TRUE., numout ) 782 IF ( nproftypes > 0 ) THEN 783 784 DO jtype = 1, nproftypes 785 786 IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN 787 788 ! For velocity data, rotate the model velocities to N/S, E/W 789 ! using the compressed data structure. 790 ALLOCATE( & 791 & zu(profdataqc(jtype)%nvprot(1)), & 792 & zv(profdataqc(jtype)%nvprot(2)) & 793 & ) 794 795 CALL obs_rotvel( profdataqc(jtype), nn_2dint, zu, zv ) 796 797 DO jo = 1, profdataqc(jtype)%nprof 798 DO jvar = 1, 2 799 DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar) 800 801 IF ( jvar == 1 ) THEN 802 profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk) 803 ELSE 804 profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk) 805 ENDIF 806 807 END DO 808 END DO 809 END DO 810 811 DEALLOCATE( zu ) 812 DEALLOCATE( zv ) 813 814 END IF 815 816 CALL obs_prof_decompress( profdataqc(jtype), & 817 & profdata(jtype), .TRUE., numout ) 818 819 CALL obs_wri_prof( profdata(jtype) ) 1217 820 1218 821 END DO 1219 822 1220 ! Write the profiles. 1221 1222 jprofset = 0 1223 1224 ! ENACT insitu data 1225 1226 IF ( ln_ena ) THEN 1227 1228 jprofset = jprofset + 1 1229 1230 CALL obs_wri_p3d( 'enact', profdata(jprofset) ) 1231 1232 ENDIF 1233 1234 ! Coriolis insitu data 1235 1236 IF ( ln_cor ) THEN 1237 1238 jprofset = jprofset + 1 1239 1240 CALL obs_wri_p3d( 'corio', profdata(jprofset) ) 1241 1242 ENDIF 1243 1244 ! Feedback insitu data 1245 1246 IF ( ln_profb ) THEN 1247 1248 jfbini = jprofset + 1 1249 1250 DO jprofset = jfbini, nprofsets 1251 1252 jset = jprofset - jfbini + 1 1253 WRITE(cdtmp,'(A,I2.2)')'profb_',jset 1254 CALL obs_wri_p3d( cdtmp, profdata(jprofset) ) 1255 1256 END DO 1257 1258 ENDIF 1259 1260 ENDIF 1261 1262 ! - Sea surface anomaly 1263 IF ( ln_sla ) THEN 1264 1265 ! Copy data from sladatqc to sladata structures 1266 DO jslaset = 1, nslasets 1267 1268 CALL obs_surf_decompress( sladatqc(jslaset), & 1269 & sladata(jslaset), .TRUE., numout ) 823 ENDIF 824 825 IF ( nsurftypes > 0 ) THEN 826 827 DO jtype = 1, nsurftypes 828 829 CALL obs_surf_decompress( surfdataqc(jtype), & 830 & surfdata(jtype), .TRUE., numout ) 831 832 CALL obs_wri_surf( surfdata(jtype) ) 1270 833 1271 834 END DO 1272 835 1273 jslaset = 01274 1275 ! Write the AVISO SLA data1276 1277 IF ( ln_sladt ) THEN1278 1279 jslaset = 11280 CALL obs_wri_sla( 'aviso_act', sladata(jslaset) )1281 jslaset = 21282 CALL obs_wri_sla( 'aviso_pas', sladata(jslaset) )1283 1284 ENDIF1285 1286 IF ( ln_slafb ) THEN1287 1288 jfbini = jslaset + 11289 1290 DO jslaset = jfbini, nslasets1291 1292 jset = jslaset - jfbini + 11293 WRITE(cdtmp,'(A,I2.2)')'slafb_',jset1294 CALL obs_wri_sla( cdtmp, sladata(jslaset) )1295 1296 END DO1297 1298 ENDIF1299 1300 ENDIF1301 1302 ! - Sea surface temperature1303 IF ( ln_sst ) THEN1304 1305 ! Copy data from sstdatqc to sstdata structures1306 DO jsstset = 1, nsstsets1307 1308 CALL obs_surf_decompress( sstdatqc(jsstset), &1309 & sstdata(jsstset), .TRUE., numout )1310 1311 END DO1312 1313 jsstset = 01314 1315 ! Write the AVISO SST data1316 1317 IF ( ln_reysst ) THEN1318 1319 jsstset = jsstset + 11320 CALL obs_wri_sst( 'reynolds', sstdata(jsstset) )1321 1322 ENDIF1323 1324 IF ( ln_ghrsst ) THEN1325 1326 jsstset = jsstset + 11327 CALL obs_wri_sst( 'ghr', sstdata(jsstset) )1328 1329 ENDIF1330 1331 IF ( ln_sstfb ) THEN1332 1333 jfbini = jsstset + 11334 1335 DO jsstset = jfbini, nsstsets1336 1337 jset = jsstset - jfbini + 11338 WRITE(cdtmp,'(A,I2.2)')'sstfb_',jset1339 CALL obs_wri_sst( cdtmp, sstdata(jsstset) )1340 1341 END DO1342 1343 ENDIF1344 1345 ENDIF1346 1347 ! - Sea surface salinity1348 IF ( ln_sss ) THEN1349 IF(lwp) WRITE(numout,*) ' SSS currently not available'1350 ENDIF1351 1352 ! - Sea Ice Concentration1353 IF ( ln_seaice ) THEN1354 1355 ! Copy data from seaicedatqc to seaicedata structures1356 DO jseaiceset = 1, nseaicesets1357 1358 CALL obs_surf_decompress( seaicedatqc(jseaiceset), &1359 & seaicedata(jseaiceset), .TRUE., numout )1360 1361 END DO1362 1363 ! Write the Sea Ice data1364 DO jseaiceset = 1, nseaicesets1365 1366 WRITE(cdtmp,'(A,I2.2)')'seaicefb_',jseaiceset1367 CALL obs_wri_seaice( cdtmp, seaicedata(jseaiceset) )1368 1369 END DO1370 1371 ENDIF1372 1373 ! Velocity data1374 IF( ln_vel3d ) THEN1375 1376 ! Copy data from veldatqc to velodata structures1377 DO jveloset = 1, nvelosets1378 1379 CALL obs_prof_decompress( veldatqc(jveloset), &1380 & velodata(jveloset), .TRUE., numout )1381 1382 END DO1383 1384 ! Write the profiles.1385 1386 jveloset = 01387 1388 ! Daily averaged data1389 1390 IF ( ln_velavcur ) THEN1391 1392 jveloset = jveloset + 11393 1394 CALL obs_wri_vel( 'velavcurr', velodata(jveloset), n2dint )1395 1396 ENDIF1397 1398 ! High frequency data1399 1400 IF ( ln_velhrcur ) THEN1401 1402 jveloset = jveloset + 11403 1404 CALL obs_wri_vel( 'velhrcurr', velodata(jveloset), n2dint )1405 1406 ENDIF1407 1408 ! Daily averaged data1409 1410 IF ( ln_velavadcp ) THEN1411 1412 jveloset = jveloset + 11413 1414 CALL obs_wri_vel( 'velavadcp', velodata(jveloset), n2dint )1415 1416 ENDIF1417 1418 ! High frequency data1419 1420 IF ( ln_velhradcp ) THEN1421 1422 jveloset = jveloset + 11423 1424 CALL obs_wri_vel( 'velhradcp', velodata(jveloset), n2dint )1425 1426 ENDIF1427 1428 ! Feedback velocity data1429 1430 IF ( ln_velfb ) THEN1431 1432 jfbini = jveloset + 11433 1434 DO jveloset = jfbini, nvelosets1435 1436 jset = jveloset - jfbini + 11437 WRITE(cdtmp,'(A,I2.2)')'velfb_',jset1438 CALL obs_wri_vel( cdtmp, velodata(jveloset), n2dint )1439 1440 END DO1441 1442 ENDIF1443 1444 836 ENDIF 1445 837 … … 1459 851 !! 1460 852 !!---------------------------------------------------------------------- 1461 ! !obs_grid deallocation853 ! obs_grid deallocation 1462 854 CALL obs_grid_deallocate 1463 855 1464 !! diaobs deallocation 1465 IF ( nprofsets > 0 ) THEN 1466 DEALLOCATE(ld_enact, & 1467 & profdata, & 1468 & prodatqc) 1469 END IF 1470 IF ( ln_sla ) THEN 1471 DEALLOCATE(sladata, & 1472 & sladatqc) 1473 END IF 1474 IF ( ln_seaice ) THEN 1475 DEALLOCATE(sladata, & 1476 & sladatqc) 1477 END IF 1478 IF ( ln_sst ) THEN 1479 DEALLOCATE(sstdata, & 1480 & sstdatqc) 1481 END IF 1482 IF ( ln_vel3d ) THEN 1483 DEALLOCATE(ld_velav, & 1484 & velodata, & 1485 & veldatqc) 1486 END IF 856 ! diaobs deallocation 857 IF ( nproftypes > 0 ) & 858 & DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof ) 859 860 IF ( nsurftypes > 0 ) & 861 & DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf ) 862 1487 863 END SUBROUTINE dia_obs_dealloc 1488 864 … … 1496 872 !! 1497 873 !! ** Action : Get date in double precision YYYYMMDD.HHMMSS format 874 !! 875 !! ** Action : Get initial date in double precision YYYYMMDD.HHMMSS format 1498 876 !! 1499 877 !! History : … … 1507 885 USE phycst, ONLY : & ! Physical constants 1508 886 & rday 1509 ! USE daymod, ONLY : & ! Time variables1510 ! & nmonth_len1511 887 USE dom_oce, ONLY : & ! Ocean space and time domain variables 1512 888 & rdt … … 1524 900 INTEGER :: ihou 1525 901 INTEGER :: imin 1526 INTEGER :: imday ! Number of days in month. 1527 REAL(KIND=wp) :: zdayfrc ! Fraction of day 902 INTEGER :: imday ! Number of days in month. 903 INTEGER, DIMENSION(12) :: & 904 & imonth_len ! Length in days of the months of the current year 905 REAL(wp) :: zdayfrc ! Fraction of day 1528 906 1529 907 INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year … … 1554 932 iday = iday + kstp * rdt / rday 1555 933 1556 ! !-----------------------------------------------------------------------1557 ! !Convert number of days (iday) into a real date1558 ! !----------------------------------------------------------------------934 !----------------------------------------------------------------------- 935 ! Convert number of days (iday) into a real date 936 !---------------------------------------------------------------------- 1559 937 1560 938 CALL calc_month_len( iyea, imonth_len ) 1561 939 1562 940 DO WHILE ( iday > imonth_len(imon) ) 1563 941 iday = iday - imonth_len(imon) … … 1570 948 END DO 1571 949 1572 ! !----------------------------------------------------------------------1573 ! !Convert it into YYYYMMDD.HHMMSS format.1574 ! !----------------------------------------------------------------------950 !---------------------------------------------------------------------- 951 ! Convert it into YYYYMMDD.HHMMSS format. 952 !---------------------------------------------------------------------- 1575 953 ddobs = iyea * 10000_dp + imon * 100_dp + & 1576 954 & iday + ihou * 0.01_dp + imin * 0.0001_dp … … 1626 1004 1627 1005 !! * Arguments 1628 REAL( KIND=dp), INTENT(OUT) :: ddobsfin! Final date in YYYYMMDD.HHMMSS1006 REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS 1629 1007 1630 1008 CALL calc_date( nitend, ddobsfin )
Note: See TracChangeset
for help on using the changeset viewer.
