- Timestamp:
- 2015-12-16T16:44:35+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
r4990 r6069 7 7 8 8 !!---------------------------------------------------------------------- 9 !! 'key_diaobs' : Switch on the observation diagnostic computation10 !!----------------------------------------------------------------------11 9 !! dia_obs_init : Reading and prepare observations 12 10 !! dia_obs : Compute model equivalent to observations 13 11 !! dia_obs_wri : Write observational diagnostics 12 !! calc_date : Compute the date of timestep in YYYYMMDD.HHMMSS format 14 13 !! ini_date : Compute the initial date YYYYMMDD.HHMMSS 15 14 !! fin_date : Compute the final date YYYYMMDD.HHMMSS 16 15 !!---------------------------------------------------------------------- 17 !! * Modules used 16 !! * Modules used 18 17 USE wrk_nemo ! Memory Allocation 19 18 USE par_kind ! Precision variables … … 21 20 USE par_oce 22 21 USE dom_oce ! Ocean space and time domain variables 23 USE obs_fbm, ONLY: ln_cl4 ! Class 4 diagnostic switch 24 USE obs_read_prof ! Reading and allocation of observations (Coriolis) 25 USE obs_read_sla ! Reading and allocation of SLA observations 26 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 24 USE obs_sstbias ! Bias correction routine for SST 27 25 USE obs_readmdt ! Reading and allocation of MDT for SLA. 28 USE obs_read_seaice ! Reading and allocation of Sea Ice observations29 USE obs_read_vel ! Reading and allocation of velocity component observations30 26 USE obs_prep ! Preparation of obs. (grid search etc). 31 27 USE obs_oper ! Observation operators … … 34 30 USE obs_read_altbias ! Bias treatment for altimeter 35 31 USE obs_profiles_def ! Profile data definitions 36 USE obs_profiles ! Profile data storage37 32 USE obs_surf_def ! Surface data definitions 38 USE obs_sla ! SLA data storage39 USE obs_sst ! SST data storage40 USE obs_seaice ! Sea Ice data storage41 33 USE obs_types ! Definitions for observation types 42 34 USE mpp_map ! MPP mapping … … 50 42 & dia_obs, & ! Compute model equivalent to observations 51 43 & dia_obs_wri, & ! Write model equivalent to observations 52 & dia_obs_dealloc ! Deallocate dia_obs data 53 54 !! * Shared Module variables 55 LOGICAL, PUBLIC, PARAMETER :: & 56 #if defined key_diaobs 57 & lk_diaobs = .TRUE. !: Logical switch for observation diangostics 58 #else 59 & lk_diaobs = .FALSE. !: Logical switch for observation diangostics 60 #endif 44 & dia_obs_dealloc, & ! Deallocate dia_obs data 45 & calc_date ! Compute the date of a timestep 61 46 62 47 !! * Module variables 63 LOGICAL, PUBLIC :: ln_t3d !: Logical switch for temperature profiles 64 LOGICAL, PUBLIC :: ln_s3d !: Logical switch for salinity profiles 65 LOGICAL, PUBLIC :: ln_ena !: Logical switch for the ENACT data set 66 LOGICAL, PUBLIC :: ln_cor !: Logical switch for the Coriolis data set 67 LOGICAL, PUBLIC :: ln_profb !: Logical switch for profile feedback datafiles 68 LOGICAL, PUBLIC :: ln_sla !: Logical switch for sea level anomalies 69 LOGICAL, PUBLIC :: ln_sladt !: Logical switch for SLA from AVISO files 70 LOGICAL, PUBLIC :: ln_slafb !: Logical switch for SLA from feedback files 71 LOGICAL, PUBLIC :: ln_sst !: Logical switch for sea surface temperature 72 LOGICAL, PUBLIC :: ln_reysst !: Logical switch for Reynolds sea surface temperature 73 LOGICAL, PUBLIC :: ln_ghrsst !: Logical switch for GHRSST data 74 LOGICAL, PUBLIC :: ln_sstfb !: Logical switch for SST from feedback files 75 LOGICAL, PUBLIC :: ln_seaice !: Logical switch for sea ice concentration 76 LOGICAL, PUBLIC :: ln_vel3d !: Logical switch for velocity component (u,v) observations 77 LOGICAL, PUBLIC :: ln_velavcur !: Logical switch for raw daily averaged netCDF current meter vel. data 78 LOGICAL, PUBLIC :: ln_velhrcur !: Logical switch for raw high freq netCDF current meter vel. data 79 LOGICAL, PUBLIC :: ln_velavadcp !: Logical switch for raw daily averaged netCDF ADCP vel. data 80 LOGICAL, PUBLIC :: ln_velhradcp !: Logical switch for raw high freq netCDF ADCP vel. data 81 LOGICAL, PUBLIC :: ln_velfb !: Logical switch for velocities from feedback files 82 LOGICAL, PUBLIC :: ln_ssh !: Logical switch for sea surface height 83 LOGICAL, PUBLIC :: ln_sss !: Logical switch for sea surface salinity 84 LOGICAL, PUBLIC :: ln_sstnight !: Logical switch for night mean SST observations 85 LOGICAL, PUBLIC :: ln_nea !: Remove observations near land 86 LOGICAL, PUBLIC :: ln_altbias !: Logical switch for altimeter bias 87 LOGICAL, PUBLIC :: ln_ignmis !: Logical switch for ignoring missing files 88 LOGICAL, PUBLIC :: ln_s_at_t !: Logical switch to compute model S at T observations 89 90 REAL(KIND=dp), PUBLIC :: dobsini !: Observation window start date YYYYMMDD.HHMMSS 91 REAL(KIND=dp), PUBLIC :: dobsend !: Observation window end date YYYYMMDD.HHMMSS 92 93 INTEGER, PUBLIC :: n1dint !: Vertical interpolation method 94 INTEGER, PUBLIC :: n2dint !: Horizontal interpolation method 95 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 96 53 INTEGER, DIMENSION(imaxavtypes) :: & 97 & endailyavtypes !: ENACT data types which are daily average 98 99 INTEGER, PARAMETER :: MaxNumFiles = 1000 100 LOGICAL, DIMENSION(MaxNumFiles) :: & 101 & ln_profb_ena, & !: Is the feedback files from ENACT data ? 102 ! !: If so use endailyavtypes 103 & ln_profb_enatim !: Change tim for 820 enact data set. 104 105 LOGICAL, DIMENSION(MaxNumFiles) :: & 106 & ln_velfb_av !: Is the velocity feedback files daily average? 107 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 108 & ld_enact !: Profile data is ENACT so use endailyavtypes 109 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 110 & ld_velav !: Velocity data is daily averaged 111 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 112 & 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 113 74 114 75 !!---------------------------------------------------------------------- … … 118 79 !!---------------------------------------------------------------------- 119 80 81 !! * Substitutions 82 # include "domzgr_substitute.h90" 120 83 CONTAINS 121 84 … … 135 98 !! ! 06-10 (A. Weaver) Cleaning and add controls 136 99 !! ! 07-03 (K. Mogensen) General handling of profiles 100 !! ! 14-08 (J.While) Incorporated SST bias correction 101 !! ! 15-02 (M. Martin) Simplification of namelist and code 137 102 !!---------------------------------------------------------------------- 138 103 … … 140 105 141 106 !! * Local declarations 142 CHARACTER(len=128) :: enactfiles(MaxNumFiles) 143 CHARACTER(len=128) :: coriofiles(MaxNumFiles) 144 CHARACTER(len=128) :: profbfiles(MaxNumFiles) 145 CHARACTER(len=128) :: sstfiles(MaxNumFiles) 146 CHARACTER(len=128) :: sstfbfiles(MaxNumFiles) 147 CHARACTER(len=128) :: slafilesact(MaxNumFiles) 148 CHARACTER(len=128) :: slafilespas(MaxNumFiles) 149 CHARACTER(len=128) :: slafbfiles(MaxNumFiles) 150 CHARACTER(len=128) :: seaicefiles(MaxNumFiles) 151 CHARACTER(len=128) :: velcurfiles(MaxNumFiles) 152 CHARACTER(len=128) :: veladcpfiles(MaxNumFiles) 153 CHARACTER(len=128) :: velavcurfiles(MaxNumFiles) 154 CHARACTER(len=128) :: velhrcurfiles(MaxNumFiles) 155 CHARACTER(len=128) :: velavadcpfiles(MaxNumFiles) 156 CHARACTER(len=128) :: velhradcpfiles(MaxNumFiles) 157 CHARACTER(len=128) :: velfbfiles(MaxNumFiles) 158 CHARACTER(LEN=128) :: reysstname 159 CHARACTER(LEN=12) :: reysstfmt 160 CHARACTER(LEN=128) :: bias_file 161 CHARACTER(LEN=20) :: datestr=" ", timestr=" " 162 NAMELIST/namobs/ln_ena, ln_cor, ln_profb, ln_t3d, ln_s3d, & 163 & ln_sla, ln_sladt, ln_slafb, & 164 & ln_ssh, ln_sst, ln_sstfb, ln_sss, ln_nea, & 165 & enactfiles, coriofiles, profbfiles, & 166 & slafilesact, slafilespas, slafbfiles, & 167 & sstfiles, sstfbfiles, & 168 & ln_seaice, seaicefiles, & 169 & dobsini, dobsend, n1dint, n2dint, & 170 & nmsshc, mdtcorr, mdtcutoff, & 171 & ln_reysst, ln_ghrsst, reysstname, reysstfmt, & 172 & 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 & cn_sstbias_files ! SST bias input filenames 124 CHARACTER(LEN=128) :: & 125 & cn_altbiasfile ! Altimeter bias input filename 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 LOGICAL, DIMENSION(jpmaxnfiles) :: lmask ! Used for finding number of sstbias files 145 146 REAL(dp) :: rn_dobsini ! Obs window start date YYYYMMDD.HHMMSS 147 REAL(dp) :: rn_dobsend ! Obs window end date YYYYMMDD.HHMMSS 148 REAL(wp), POINTER, DIMENSION(:,:) :: & 149 & zglam1, & ! Model longitudes for profile variable 1 150 & zglam2 ! Model longitudes for profile variable 2 151 REAL(wp), POINTER, DIMENSION(:,:) :: & 152 & zgphi1, & ! Model latitudes for profile variable 1 153 & zgphi2 ! Model latitudes for profile variable 2 154 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 155 & zmask1, & ! Model land/sea mask associated with variable 1 156 & zmask2 ! Model land/sea mask associated with variable 2 157 158 NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla, & 159 & ln_sst, ln_sic, ln_vel3d, & 160 & ln_altbias, ln_nea, ln_grid_global, & 173 161 & ln_grid_search_lookup, & 174 & grid_search_file, grid_search_res, & 175 & ln_grid_global, bias_file, ln_altbias, & 176 & endailyavtypes, ln_s_at_t, ln_profb_ena, & 177 & ln_vel3d, ln_velavcur, velavcurfiles, & 178 & ln_velhrcur, velhrcurfiles, & 179 & ln_velavadcp, velavadcpfiles, & 180 & ln_velhradcp, velhradcpfiles, & 181 & ln_velfb, velfbfiles, ln_velfb_av, & 182 & ln_profb_enatim, ln_ignmis, ln_cl4 183 184 INTEGER :: jprofset 185 INTEGER :: jveloset 186 INTEGER :: jvar 187 INTEGER :: jnumenact 188 INTEGER :: jnumcorio 189 INTEGER :: jnumprofb 190 INTEGER :: jnumslaact 191 INTEGER :: jnumslapas 192 INTEGER :: jnumslafb 193 INTEGER :: jnumsst 194 INTEGER :: jnumsstfb 195 INTEGER :: jnumseaice 196 INTEGER :: jnumvelavcur 197 INTEGER :: jnumvelhrcur 198 INTEGER :: jnumvelavadcp 199 INTEGER :: jnumvelhradcp 200 INTEGER :: jnumvelfb 201 INTEGER :: ji 202 INTEGER :: jset 203 INTEGER :: ios ! Local integer output status for namelist read 204 LOGICAL :: lmask(MaxNumFiles), ll_u3d, ll_v3d 162 & ln_ignmis, ln_s_at_t, ln_sstnight, & 163 & cn_profbfiles, cn_slafbfiles, & 164 & cn_sstfbfiles, cn_sicfbfiles, & 165 & cn_velfbfiles, cn_altbiasfile, & 166 & cn_gridsearchfile, rn_gridsearchres, & 167 & rn_dobsini, rn_dobsend, nn_1dint, nn_2dint, & 168 & nn_msshc, rn_mdtcorr, rn_mdtcutoff, & 169 & nn_profdavtypes, ln_sstbias, cn_sstbias_files 170 171 INTEGER :: jnumsstbias 172 CALL wrk_alloc( jpi, jpj, zglam1 ) 173 CALL wrk_alloc( jpi, jpj, zglam2 ) 174 CALL wrk_alloc( jpi, jpj, zgphi1 ) 175 CALL wrk_alloc( jpi, jpj, zgphi2 ) 176 CALL wrk_alloc( jpi, jpj, jpk, zmask1 ) 177 CALL wrk_alloc( jpi, jpj, jpk, zmask2 ) 205 178 206 179 !----------------------------------------------------------------------- 207 180 ! Read namelist parameters 208 181 !----------------------------------------------------------------------- 209 210 enactfiles(:) = ''211 coriofiles(:) = ''212 profbfiles(:) = ''213 slafilesact(:) = ''214 slafilespas(:) = ''215 slafbfiles(:) = ''216 sstfiles(:) = ''217 sstfbfiles(:) = ''218 seaicefiles(:) = ''219 velcurfiles(:) = ''220 veladcpfiles(:) = ''221 velavcurfiles(:) = ''222 velhrcurfiles(:) = ''223 velavadcpfiles(:) = ''224 velhradcpfiles(:) = ''225 velfbfiles(:) = ''226 velcurfiles(:) = ''227 veladcpfiles(:) = ''228 endailyavtypes(:) = -1229 endailyavtypes(1) = 820230 ln_profb_ena(:) = .FALSE.231 ln_profb_enatim(:) = .TRUE.232 ln_velfb_av(:) = .FALSE.233 ln_ignmis = .FALSE.234 182 235 CALL ini_date( dobsini ) 236 CALL fin_date( dobsend ) 237 238 ! Read Namelist namobs : control observation diagnostics 239 REWIND( numnam_ref ) ! Namelist namobs in reference namelist : Diagnostic: control observation 183 !Initalise all values in namelist arrays 184 ALLOCATE(sstbias_type(jpmaxnfiles)) 185 ! Some namelist arrays need initialising 186 cn_profbfiles(:) = '' 187 cn_slafbfiles(:) = '' 188 cn_sstfbfiles(:) = '' 189 cn_sicfbfiles(:) = '' 190 cn_velfbfiles(:) = '' 191 cn_sstbias_files(:) = '' 192 nn_profdavtypes(:) = -1 193 194 CALL ini_date( rn_dobsini ) 195 CALL fin_date( rn_dobsend ) 196 197 ! Read namelist namobs : control observation diagnostics 198 REWIND( numnam_ref ) ! Namelist namobs in reference namelist 240 199 READ ( numnam_ref, namobs, IOSTAT = ios, ERR = 901) 241 200 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp ) 242 201 243 REWIND( numnam_cfg ) ! Namelist namobs in configuration namelist : Diagnostic: control observation202 REWIND( numnam_cfg ) ! Namelist namobs in configuration namelist 244 203 READ ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 ) 245 204 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp ) 246 205 IF(lwm) WRITE ( numond, namobs ) 247 206 248 ! Count number of files for each type 249 IF (ln_ena) THEN 250 lmask(:) = .FALSE. 251 WHERE (enactfiles(:) /= '') lmask(:) = .TRUE. 252 jnumenact = COUNT(lmask) 253 ENDIF 254 IF (ln_cor) THEN 255 lmask(:) = .FALSE. 256 WHERE (coriofiles(:) /= '') lmask(:) = .TRUE. 257 jnumcorio = COUNT(lmask) 258 ENDIF 259 IF (ln_profb) THEN 260 lmask(:) = .FALSE. 261 WHERE (profbfiles(:) /= '') lmask(:) = .TRUE. 262 jnumprofb = COUNT(lmask) 263 ENDIF 264 IF (ln_sladt) THEN 265 lmask(:) = .FALSE. 266 WHERE (slafilesact(:) /= '') lmask(:) = .TRUE. 267 jnumslaact = COUNT(lmask) 268 lmask(:) = .FALSE. 269 WHERE (slafilespas(:) /= '') lmask(:) = .TRUE. 270 jnumslapas = COUNT(lmask) 271 ENDIF 272 IF (ln_slafb) THEN 273 lmask(:) = .FALSE. 274 WHERE (slafbfiles(:) /= '') lmask(:) = .TRUE. 275 jnumslafb = COUNT(lmask) 276 lmask(:) = .FALSE. 277 ENDIF 278 IF (ln_ghrsst) THEN 279 lmask(:) = .FALSE. 280 WHERE (sstfiles(:) /= '') lmask(:) = .TRUE. 281 jnumsst = COUNT(lmask) 207 IF ( .NOT. ln_diaobs ) THEN 208 IF(lwp) WRITE(numout,cform_war) 209 IF(lwp) WRITE(numout,*)' ln_diaobs is set to false so not calling dia_obs' 210 RETURN 211 ENDIF 212 213 !----------------------------------------------------------------------- 214 ! Set up list of observation types to be used 215 ! and the files associated with each type 216 !----------------------------------------------------------------------- 217 218 nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 219 nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic /) ) 220 221 IF (ln_sstbias) THEN 222 lmask(:) = .FALSE. 223 WHERE (cn_sstbias_files(:) /= '') lmask(:) = .TRUE. 224 jnumsstbias = COUNT(lmask) 225 lmask(:) = .FALSE. 282 226 ENDIF 283 IF (ln_sstfb) THEN 284 lmask(:) = .FALSE. 285 WHERE (sstfbfiles(:) /= '') lmask(:) = .TRUE. 286 jnumsstfb = COUNT(lmask) 287 lmask(:) = .FALSE. 288 ENDIF 289 IF (ln_seaice) THEN 290 lmask(:) = .FALSE. 291 WHERE (seaicefiles(:) /= '') lmask(:) = .TRUE. 292 jnumseaice = COUNT(lmask) 293 ENDIF 294 IF (ln_velavcur) THEN 295 lmask(:) = .FALSE. 296 WHERE (velavcurfiles(:) /= '') lmask(:) = .TRUE. 297 jnumvelavcur = COUNT(lmask) 298 ENDIF 299 IF (ln_velhrcur) THEN 300 lmask(:) = .FALSE. 301 WHERE (velhrcurfiles(:) /= '') lmask(:) = .TRUE. 302 jnumvelhrcur = COUNT(lmask) 303 ENDIF 304 IF (ln_velavadcp) THEN 305 lmask(:) = .FALSE. 306 WHERE (velavadcpfiles(:) /= '') lmask(:) = .TRUE. 307 jnumvelavadcp = COUNT(lmask) 308 ENDIF 309 IF (ln_velhradcp) THEN 310 lmask(:) = .FALSE. 311 WHERE (velhradcpfiles(:) /= '') lmask(:) = .TRUE. 312 jnumvelhradcp = COUNT(lmask) 313 ENDIF 314 IF (ln_velfb) THEN 315 lmask(:) = .FALSE. 316 WHERE (velfbfiles(:) /= '') lmask(:) = .TRUE. 317 jnumvelfb = COUNT(lmask) 318 lmask(:) = .FALSE. 319 ENDIF 320 321 ! Control print 227 228 IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 229 IF(lwp) WRITE(numout,cform_war) 230 IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 231 & ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', & 232 & ' are set to .FALSE. so turning off calls to dia_obs' 233 nwarn = nwarn + 1 234 ln_diaobs = .FALSE. 235 RETURN 236 ENDIF 237 238 IF ( nproftypes > 0 ) THEN 239 240 ALLOCATE( cobstypesprof(nproftypes) ) 241 ALLOCATE( ifilesprof(nproftypes) ) 242 ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 243 244 jtype = 0 245 IF (ln_t3d .OR. ln_s3d) THEN 246 jtype = jtype + 1 247 clproffiles(jtype,:) = cn_profbfiles(:) 248 cobstypesprof(jtype) = 'prof ' 249 ifilesprof(jtype) = 0 250 DO jfile = 1, jpmaxnfiles 251 IF ( trim(clproffiles(jtype,jfile)) /= '' ) & 252 ifilesprof(jtype) = ifilesprof(jtype) + 1 253 END DO 254 ENDIF 255 IF (ln_vel3d) THEN 256 jtype = jtype + 1 257 clproffiles(jtype,:) = cn_velfbfiles(:) 258 cobstypesprof(jtype) = 'vel ' 259 ifilesprof(jtype) = 0 260 DO jfile = 1, jpmaxnfiles 261 IF ( trim(clproffiles(jtype,jfile)) /= '' ) & 262 ifilesprof(jtype) = ifilesprof(jtype) + 1 263 END DO 264 ENDIF 265 266 ENDIF 267 268 IF ( nsurftypes > 0 ) THEN 269 270 ALLOCATE( cobstypessurf(nsurftypes) ) 271 ALLOCATE( ifilessurf(nsurftypes) ) 272 ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 273 274 jtype = 0 275 IF (ln_sla) THEN 276 jtype = jtype + 1 277 clsurffiles(jtype,:) = cn_slafbfiles(:) 278 cobstypessurf(jtype) = 'sla ' 279 ifilessurf(jtype) = 0 280 DO jfile = 1, jpmaxnfiles 281 IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 282 ifilessurf(jtype) = ifilessurf(jtype) + 1 283 END DO 284 ENDIF 285 IF (ln_sst) THEN 286 jtype = jtype + 1 287 clsurffiles(jtype,:) = cn_sstfbfiles(:) 288 cobstypessurf(jtype) = 'sst ' 289 ifilessurf(jtype) = 0 290 DO jfile = 1, jpmaxnfiles 291 IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 292 ifilessurf(jtype) = ifilessurf(jtype) + 1 293 END DO 294 ENDIF 295 #if defined key_lim2 || defined key_lim3 296 IF (ln_sic) THEN 297 jtype = jtype + 1 298 clsurffiles(jtype,:) = cn_sicfbfiles(:) 299 cobstypessurf(jtype) = 'sic ' 300 ifilessurf(jtype) = 0 301 DO jfile = 1, jpmaxnfiles 302 IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 303 ifilessurf(jtype) = ifilessurf(jtype) + 1 304 END DO 305 ENDIF 306 #endif 307 308 ENDIF 309 310 !Write namelist settings to stdout 322 311 IF(lwp) THEN 323 312 WRITE(numout,*) … … 325 314 WRITE(numout,*) '~~~~~~~~~~~~' 326 315 WRITE(numout,*) ' Namelist namobs : set observation diagnostic parameters' 327 WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d 328 WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d 329 WRITE(numout,*) ' Logical switch for ENACT insitu data set ln_ena = ', ln_ena 330 WRITE(numout,*) ' Logical switch for Coriolis insitu data set ln_cor = ', ln_cor 331 WRITE(numout,*) ' Logical switch for feedback insitu data set ln_profb = ', ln_profb 332 WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla 333 WRITE(numout,*) ' Logical switch for AVISO SLA data ln_sladt = ', ln_sladt 334 WRITE(numout,*) ' Logical switch for feedback SLA data ln_slafb = ', ln_slafb 335 WRITE(numout,*) ' Logical switch for SSH observations ln_ssh = ', ln_ssh 336 WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst 337 WRITE(numout,*) ' Logical switch for Reynolds observations ln_reysst = ', ln_reysst 338 WRITE(numout,*) ' Logical switch for GHRSST observations ln_ghrsst = ', ln_ghrsst 339 WRITE(numout,*) ' Logical switch for feedback SST data ln_sstfb = ', ln_sstfb 340 WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight 341 WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss 342 WRITE(numout,*) ' Logical switch for Sea Ice observations ln_seaice = ', ln_seaice 343 WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d 344 WRITE(numout,*) ' Logical switch for velocity daily av. cur. ln_velavcur = ', ln_velavcur 345 WRITE(numout,*) ' Logical switch for velocity high freq. cur. ln_velhrcur = ', ln_velhrcur 346 WRITE(numout,*) ' Logical switch for velocity daily av. ADCP ln_velavadcp = ', ln_velavadcp 347 WRITE(numout,*) ' Logical switch for velocity high freq. ADCP ln_velhradcp = ', ln_velhradcp 348 WRITE(numout,*) ' Logical switch for feedback velocity data ln_velfb = ', ln_velfb 349 WRITE(numout,*) ' Global distribtion of observations ln_grid_global = ',ln_grid_global 350 WRITE(numout,*) & 351 ' Logical switch for obs grid search w/lookup table ln_grid_search_lookup = ',ln_grid_search_lookup 316 WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d 317 WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d 318 WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla 319 WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst 320 WRITE(numout,*) ' Logical switch for Sea Ice observations ln_sic = ', ln_sic 321 WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d 322 WRITE(numout,*) ' Global distribution of observations ln_grid_global = ',ln_grid_global 323 WRITE(numout,*) ' Logical switch for SST bias correction ln_sstbias = ', ln_sstbias 324 WRITE(numout,*) ' Logical switch for obs grid search lookup ln_grid_search_lookup = ',ln_grid_search_lookup 352 325 IF (ln_grid_search_lookup) & 353 WRITE(numout,*) ' Grid search lookup file header grid_search_file = ', grid_search_file 354 IF (ln_ena) THEN 355 DO ji = 1, jnumenact 356 WRITE(numout,'(1X,2A)') ' ENACT input observation file name enactfiles = ', & 357 TRIM(enactfiles(ji)) 326 WRITE(numout,*) ' Grid search lookup file header cn_gridsearchfile = ', cn_gridsearchfile 327 WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS rn_dobsini = ', rn_dobsini 328 WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS rn_dobsend = ', rn_dobsend 329 WRITE(numout,*) ' Type of vertical interpolation method nn_1dint = ', nn_1dint 330 WRITE(numout,*) ' Type of horizontal interpolation method nn_2dint = ', nn_2dint 331 WRITE(numout,*) ' Rejection of observations near land switch ln_nea = ', ln_nea 332 WRITE(numout,*) ' MSSH correction scheme nn_msshc = ', nn_msshc 333 WRITE(numout,*) ' MDT correction rn_mdtcorr = ', rn_mdtcorr 334 WRITE(numout,*) ' MDT cutoff for computed correction rn_mdtcutoff = ', rn_mdtcutoff 335 WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias 336 WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis 337 WRITE(numout,*) ' Daily average types nn_profdavtypes = ', nn_profdavtypes 338 WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight 339 WRITE(numout,*) ' Number of profile obs types: ',nproftypes 340 341 IF ( nproftypes > 0 ) THEN 342 DO jtype = 1, nproftypes 343 DO jfile = 1, ifilesprof(jtype) 344 WRITE(numout,'(1X,2A)') ' '//cobstypesprof(jtype)//' input observation file names = ', & 345 TRIM(clproffiles(jtype,jfile)) 346 END DO 358 347 END DO 359 348 ENDIF 360 IF (ln_cor) THEN 361 DO ji = 1, jnumcorio 362 WRITE(numout,'(1X,2A)') ' Coriolis input observation file name coriofiles = ', & 363 TRIM(coriofiles(ji)) 349 350 WRITE(numout,*)' Number of surface obs types: ',nsurftypes 351 IF ( nsurftypes > 0 ) THEN 352 DO jtype = 1, nsurftypes 353 DO jfile = 1, ifilessurf(jtype) 354 WRITE(numout,'(1X,2A)') ' '//cobstypessurf(jtype)//' input observation file names = ', & 355 TRIM(clsurffiles(jtype,jfile)) 356 END DO 364 357 END DO 365 358 ENDIF 366 IF (ln_profb) THEN 367 DO ji = 1, jnumprofb 368 IF (ln_profb_ena(ji)) THEN 369 WRITE(numout,'(1X,2A)') ' Enact feedback input observation file name profbfiles = ', & 370 TRIM(profbfiles(ji)) 371 ELSE 372 WRITE(numout,'(1X,2A)') ' Feedback input observation file name profbfiles = ', & 373 TRIM(profbfiles(ji)) 374 ENDIF 375 WRITE(numout,'(1X,2A)') ' Enact feedback input time setting switch ln_profb_enatim = ', ln_profb_enatim(ji) 376 END DO 377 ENDIF 378 IF (ln_sladt) THEN 379 DO ji = 1, jnumslaact 380 WRITE(numout,'(1X,2A)') ' Active SLA input observation file name slafilesact = ', & 381 TRIM(slafilesact(ji)) 382 END DO 383 DO ji = 1, jnumslapas 384 WRITE(numout,'(1X,2A)') ' Passive SLA input observation file name slafilespas = ', & 385 TRIM(slafilespas(ji)) 386 END DO 387 ENDIF 388 IF (ln_slafb) THEN 389 DO ji = 1, jnumslafb 390 WRITE(numout,'(1X,2A)') ' Feedback SLA input observation file name slafbfiles = ', & 391 TRIM(slafbfiles(ji)) 392 END DO 393 ENDIF 394 IF (ln_ghrsst) THEN 395 DO ji = 1, jnumsst 396 WRITE(numout,'(1X,2A)') ' GHRSST input observation file name sstfiles = ', & 397 TRIM(sstfiles(ji)) 398 END DO 399 ENDIF 400 IF (ln_sstfb) THEN 401 DO ji = 1, jnumsstfb 402 WRITE(numout,'(1X,2A)') ' Feedback SST input observation file name sstfbfiles = ', & 403 TRIM(sstfbfiles(ji)) 404 END DO 405 ENDIF 406 IF (ln_seaice) THEN 407 DO ji = 1, jnumseaice 408 WRITE(numout,'(1X,2A)') ' Sea Ice input observation file name seaicefiles = ', & 409 TRIM(seaicefiles(ji)) 410 END DO 411 ENDIF 412 IF (ln_velavcur) THEN 413 DO ji = 1, jnumvelavcur 414 WRITE(numout,'(1X,2A)') ' Vel. cur. daily av. input file name velavcurfiles = ', & 415 TRIM(velavcurfiles(ji)) 416 END DO 417 ENDIF 418 IF (ln_velhrcur) THEN 419 DO ji = 1, jnumvelhrcur 420 WRITE(numout,'(1X,2A)') ' Vel. cur. high freq. input file name velhvcurfiles = ', & 421 TRIM(velhrcurfiles(ji)) 422 END DO 423 ENDIF 424 IF (ln_velavadcp) THEN 425 DO ji = 1, jnumvelavadcp 426 WRITE(numout,'(1X,2A)') ' Vel. ADCP daily av. input file name velavadcpfiles = ', & 427 TRIM(velavadcpfiles(ji)) 428 END DO 429 ENDIF 430 IF (ln_velhradcp) THEN 431 DO ji = 1, jnumvelhradcp 432 WRITE(numout,'(1X,2A)') ' Vel. ADCP high freq. input file name velhvadcpfiles = ', & 433 TRIM(velhradcpfiles(ji)) 434 END DO 435 ENDIF 436 IF (ln_velfb) THEN 437 DO ji = 1, jnumvelfb 438 IF (ln_velfb_av(ji)) THEN 439 WRITE(numout,'(1X,2A)') ' Vel. feedback daily av. input file name velfbfiles = ', & 440 TRIM(velfbfiles(ji)) 441 ELSE 442 WRITE(numout,'(1X,2A)') ' Vel. feedback input observation file name velfbfiles = ', & 443 TRIM(velfbfiles(ji)) 444 ENDIF 445 END DO 446 ENDIF 447 WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS dobsini = ', dobsini 448 WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS dobsend = ', dobsend 449 WRITE(numout,*) ' Type of vertical interpolation method n1dint = ', n1dint 450 WRITE(numout,*) ' Type of horizontal interpolation method n2dint = ', n2dint 451 WRITE(numout,*) ' Rejection of observations near land swithch ln_nea = ', ln_nea 452 WRITE(numout,*) ' MSSH correction scheme nmsshc = ', nmsshc 453 WRITE(numout,*) ' MDT correction mdtcorr = ', mdtcorr 454 WRITE(numout,*) ' MDT cutoff for computed correction mdtcutoff = ', mdtcutoff 455 WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias 456 WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis 457 WRITE(numout,*) ' ENACT daily average types = ',endailyavtypes 458 459 ENDIF 460 359 WRITE(numout,*) '~~~~~~~~~~~~' 360 361 ENDIF 362 363 !----------------------------------------------------------------------- 364 ! Obs operator parameter checking and initialisations 365 !----------------------------------------------------------------------- 366 461 367 IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN 462 368 CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' ) … … 464 370 ENDIF 465 371 466 CALL obs_typ_init 467 468 CALL mppmap_init 469 470 ! Parameter control 471 #if defined key_diaobs 472 IF ( ( .NOT. ln_t3d ).AND.( .NOT. ln_s3d ).AND.( .NOT. ln_sla ).AND. & 473 & ( .NOT. ln_vel3d ).AND. & 474 & ( .NOT. ln_ssh ).AND.( .NOT. ln_sst ).AND.( .NOT. ln_sss ).AND. & 475 & ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ) ) THEN 476 IF(lwp) WRITE(numout,cform_war) 477 IF(lwp) WRITE(numout,*) ' key_diaobs is activated but logical flags', & 478 & ' ln_t3d, ln_s3d, ln_sla, ln_ssh, ln_sst, ln_sss, ln_seaice, ln_vel3d are all set to .FALSE.' 479 nwarn = nwarn + 1 480 ENDIF 481 #endif 482 483 CALL obs_grid_setup( ) 484 IF ( ( n1dint < 0 ).OR.( n1dint > 1 ) ) THEN 372 IF ( ln_grid_global ) THEN 373 CALL ctl_warn( 'ln_grid_global=T may cause memory issues when used with a large number of processors' ) 374 ENDIF 375 376 IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN 485 377 CALL ctl_stop(' Choice of vertical (1D) interpolation method', & 486 378 & ' is not available') 487 379 ENDIF 488 IF ( ( n2dint < 0 ).OR.( n2dint > 4 ) ) THEN 380 381 IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 4 ) ) THEN 489 382 CALL ctl_stop(' Choice of horizontal (2D) interpolation method', & 490 383 & ' is not available') 491 384 ENDIF 492 385 386 CALL obs_typ_init 387 IF(ln_grid_global) THEN 388 CALL mppmap_init 389 ENDIF 390 391 CALL obs_grid_setup( ) 392 493 393 !----------------------------------------------------------------------- 494 394 ! Depending on switches read the various observation types 495 395 !----------------------------------------------------------------------- 496 ! - Temperature/salinity profiles 497 498 IF ( ln_t3d .OR. ln_s3d ) THEN 499 500 ! Set the number of variables for profiles to 2 (T and S) 501 nprofvars = 2 502 ! Set the number of extra variables for profiles to 1 (insitu temp). 503 nprofextr = 1 504 505 ! Count how may insitu data sets we have and allocate data. 506 jprofset = 0 507 IF ( ln_ena ) jprofset = jprofset + 1 508 IF ( ln_cor ) jprofset = jprofset + 1 509 IF ( ln_profb ) jprofset = jprofset + jnumprofb 510 nprofsets = jprofset 511 IF ( nprofsets > 0 ) THEN 512 ALLOCATE(ld_enact(nprofsets)) 513 ALLOCATE(profdata(nprofsets)) 514 ALLOCATE(prodatqc(nprofsets)) 515 ENDIF 516 517 jprofset = 0 518 519 ! ENACT insitu data 520 521 IF ( ln_ena ) THEN 522 523 jprofset = jprofset + 1 524 525 ld_enact(jprofset) = .TRUE. 526 527 CALL obs_rea_pro_dri( 1, profdata(jprofset), & 528 & jnumenact, enactfiles(1:jnumenact), & 529 & nprofvars, nprofextr, & 530 & nitend-nit000+2, & 531 & dobsini, dobsend, ln_t3d, ln_s3d, & 532 & ln_ignmis, ln_s_at_t, .TRUE., .FALSE., & 533 & kdailyavtypes = endailyavtypes ) 534 535 DO jvar = 1, 2 536 537 CALL obs_prof_staend( profdata(jprofset), jvar ) 538 396 397 IF ( nproftypes > 0 ) THEN 398 399 ALLOCATE(profdata(nproftypes)) 400 ALLOCATE(profdataqc(nproftypes)) 401 ALLOCATE(nvarsprof(nproftypes)) 402 ALLOCATE(nextrprof(nproftypes)) 403 404 DO jtype = 1, nproftypes 405 406 nvarsprof(jtype) = 2 407 IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN 408 nextrprof(jtype) = 1 409 llvar1 = ln_t3d 410 llvar2 = ln_s3d 411 zglam1 = glamt 412 zgphi1 = gphit 413 zmask1 = tmask 414 zglam2 = glamt 415 zgphi2 = gphit 416 zmask2 = tmask 417 ENDIF 418 IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN 419 nextrprof(jtype) = 2 420 llvar1 = ln_vel3d 421 llvar2 = ln_vel3d 422 zglam1 = glamu 423 zgphi1 = gphiu 424 zmask1 = umask 425 zglam2 = glamv 426 zgphi2 = gphiv 427 zmask2 = vmask 428 ENDIF 429 430 !Read in profile or profile obs types 431 CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype), & 432 & clproffiles(jtype,1:ifilesprof(jtype)), & 433 & nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, & 434 & rn_dobsini, rn_dobsend, llvar1, llvar2, & 435 & ln_ignmis, ln_s_at_t, .FALSE., & 436 & kdailyavtypes = nn_profdavtypes ) 437 438 DO jvar = 1, nvarsprof(jtype) 439 CALL obs_prof_staend( profdata(jtype), jvar ) 539 440 END DO 540 441 541 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 542 & ln_t3d, ln_s3d, ln_nea, & 543 & kdailyavtypes=endailyavtypes ) 544 545 ENDIF 546 547 ! Coriolis insitu data 548 549 IF ( ln_cor ) THEN 550 551 jprofset = jprofset + 1 552 553 ld_enact(jprofset) = .FALSE. 554 555 CALL obs_rea_pro_dri( 2, profdata(jprofset), & 556 & jnumcorio, coriofiles(1:jnumcorio), & 557 & nprofvars, nprofextr, & 558 & nitend-nit000+2, & 559 & dobsini, dobsend, ln_t3d, ln_s3d, & 560 & ln_ignmis, ln_s_at_t, .FALSE., .FALSE. ) 561 562 DO jvar = 1, 2 563 564 CALL obs_prof_staend( profdata(jprofset), jvar ) 565 566 END DO 567 568 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 569 & ln_t3d, ln_s3d, ln_nea ) 570 571 ENDIF 572 573 ! Feedback insitu data 574 575 IF ( ln_profb ) THEN 576 577 DO jset = 1, jnumprofb 578 579 jprofset = jprofset + 1 580 ld_enact (jprofset) = ln_profb_ena(jset) 581 582 CALL obs_rea_pro_dri( 0, profdata(jprofset), & 583 & 1, profbfiles(jset:jset), & 584 & nprofvars, nprofextr, & 585 & nitend-nit000+2, & 586 & dobsini, dobsend, ln_t3d, ln_s3d, & 587 & ln_ignmis, ln_s_at_t, & 588 & ld_enact(jprofset).AND.& 589 & ln_profb_enatim(jset), & 590 & .FALSE., kdailyavtypes = endailyavtypes ) 591 592 DO jvar = 1, 2 593 594 CALL obs_prof_staend( profdata(jprofset), jvar ) 595 596 END DO 597 598 IF ( ld_enact(jprofset) ) THEN 599 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 600 & ln_t3d, ln_s3d, ln_nea, & 601 & kdailyavtypes = endailyavtypes ) 602 ELSE 603 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 604 & ln_t3d, ln_s3d, ln_nea ) 605 ENDIF 606 607 END DO 608 609 ENDIF 610 611 ENDIF 612 613 ! - Sea level anomalies 614 IF ( ln_sla ) THEN 615 ! Set the number of variables for sla to 1 616 nslavars = 1 617 618 ! Set the number of extra variables for sla to 2 619 nslaextr = 2 442 CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), & 443 & llvar1, llvar2, & 444 & jpi, jpj, jpk, & 445 & zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2, & 446 & ln_nea, kdailyavtypes = nn_profdavtypes ) 447 448 END DO 449 450 DEALLOCATE( ifilesprof, clproffiles ) 451 452 ENDIF 453 454 IF ( nsurftypes > 0 ) THEN 455 456 ALLOCATE(surfdata(nsurftypes)) 457 ALLOCATE(surfdataqc(nsurftypes)) 458 ALLOCATE(nvarssurf(nsurftypes)) 459 ALLOCATE(nextrsurf(nsurftypes)) 460 461 DO jtype = 1, nsurftypes 462 463 nvarssurf(jtype) = 1 464 nextrsurf(jtype) = 0 465 llnightav = .FALSE. 466 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2 467 IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav = ln_sstnight 468 469 !Read in surface obs types 470 CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), & 471 & clsurffiles(jtype,1:ifilessurf(jtype)), & 472 & nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 473 & rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav ) 620 474 621 ! Set the number of sla data sets to 2622 nslasets = 0623 IF ( ln_sladt ) THEN624 nslasets = nslasets + 2625 ENDIF626 IF ( ln_slafb ) THEN627 nslasets = nslasets + jnumslafb628 ENDIF629 475 630 ALLOCATE(sladata(nslasets)) 631 ALLOCATE(sladatqc(nslasets)) 632 sladata(:)%nsurf=0 633 sladatqc(:)%nsurf=0 634 635 nslasets = 0 636 637 ! AVISO SLA data 638 639 IF ( ln_sladt ) THEN 640 641 ! Active SLA observations 642 643 nslasets = nslasets + 1 644 645 CALL obs_rea_sla( 1, sladata(nslasets), jnumslaact, & 646 & slafilesact(1:jnumslaact), & 647 & nslavars, nslaextr, nitend-nit000+2, & 648 & dobsini, dobsend, ln_ignmis, .FALSE. ) 649 CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), & 650 & ln_sla, ln_nea ) 651 652 ! Passive SLA observations 653 654 nslasets = nslasets + 1 655 656 CALL obs_rea_sla( 1, sladata(nslasets), jnumslapas, & 657 & slafilespas(1:jnumslapas), & 658 & nslavars, nslaextr, nitend-nit000+2, & 659 & dobsini, dobsend, ln_ignmis, .FALSE. ) 660 661 CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), & 662 & ln_sla, ln_nea ) 663 664 ENDIF 665 666 ! Feedback SLA data 667 668 IF ( ln_slafb ) THEN 669 670 DO jset = 1, jnumslafb 671 672 nslasets = nslasets + 1 673 674 CALL obs_rea_sla( 0, sladata(nslasets), 1, & 675 & slafbfiles(jset:jset), & 676 & nslavars, nslaextr, nitend-nit000+2, & 677 & dobsini, dobsend, ln_ignmis, .FALSE. ) 678 CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), & 679 & ln_sla, ln_nea ) 680 681 END DO 682 683 ENDIF 684 685 CALL obs_rea_mdt( nslasets, sladatqc, n2dint ) 686 687 ! read in altimeter bias 688 689 IF ( ln_altbias ) THEN 690 CALL obs_rea_altbias ( nslasets, sladatqc, n2dint, bias_file ) 691 ENDIF 692 693 ENDIF 694 695 ! - Sea surface height 696 IF ( ln_ssh ) THEN 697 IF(lwp) WRITE(numout,*) ' SSH currently not available' 698 ENDIF 699 700 ! - Sea surface temperature 701 IF ( ln_sst ) THEN 702 703 ! Set the number of variables for sst to 1 704 nsstvars = 1 705 706 ! Set the number of extra variables for sst to 0 707 nsstextr = 0 708 709 nsstsets = 0 710 711 IF (ln_reysst) nsstsets = nsstsets + 1 712 IF (ln_ghrsst) nsstsets = nsstsets + 1 713 IF ( ln_sstfb ) THEN 714 nsstsets = nsstsets + jnumsstfb 715 ENDIF 716 717 ALLOCATE(sstdata(nsstsets)) 718 ALLOCATE(sstdatqc(nsstsets)) 719 ALLOCATE(ld_sstnight(nsstsets)) 720 sstdata(:)%nsurf=0 721 sstdatqc(:)%nsurf=0 722 ld_sstnight(:)=.false. 723 724 nsstsets = 0 725 726 IF (ln_reysst) THEN 727 728 nsstsets = nsstsets + 1 729 730 ld_sstnight(nsstsets) = ln_sstnight 731 732 CALL obs_rea_sst_rey( reysstname, reysstfmt, sstdata(nsstsets), & 733 & nsstvars, nsstextr, & 734 & nitend-nit000+2, dobsini, dobsend ) 735 CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, & 736 & ln_nea ) 737 738 ENDIF 739 740 IF (ln_ghrsst) THEN 741 742 nsstsets = nsstsets + 1 743 744 ld_sstnight(nsstsets) = ln_sstnight 745 746 CALL obs_rea_sst( 1, sstdata(nsstsets), jnumsst, & 747 & sstfiles(1:jnumsst), & 748 & nsstvars, nsstextr, nitend-nit000+2, & 749 & dobsini, dobsend, ln_ignmis, .FALSE. ) 750 CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, & 751 & ln_nea ) 752 753 ENDIF 754 755 ! Feedback SST data 756 757 IF ( ln_sstfb ) THEN 758 759 DO jset = 1, jnumsstfb 760 761 nsstsets = nsstsets + 1 762 763 ld_sstnight(nsstsets) = ln_sstnight 764 765 CALL obs_rea_sst( 0, sstdata(nsstsets), 1, & 766 & sstfbfiles(jset:jset), & 767 & nsstvars, nsstextr, nitend-nit000+2, & 768 & dobsini, dobsend, ln_ignmis, .FALSE. ) 769 CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), & 770 & ln_sst, ln_nea ) 771 772 END DO 773 774 ENDIF 775 776 ENDIF 777 778 ! - Sea surface salinity 779 IF ( ln_sss ) THEN 780 IF(lwp) WRITE(numout,*) ' SSS currently not available' 781 ENDIF 782 783 ! - Sea Ice Concentration 784 785 IF ( ln_seaice ) THEN 786 787 ! Set the number of variables for seaice to 1 788 nseaicevars = 1 789 790 ! Set the number of extra variables for seaice to 0 791 nseaiceextr = 0 792 793 ! Set the number of data sets to 1 794 nseaicesets = 1 795 796 ALLOCATE(seaicedata(nseaicesets)) 797 ALLOCATE(seaicedatqc(nseaicesets)) 798 seaicedata(:)%nsurf=0 799 seaicedatqc(:)%nsurf=0 800 801 CALL obs_rea_seaice( 1, seaicedata(nseaicesets), jnumseaice, & 802 & seaicefiles(1:jnumseaice), & 803 & nseaicevars, nseaiceextr, nitend-nit000+2, & 804 & dobsini, dobsend, ln_ignmis, .FALSE. ) 805 806 CALL obs_pre_seaice( seaicedata(nseaicesets), seaicedatqc(nseaicesets), & 807 & ln_seaice, ln_nea ) 808 809 ENDIF 810 811 IF (ln_vel3d) THEN 812 813 ! Set the number of variables for profiles to 2 (U and V) 814 nvelovars = 2 815 816 ! Set the number of extra variables for profiles to 2 to store 817 ! rotation parameters 818 nveloextr = 2 819 820 jveloset = 0 821 822 IF ( ln_velavcur ) jveloset = jveloset + 1 823 IF ( ln_velhrcur ) jveloset = jveloset + 1 824 IF ( ln_velavadcp ) jveloset = jveloset + 1 825 IF ( ln_velhradcp ) jveloset = jveloset + 1 826 IF (ln_velfb) jveloset = jveloset + jnumvelfb 827 828 nvelosets = jveloset 829 IF ( nvelosets > 0 ) THEN 830 ALLOCATE( velodata(nvelosets) ) 831 ALLOCATE( veldatqc(nvelosets) ) 832 ALLOCATE( ld_velav(nvelosets) ) 833 ENDIF 834 835 jveloset = 0 836 837 ! Daily averaged data 838 839 IF ( ln_velavcur ) THEN 840 841 jveloset = jveloset + 1 842 843 ld_velav(jveloset) = .TRUE. 844 845 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavcur, & 846 & velavcurfiles(1:jnumvelavcur), & 847 & nvelovars, nveloextr, & 848 & nitend-nit000+2, & 849 & dobsini, dobsend, ln_ignmis, & 850 & ld_velav(jveloset), & 851 & .FALSE. ) 852 853 DO jvar = 1, 2 854 CALL obs_prof_staend( velodata(jveloset), jvar ) 855 END DO 856 857 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 858 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 859 860 ENDIF 861 862 ! High frequency data 863 864 IF ( ln_velhrcur ) THEN 865 866 jveloset = jveloset + 1 867 868 ld_velav(jveloset) = .FALSE. 869 870 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhrcur, & 871 & velhrcurfiles(1:jnumvelhrcur), & 872 & nvelovars, nveloextr, & 873 & nitend-nit000+2, & 874 & dobsini, dobsend, ln_ignmis, & 875 & ld_velav(jveloset), & 876 & .FALSE. ) 877 878 DO jvar = 1, 2 879 CALL obs_prof_staend( velodata(jveloset), jvar ) 880 END DO 881 882 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 883 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 884 885 ENDIF 886 887 ! Daily averaged data 888 889 IF ( ln_velavadcp ) THEN 890 891 jveloset = jveloset + 1 892 893 ld_velav(jveloset) = .TRUE. 894 895 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavadcp, & 896 & velavadcpfiles(1:jnumvelavadcp), & 897 & nvelovars, nveloextr, & 898 & nitend-nit000+2, & 899 & dobsini, dobsend, ln_ignmis, & 900 & ld_velav(jveloset), & 901 & .FALSE. ) 902 903 DO jvar = 1, 2 904 CALL obs_prof_staend( velodata(jveloset), jvar ) 905 END DO 906 907 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 908 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 909 910 ENDIF 911 912 ! High frequency data 913 914 IF ( ln_velhradcp ) THEN 915 916 jveloset = jveloset + 1 917 918 ld_velav(jveloset) = .FALSE. 919 920 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhradcp, & 921 & velhradcpfiles(1:jnumvelhradcp), & 922 & nvelovars, nveloextr, & 923 & nitend-nit000+2, & 924 & dobsini, dobsend, ln_ignmis, & 925 & ld_velav(jveloset), & 926 & .FALSE. ) 927 928 DO jvar = 1, 2 929 CALL obs_prof_staend( velodata(jveloset), jvar ) 930 END DO 931 932 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 933 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 934 935 ENDIF 936 937 IF ( ln_velfb ) THEN 938 939 DO jset = 1, jnumvelfb 940 941 jveloset = jveloset + 1 942 943 ld_velav(jveloset) = ln_velfb_av(jset) 944 945 CALL obs_rea_vel_dri( 0, velodata(jveloset), 1, & 946 & velfbfiles(jset:jset), & 947 & nvelovars, nveloextr, & 948 & nitend-nit000+2, & 949 & dobsini, dobsend, ln_ignmis, & 950 & ld_velav(jveloset), & 951 & .FALSE. ) 952 953 DO jvar = 1, 2 954 CALL obs_prof_staend( velodata(jveloset), jvar ) 955 END DO 956 957 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 958 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 959 960 961 END DO 962 963 ENDIF 964 965 ENDIF 966 476 CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea ) 477 478 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 479 CALL obs_rea_mdt( surfdataqc(jtype), nn_2dint ) 480 IF ( ln_altbias ) CALL obs_rea_altbias ( surfdataqc(jtype), nn_2dint, cn_altbiasfile ) 481 ENDIF 482 IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN 483 !Read in bias field and correct SST. 484 IF ( jnumsstbias == 0 ) CALL ctl_stop("ln_sstbias set,"// & 485 " but no bias"// & 486 " files to read in") 487 CALL obs_app_sstbias( surfdataqc(jtype), nn_2dint, & 488 jnumsstbias, cn_sstbias_files(1:jnumsstbias) ) 489 ENDIF 490 END DO 491 492 DEALLOCATE( ifilessurf, clsurffiles ) 493 494 ENDIF 495 496 CALL wrk_dealloc( jpi, jpj, zglam1 ) 497 CALL wrk_dealloc( jpi, jpj, zglam2 ) 498 CALL wrk_dealloc( jpi, jpj, zgphi1 ) 499 CALL wrk_dealloc( jpi, jpj, zgphi2 ) 500 CALL wrk_dealloc( jpi, jpj, jpk, zmask1 ) 501 CALL wrk_dealloc( jpi, jpj, jpk, zmask2 ) 502 967 503 END SUBROUTINE dia_obs_init 968 504 … … 974 510 !! 975 511 !! ** Method : Call the observation operators on each time step to 976 !! compute the model equivalent of the following date: 977 !! - T profiles 978 !! - S profiles 979 !! - Sea surface height (referenced to a mean) 980 !! - Sea surface temperature 981 !! - Sea surface salinity 982 !! - Velocity component (U,V) profiles 983 !! 984 !! ** Action : 512 !! compute the model equivalent of the following data: 513 !! - Profile data, currently T/S or U/V 514 !! - Surface data, currently SST, SLA or sea-ice concentration. 515 !! 516 !! ** Action : 985 517 !! 986 518 !! History : … … 991 523 !! ! 07-04 (G. Smith) Generalized surface operators 992 524 !! ! 08-10 (M. Valdivieso) obs operator for velocity profiles 525 !! ! 14-08 (J. While) observation operator for profiles in 526 !! generalised vertical coordinates 527 !! ! 15-08 (M. Martin) Combined surface/profile routines. 993 528 !!---------------------------------------------------------------------- 994 529 !! * Modules used 995 530 USE dom_oce, ONLY : & ! Ocean space and time domain variables 996 & rdt, & 997 & gdept_1d, & 998 & tmask, umask, vmask 531 #if defined key_vvl 532 & gdept_n 533 #else 534 & gdept_1d 535 #endif 999 536 USE phycst, ONLY : & ! Physical constants 1000 537 & rday 1001 538 USE oce, ONLY : & ! Ocean dynamics and tracers variables 1002 539 & tsn, & 1003 & un, vn, & 1004 & sshn 540 & un, vn, & 541 & sshn 542 USE phycst, ONLY : & ! Physical constants 543 & rday 1005 544 #if defined key_lim3 1006 USE ice, ONLY : & ! LIMIce model variables545 USE ice, ONLY : & ! LIM3 Ice model variables 1007 546 & frld 1008 547 #endif 1009 548 #if defined key_lim2 1010 USE ice_2, ONLY : & ! LIMIce model variables549 USE ice_2, ONLY : & ! LIM2 Ice model variables 1011 550 & frld 1012 551 #endif … … 1014 553 1015 554 !! * Arguments 1016 INTEGER, INTENT(IN) :: kstp 555 INTEGER, INTENT(IN) :: kstp ! Current timestep 1017 556 !! * Local declarations 1018 INTEGER :: idaystp ! Number of timesteps per day 1019 INTEGER :: jprofset ! Profile data set loop variable 1020 INTEGER :: jslaset ! SLA data set loop variable 1021 INTEGER :: jsstset ! SST data set loop variable 1022 INTEGER :: jseaiceset ! sea ice data set loop variable 1023 INTEGER :: jveloset ! velocity profile data loop variable 1024 INTEGER :: jvar ! Variable number 557 INTEGER :: idaystp ! Number of timesteps per day 558 INTEGER :: jtype ! Data loop variable 559 INTEGER :: jvar ! Variable number 560 INTEGER :: ji, jj ! Loop counters 561 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 562 & zprofvar1, & ! Model values for 1st variable in a prof ob 563 & zprofvar2 ! Model values for 2nd variable in a prof ob 564 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 565 & zprofmask1, & ! Mask associated with zprofvar1 566 & zprofmask2 ! Mask associated with zprofvar2 567 REAL(wp), POINTER, DIMENSION(:,:) :: & 568 & zsurfvar ! Model values equivalent to surface ob. 569 REAL(wp), POINTER, DIMENSION(:,:) :: & 570 & zglam1, & ! Model longitudes for prof variable 1 571 & zglam2, & ! Model longitudes for prof variable 2 572 & zgphi1, & ! Model latitudes for prof variable 1 573 & zgphi2 ! Model latitudes for prof variable 2 1025 574 #if ! defined key_lim2 && ! defined key_lim3 1026 REAL(wp), POINTER, DIMENSION(:,:) :: frld 575 REAL(wp), POINTER, DIMENSION(:,:) :: frld 1027 576 #endif 1028 CHARACTER(LEN=20) :: datestr=" ",timestr=" " 1029 577 LOGICAL :: llnightav ! Logical for calculating night-time average 578 579 !Allocate local work arrays 580 CALL wrk_alloc( jpi, jpj, jpk, zprofvar1 ) 581 CALL wrk_alloc( jpi, jpj, jpk, zprofvar2 ) 582 CALL wrk_alloc( jpi, jpj, jpk, zprofmask1 ) 583 CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 ) 584 CALL wrk_alloc( jpi, jpj, zsurfvar ) 585 CALL wrk_alloc( jpi, jpj, zglam1 ) 586 CALL wrk_alloc( jpi, jpj, zglam2 ) 587 CALL wrk_alloc( jpi, jpj, zgphi1 ) 588 CALL wrk_alloc( jpi, jpj, zgphi2 ) 1030 589 #if ! defined key_lim2 && ! defined key_lim3 1031 590 CALL wrk_alloc(jpi,jpj,frld) … … 1047 606 #endif 1048 607 !----------------------------------------------------------------------- 1049 ! Depending on switches call various observation operators 1050 !----------------------------------------------------------------------- 1051 1052 ! - Temperature/salinity profiles 1053 IF ( ln_t3d .OR. ln_s3d ) THEN 1054 DO jprofset = 1, nprofsets 1055 IF ( ld_enact(jprofset) ) THEN 1056 CALL obs_pro_opt( prodatqc(jprofset), & 1057 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1058 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1059 & gdept_1d, tmask, n1dint, n2dint, & 1060 & kdailyavtypes = endailyavtypes ) 608 ! Call the profile and surface observation operators 609 !----------------------------------------------------------------------- 610 611 IF ( nproftypes > 0 ) THEN 612 613 DO jtype = 1, nproftypes 614 615 SELECT CASE ( TRIM(cobstypesprof(jtype)) ) 616 CASE('prof') 617 zprofvar1(:,:,:) = tsn(:,:,:,jp_tem) 618 zprofvar2(:,:,:) = tsn(:,:,:,jp_sal) 619 zprofmask1(:,:,:) = tmask(:,:,:) 620 zprofmask2(:,:,:) = tmask(:,:,:) 621 zglam1(:,:) = glamt(:,:) 622 zglam2(:,:) = glamt(:,:) 623 zgphi1(:,:) = gphit(:,:) 624 zgphi2(:,:) = gphit(:,:) 625 CASE('vel') 626 zprofvar1(:,:,:) = un(:,:,:) 627 zprofvar2(:,:,:) = vn(:,:,:) 628 zprofmask1(:,:,:) = umask(:,:,:) 629 zprofmask2(:,:,:) = vmask(:,:,:) 630 zglam1(:,:) = glamu(:,:) 631 zglam2(:,:) = glamv(:,:) 632 zgphi1(:,:) = gphiu(:,:) 633 zgphi2(:,:) = gphiv(:,:) 634 END SELECT 635 636 IF( ln_zco .OR. ln_zps ) THEN 637 CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk, & 638 & nit000, idaystp, & 639 & zprofvar1, zprofvar2, & 640 & gdept_1d, zprofmask1, zprofmask2, & 641 & zglam1, zglam2, zgphi1, zgphi2, & 642 & nn_1dint, nn_2dint, & 643 & kdailyavtypes = nn_profdavtypes ) 644 ELSE IF(TRIM(cobstypesprof(jtype)) == 'prof') THEN 645 !TG - THIS NEEDS MODIFICATION TO MATCH SIMPLIFICATION 646 CALL obs_pro_sco_opt( profdataqc(jtype), & 647 & kstp, jpi, jpj, jpk, nit000, idaystp, & 648 & zprofvar1, zprofvar2, & 649 & fsdept(:,:,:), fsdepw(:,:,:), & 650 & tmask, nn_1dint, nn_2dint, & 651 & kdailyavtypes = nn_profdavtypes ) 1061 652 ELSE 1062 CALL obs_pro_opt( prodatqc(jprofset), & 1063 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1064 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1065 & gdept_1d, tmask, n1dint, n2dint ) 653 CALL ctl_stop('DIA_OBS: Generalised vertical interpolation not'// & 654 'yet working for velocity data (turn off velocity observations') 1066 655 ENDIF 656 1067 657 END DO 1068 ENDIF 1069 1070 ! - Sea surface anomaly 1071 IF ( ln_sla ) THEN 1072 DO jslaset = 1, nslasets 1073 CALL obs_sla_opt( sladatqc(jslaset), & 1074 & kstp, jpi, jpj, nit000, sshn, & 1075 & tmask(:,:,1), n2dint ) 1076 END DO 1077 ENDIF 1078 1079 ! - Sea surface temperature 1080 IF ( ln_sst ) THEN 1081 DO jsstset = 1, nsstsets 1082 CALL obs_sst_opt( sstdatqc(jsstset), & 1083 & kstp, jpi, jpj, nit000, idaystp, & 1084 & tsn(:,:,1,jp_tem), tmask(:,:,1), & 1085 & n2dint, ld_sstnight(jsstset) ) 658 659 ENDIF 660 661 IF ( nsurftypes > 0 ) THEN 662 663 DO jtype = 1, nsurftypes 664 665 SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 666 CASE('sst') 667 zsurfvar(:,:) = tsn(:,:,1,jp_tem) 668 llnightav = ln_sstnight 669 CASE('sla') 670 zsurfvar(:,:) = sshn(:,:) 671 llnightav = .FALSE. 672 #if defined key_lim2 || defined key_lim3 673 CASE('sic') 674 IF ( kstp == 0 ) THEN 675 IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN 676 CALL ctl_warn( 'Sea-ice not initialised on zeroth '// & 677 & 'time-step but some obs are valid then.' ) 678 WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & 679 & ' sea-ice obs will be missed' 680 ENDIF 681 surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + & 682 & surfdataqc(jtype)%nsstp(1) 683 CYCLE 684 ELSE 685 zsurfvar(:,:) = 1._wp - frld(:,:) 686 ENDIF 687 688 llnightav = .FALSE. 689 #endif 690 END SELECT 691 692 CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj, & 693 & nit000, idaystp, zsurfvar, tmask(:,:,1), & 694 & nn_2dint, llnightav ) 695 1086 696 END DO 1087 ENDIF 1088 1089 ! - Sea surface salinity 1090 IF ( ln_sss ) THEN 1091 IF(lwp) WRITE(numout,*) ' SSS currently not available' 1092 ENDIF 1093 1094 #if defined key_lim2 || defined key_lim3 1095 IF ( ln_seaice ) THEN 1096 DO jseaiceset = 1, nseaicesets 1097 CALL obs_seaice_opt( seaicedatqc(jseaiceset), & 1098 & kstp, jpi, jpj, nit000, 1.-frld, & 1099 & tmask(:,:,1), n2dint ) 1100 END DO 1101 ENDIF 697 698 ENDIF 699 700 CALL wrk_dealloc( jpi, jpj, jpk, zprofvar1 ) 701 CALL wrk_dealloc( jpi, jpj, jpk, zprofvar2 ) 702 CALL wrk_dealloc( jpi, jpj, jpk, zprofmask1 ) 703 CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 ) 704 CALL wrk_dealloc( jpi, jpj, zsurfvar ) 705 CALL wrk_dealloc( jpi, jpj, zglam1 ) 706 CALL wrk_dealloc( jpi, jpj, zglam2 ) 707 CALL wrk_dealloc( jpi, jpj, zgphi1 ) 708 CALL wrk_dealloc( jpi, jpj, zgphi2 ) 709 #if ! defined key_lim2 && ! defined key_lim3 710 CALL wrk_dealloc(jpi,jpj,frld) 1102 711 #endif 1103 712 1104 ! - Velocity profiles1105 IF ( ln_vel3d ) THEN1106 DO jveloset = 1, nvelosets1107 ! zonal component of velocity1108 CALL obs_vel_opt( veldatqc(jveloset), kstp, jpi, jpj, jpk, &1109 & nit000, idaystp, un, vn, gdept_1d, umask, vmask, &1110 n1dint, n2dint, ld_velav(jveloset) )1111 END DO1112 ENDIF1113 1114 #if ! defined key_lim2 && ! defined key_lim31115 CALL wrk_dealloc(jpi,jpj,frld)1116 #endif1117 1118 713 END SUBROUTINE dia_obs 1119 1120 SUBROUTINE dia_obs_wri 714 715 SUBROUTINE dia_obs_wri 1121 716 !!---------------------------------------------------------------------- 1122 717 !! *** ROUTINE dia_obs_wri *** … … 1126 721 !! ** Method : Call observation diagnostic output routines 1127 722 !! 1128 !! ** Action : 723 !! ** Action : 1129 724 !! 1130 725 !! History : … … 1134 729 !! ! 07-03 (K. Mogensen) General handling of profiles 1135 730 !! ! 08-09 (M. Valdivieso) Velocity component (U,V) profiles 1136 !!---------------------------------------------------------------------- 731 !! ! 15-08 (M. Martin) Combined writing for prof and surf types 732 !!---------------------------------------------------------------------- 733 !! * Modules used 734 USE obs_rot_vel ! Rotation of velocities 735 1137 736 IMPLICIT NONE 1138 737 1139 738 !! * Local declarations 1140 1141 INTEGER :: jprofset ! Profile data set loop variable 1142 INTEGER :: jveloset ! Velocity data set loop variable 1143 INTEGER :: jslaset ! SLA data set loop variable 1144 INTEGER :: jsstset ! SST data set loop variable 1145 INTEGER :: jseaiceset ! Sea Ice data set loop variable 1146 INTEGER :: jset 1147 INTEGER :: jfbini 1148 CHARACTER(LEN=20) :: datestr=" ",timestr=" " 1149 CHARACTER(LEN=10) :: cdtmp 739 INTEGER :: jtype ! Data set loop variable 740 INTEGER :: jo, jvar, jk 741 REAL(wp), DIMENSION(:), ALLOCATABLE :: & 742 & zu, & 743 & zv 744 1150 745 !----------------------------------------------------------------------- 1151 746 ! Depending on switches call various observation output routines 1152 747 !----------------------------------------------------------------------- 1153 748 1154 ! - Temperature/salinity profiles 1155 1156 IF( ln_t3d .OR. ln_s3d ) THEN 1157 1158 ! Copy data from prodatqc to profdata structures 1159 DO jprofset = 1, nprofsets 1160 1161 CALL obs_prof_decompress( prodatqc(jprofset), & 1162 & profdata(jprofset), .TRUE., numout ) 749 IF ( nproftypes > 0 ) THEN 750 751 DO jtype = 1, nproftypes 752 753 IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN 754 755 ! For velocity data, rotate the model velocities to N/S, E/W 756 ! using the compressed data structure. 757 ALLOCATE( & 758 & zu(profdataqc(jtype)%nvprot(1)), & 759 & zv(profdataqc(jtype)%nvprot(2)) & 760 & ) 761 762 CALL obs_rotvel( profdataqc(jtype), nn_2dint, zu, zv ) 763 764 DO jo = 1, profdataqc(jtype)%nprof 765 DO jvar = 1, 2 766 DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar) 767 768 IF ( jvar == 1 ) THEN 769 profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk) 770 ELSE 771 profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk) 772 ENDIF 773 774 END DO 775 END DO 776 END DO 777 778 DEALLOCATE( zu ) 779 DEALLOCATE( zv ) 780 781 END IF 782 783 CALL obs_prof_decompress( profdataqc(jtype), & 784 & profdata(jtype), .TRUE., numout ) 785 786 CALL obs_wri_prof( profdata(jtype) ) 1163 787 1164 788 END DO 1165 789 1166 ! Write the profiles. 1167 1168 jprofset = 0 1169 1170 ! ENACT insitu data 1171 1172 IF ( ln_ena ) THEN 1173 1174 jprofset = jprofset + 1 1175 1176 CALL obs_wri_p3d( 'enact', profdata(jprofset) ) 1177 1178 ENDIF 1179 1180 ! Coriolis insitu data 1181 1182 IF ( ln_cor ) THEN 1183 1184 jprofset = jprofset + 1 1185 1186 CALL obs_wri_p3d( 'corio', profdata(jprofset) ) 1187 1188 ENDIF 1189 1190 ! Feedback insitu data 1191 1192 IF ( ln_profb ) THEN 1193 1194 jfbini = jprofset + 1 1195 1196 DO jprofset = jfbini, nprofsets 1197 1198 jset = jprofset - jfbini + 1 1199 WRITE(cdtmp,'(A,I2.2)')'profb_',jset 1200 CALL obs_wri_p3d( cdtmp, profdata(jprofset) ) 1201 1202 END DO 1203 1204 ENDIF 1205 1206 ENDIF 1207 1208 ! - Sea surface anomaly 1209 IF ( ln_sla ) THEN 1210 1211 ! Copy data from sladatqc to sladata structures 1212 DO jslaset = 1, nslasets 1213 1214 CALL obs_surf_decompress( sladatqc(jslaset), & 1215 & sladata(jslaset), .TRUE., numout ) 790 ENDIF 791 792 IF ( nsurftypes > 0 ) THEN 793 794 DO jtype = 1, nsurftypes 795 796 CALL obs_surf_decompress( surfdataqc(jtype), & 797 & surfdata(jtype), .TRUE., numout ) 798 799 CALL obs_wri_surf( surfdata(jtype) ) 1216 800 1217 801 END DO 1218 802 1219 jslaset = 01220 1221 ! Write the AVISO SLA data1222 1223 IF ( ln_sladt ) THEN1224 1225 jslaset = 11226 CALL obs_wri_sla( 'aviso_act', sladata(jslaset) )1227 jslaset = 21228 CALL obs_wri_sla( 'aviso_pas', sladata(jslaset) )1229 1230 ENDIF1231 1232 IF ( ln_slafb ) THEN1233 1234 jfbini = jslaset + 11235 1236 DO jslaset = jfbini, nslasets1237 1238 jset = jslaset - jfbini + 11239 WRITE(cdtmp,'(A,I2.2)')'slafb_',jset1240 CALL obs_wri_sla( cdtmp, sladata(jslaset) )1241 1242 END DO1243 1244 ENDIF1245 1246 ENDIF1247 1248 ! - Sea surface temperature1249 IF ( ln_sst ) THEN1250 1251 ! Copy data from sstdatqc to sstdata structures1252 DO jsstset = 1, nsstsets1253 1254 CALL obs_surf_decompress( sstdatqc(jsstset), &1255 & sstdata(jsstset), .TRUE., numout )1256 1257 END DO1258 1259 jsstset = 01260 1261 ! Write the AVISO SST data1262 1263 IF ( ln_reysst ) THEN1264 1265 jsstset = jsstset + 11266 CALL obs_wri_sst( 'reynolds', sstdata(jsstset) )1267 1268 ENDIF1269 1270 IF ( ln_ghrsst ) THEN1271 1272 jsstset = jsstset + 11273 CALL obs_wri_sst( 'ghr', sstdata(jsstset) )1274 1275 ENDIF1276 1277 IF ( ln_sstfb ) THEN1278 1279 jfbini = jsstset + 11280 1281 DO jsstset = jfbini, nsstsets1282 1283 jset = jsstset - jfbini + 11284 WRITE(cdtmp,'(A,I2.2)')'sstfb_',jset1285 CALL obs_wri_sst( cdtmp, sstdata(jsstset) )1286 1287 END DO1288 1289 ENDIF1290 1291 ENDIF1292 1293 ! - Sea surface salinity1294 IF ( ln_sss ) THEN1295 IF(lwp) WRITE(numout,*) ' SSS currently not available'1296 ENDIF1297 1298 ! - Sea Ice Concentration1299 IF ( ln_seaice ) THEN1300 1301 ! Copy data from seaicedatqc to seaicedata structures1302 DO jseaiceset = 1, nseaicesets1303 1304 CALL obs_surf_decompress( seaicedatqc(jseaiceset), &1305 & seaicedata(jseaiceset), .TRUE., numout )1306 1307 END DO1308 1309 ! Write the Sea Ice data1310 DO jseaiceset = 1, nseaicesets1311 1312 WRITE(cdtmp,'(A,I2.2)')'seaicefb_',jseaiceset1313 CALL obs_wri_seaice( cdtmp, seaicedata(jseaiceset) )1314 1315 END DO1316 1317 ENDIF1318 1319 ! Velocity data1320 IF( ln_vel3d ) THEN1321 1322 ! Copy data from veldatqc to velodata structures1323 DO jveloset = 1, nvelosets1324 1325 CALL obs_prof_decompress( veldatqc(jveloset), &1326 & velodata(jveloset), .TRUE., numout )1327 1328 END DO1329 1330 ! Write the profiles.1331 1332 jveloset = 01333 1334 ! Daily averaged data1335 1336 IF ( ln_velavcur ) THEN1337 1338 jveloset = jveloset + 11339 1340 CALL obs_wri_vel( 'velavcurr', velodata(jveloset), n2dint )1341 1342 ENDIF1343 1344 ! High frequency data1345 1346 IF ( ln_velhrcur ) THEN1347 1348 jveloset = jveloset + 11349 1350 CALL obs_wri_vel( 'velhrcurr', velodata(jveloset), n2dint )1351 1352 ENDIF1353 1354 ! Daily averaged data1355 1356 IF ( ln_velavadcp ) THEN1357 1358 jveloset = jveloset + 11359 1360 CALL obs_wri_vel( 'velavadcp', velodata(jveloset), n2dint )1361 1362 ENDIF1363 1364 ! High frequency data1365 1366 IF ( ln_velhradcp ) THEN1367 1368 jveloset = jveloset + 11369 1370 CALL obs_wri_vel( 'velhradcp', velodata(jveloset), n2dint )1371 1372 ENDIF1373 1374 ! Feedback velocity data1375 1376 IF ( ln_velfb ) THEN1377 1378 jfbini = jveloset + 11379 1380 DO jveloset = jfbini, nvelosets1381 1382 jset = jveloset - jfbini + 11383 WRITE(cdtmp,'(A,I2.2)')'velfb_',jset1384 CALL obs_wri_vel( cdtmp, velodata(jveloset), n2dint )1385 1386 END DO1387 1388 ENDIF1389 1390 803 ENDIF 1391 804 … … 1405 818 !! 1406 819 !!---------------------------------------------------------------------- 1407 ! !obs_grid deallocation820 ! obs_grid deallocation 1408 821 CALL obs_grid_deallocate 1409 822 1410 !! diaobs deallocation 1411 IF ( nprofsets > 0 ) THEN 1412 DEALLOCATE(ld_enact, & 1413 & profdata, & 1414 & prodatqc) 1415 END IF 1416 IF ( ln_sla ) THEN 1417 DEALLOCATE(sladata, & 1418 & sladatqc) 1419 END IF 1420 IF ( ln_seaice ) THEN 1421 DEALLOCATE(sladata, & 1422 & sladatqc) 1423 END IF 1424 IF ( ln_sst ) THEN 1425 DEALLOCATE(sstdata, & 1426 & sstdatqc) 1427 END IF 1428 IF ( ln_vel3d ) THEN 1429 DEALLOCATE(ld_velav, & 1430 & velodata, & 1431 & veldatqc) 1432 END IF 823 ! diaobs deallocation 824 IF ( nproftypes > 0 ) & 825 & DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof ) 826 827 IF ( nsurftypes > 0 ) & 828 & DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf ) 829 1433 830 END SUBROUTINE dia_obs_dealloc 1434 831 1435 SUBROUTINE ini_date( ddobsini)1436 !!---------------------------------------------------------------------- 1437 !! *** ROUTINE ini_date ***832 SUBROUTINE calc_date( kstp, ddobs ) 833 !!---------------------------------------------------------------------- 834 !! *** ROUTINE calc_date *** 1438 835 !! 1439 !! ** Purpose : Get initial data in double precision YYYYMMDD.HHMMSS format 1440 !! 1441 !! ** Method : Get initial data in double precision YYYYMMDD.HHMMSS format 1442 !! 1443 !! ** Action : Get initial data in double precision YYYYMMDD.HHMMSS format 836 !! ** Purpose : Get date in double precision YYYYMMDD.HHMMSS format 837 !! 838 !! ** Method : Get date in double precision YYYYMMDD.HHMMSS format 839 !! 840 !! ** Action : Get date in double precision YYYYMMDD.HHMMSS format 841 !! 842 !! ** Action : Get initial date in double precision YYYYMMDD.HHMMSS format 1444 843 !! 1445 844 !! History : … … 1449 848 !! ! 06-10 (G. Smith) Calculates initial date the same as method for final date 1450 849 !! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2 850 !! ! 2014-09 (D. Lea) New generic routine now deals with arbitrary initial time of day 1451 851 !!---------------------------------------------------------------------- 1452 852 USE phycst, ONLY : & ! Physical constants 1453 853 & rday 1454 ! USE daymod, ONLY : & ! Time variables1455 ! & nmonth_len1456 854 USE dom_oce, ONLY : & ! Ocean space and time domain variables 1457 855 & rdt … … 1460 858 1461 859 !! * Arguments 1462 REAL(KIND=dp), INTENT(OUT) :: ddobsini ! Initial date in YYYYMMDD.HHMMSS 860 REAL(KIND=dp), INTENT(OUT) :: ddobs ! Date in YYYYMMDD.HHMMSS 861 INTEGER :: kstp 1463 862 1464 863 !! * Local declarations … … 1468 867 INTEGER :: ihou 1469 868 INTEGER :: imin 1470 INTEGER :: imday ! Number of days in month. 1471 REAL(KIND=wp) :: zdayfrc ! Fraction of day 869 INTEGER :: imday ! Number of days in month. 870 INTEGER, DIMENSION(12) :: & 871 & imonth_len ! Length in days of the months of the current year 872 REAL(wp) :: zdayfrc ! Fraction of day 1472 873 1473 874 INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year … … 1475 876 !!---------------------------------------------------------------------- 1476 877 !! Initial date initialization (year, month, day, hour, minute) 1477 !! (This assumes that the initial date is for 00z))1478 878 !!---------------------------------------------------------------------- 1479 879 iyea = ndate0 / 10000 1480 880 imon = ( ndate0 - iyea * 10000 ) / 100 1481 881 iday = ndate0 - iyea * 10000 - imon * 100 1482 ihou = 01483 imin = 0882 ihou = nn_time0 / 100 883 imin = ( nn_time0 - ihou * 100 ) 1484 884 1485 885 !!---------------------------------------------------------------------- 1486 886 !! Compute number of days + number of hours + min since initial time 1487 887 !!---------------------------------------------------------------------- 1488 iday = iday + ( nit000 -1 ) * rdt / rday 1489 zdayfrc = ( nit000 -1 ) * rdt / rday 888 zdayfrc = kstp * rdt / rday 1490 889 zdayfrc = zdayfrc - aint(zdayfrc) 1491 ihou = int( zdayfrc * 24 ) 1492 imin = int( (zdayfrc * 24 - ihou) * 60 ) 1493 1494 !!----------------------------------------------------------------------- 1495 !! Convert number of days (iday) into a real date 1496 !!---------------------------------------------------------------------- 890 imin = imin + int( zdayfrc * 24 * 60 ) 891 DO WHILE (imin >= 60) 892 imin=imin-60 893 ihou=ihou+1 894 END DO 895 DO WHILE (ihou >= 24) 896 ihou=ihou-24 897 iday=iday+1 898 END DO 899 iday = iday + kstp * rdt / rday 900 901 !----------------------------------------------------------------------- 902 ! Convert number of days (iday) into a real date 903 !---------------------------------------------------------------------- 1497 904 1498 905 CALL calc_month_len( iyea, imonth_len ) 1499 906 1500 907 DO WHILE ( iday > imonth_len(imon) ) 1501 908 iday = iday - imonth_len(imon) … … 1508 915 END DO 1509 916 1510 !!---------------------------------------------------------------------- 1511 !! Convert it into YYYYMMDD.HHMMSS format. 1512 !!---------------------------------------------------------------------- 1513 ddobsini = iyea * 10000_dp + imon * 100_dp + & 1514 & iday + ihou * 0.01_dp + imin * 0.0001_dp 1515 1516 1517 END SUBROUTINE ini_date 1518 1519 SUBROUTINE fin_date( ddobsfin ) 1520 !!---------------------------------------------------------------------- 1521 !! *** ROUTINE fin_date *** 917 !---------------------------------------------------------------------- 918 ! Convert it into YYYYMMDD.HHMMSS format. 919 !---------------------------------------------------------------------- 920 ddobs = iyea * 10000_dp + imon * 100_dp + & 921 & iday + ihou * 0.01_dp + imin * 0.0001_dp 922 923 END SUBROUTINE calc_date 924 925 SUBROUTINE ini_date( ddobsini ) 926 !!---------------------------------------------------------------------- 927 !! *** ROUTINE ini_date *** 1522 928 !! 1523 !! ** Purpose : Get final datain double precision YYYYMMDD.HHMMSS format1524 !! 1525 !! ** Method : Get final data in double precision YYYYMMDD.HHMMSS format1526 !! 1527 !! ** Action : Get final data in double precision YYYYMMDD.HHMMSS format929 !! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format 930 !! 931 !! ** Method : 932 !! 933 !! ** Action : 1528 934 !! 1529 935 !! History : … … 1532 938 !! ! 06-10 (A. Weaver) Cleaning 1533 939 !! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2 1534 !!---------------------------------------------------------------------- 1535 USE phycst, ONLY : & ! Physical constants 1536 & rday 1537 ! USE daymod, ONLY : & ! Time variables 1538 ! & nmonth_len 1539 USE dom_oce, ONLY : & ! Ocean space and time domain variables 1540 & rdt 940 !! ! 2014-09 (D. Lea) Change to call generic routine calc_date 941 !!---------------------------------------------------------------------- 1541 942 1542 943 IMPLICIT NONE 1543 944 1544 945 !! * Arguments 1545 REAL(KIND=dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS 1546 1547 !! * Local declarations 1548 INTEGER :: iyea ! date - (year, month, day, hour, minute) 1549 INTEGER :: imon 1550 INTEGER :: iday 1551 INTEGER :: ihou 1552 INTEGER :: imin 1553 INTEGER :: imday ! Number of days in month. 1554 REAL(KIND=wp) :: zdayfrc ! Fraction of day 1555 1556 INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year 1557 1558 !----------------------------------------------------------------------- 1559 ! Initial date initialization (year, month, day, hour, minute) 1560 ! (This assumes that the initial date is for 00z) 1561 !----------------------------------------------------------------------- 1562 iyea = ndate0 / 10000 1563 imon = ( ndate0 - iyea * 10000 ) / 100 1564 iday = ndate0 - iyea * 10000 - imon * 100 1565 ihou = 0 1566 imin = 0 1567 1568 !----------------------------------------------------------------------- 1569 ! Compute number of days + number of hours + min since initial time 1570 !----------------------------------------------------------------------- 1571 iday = iday + nitend * rdt / rday 1572 zdayfrc = nitend * rdt / rday 1573 zdayfrc = zdayfrc - AINT( zdayfrc ) 1574 ihou = INT( zdayfrc * 24 ) 1575 imin = INT( ( zdayfrc * 24 - ihou ) * 60 ) 1576 1577 !----------------------------------------------------------------------- 1578 ! Convert number of days (iday) into a real date 1579 !---------------------------------------------------------------------- 1580 1581 CALL calc_month_len( iyea, imonth_len ) 1582 1583 DO WHILE ( iday > imonth_len(imon) ) 1584 iday = iday - imonth_len(imon) 1585 imon = imon + 1 1586 IF ( imon > 12 ) THEN 1587 imon = 1 1588 iyea = iyea + 1 1589 CALL calc_month_len( iyea, imonth_len ) ! update month lengths 1590 ENDIF 1591 END DO 1592 1593 !----------------------------------------------------------------------- 1594 ! Convert it into YYYYMMDD.HHMMSS format 1595 !----------------------------------------------------------------------- 1596 ddobsfin = iyea * 10000_dp + imon * 100_dp + iday & 1597 & + ihou * 0.01_dp + imin * 0.0001_dp 1598 1599 END SUBROUTINE fin_date 1600 946 REAL(KIND=dp), INTENT(OUT) :: ddobsini ! Initial date in YYYYMMDD.HHMMSS 947 948 CALL calc_date( nit000 - 1, ddobsini ) 949 950 END SUBROUTINE ini_date 951 952 SUBROUTINE fin_date( ddobsfin ) 953 !!---------------------------------------------------------------------- 954 !! *** ROUTINE fin_date *** 955 !! 956 !! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format 957 !! 958 !! ** Method : 959 !! 960 !! ** Action : 961 !! 962 !! History : 963 !! ! 06-03 (K. Mogensen) Original code 964 !! ! 06-05 (K. Mogensen) Reformatted 965 !! ! 06-10 (A. Weaver) Cleaning 966 !! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2 967 !! ! 2014-09 (D. Lea) Change to call generic routine calc_date 968 !!---------------------------------------------------------------------- 969 970 IMPLICIT NONE 971 972 !! * Arguments 973 REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS 974 975 CALL calc_date( nitend, ddobsfin ) 976 977 END SUBROUTINE fin_date 978 1601 979 END MODULE diaobs
Note: See TracChangeset
for help on using the changeset viewer.