- Timestamp:
- 2016-01-08T10:35:19+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
r4624 r6225 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 !!---------------------------------------------------------------------- … … 135 96 !! ! 06-10 (A. Weaver) Cleaning and add controls 136 97 !! ! 07-03 (K. Mogensen) General handling of profiles 98 !! ! 14-08 (J.While) Incorporated SST bias correction 99 !! ! 15-02 (M. Martin) Simplification of namelist and code 137 100 !!---------------------------------------------------------------------- 138 101 … … 140 103 141 104 !! * 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, & 105 INTEGER, PARAMETER :: & 106 & jpmaxnfiles = 1000 ! Maximum number of files for each obs type 107 INTEGER, DIMENSION(:), ALLOCATABLE :: & 108 & ifilesprof, & ! Number of profile files 109 & ifilessurf ! Number of surface files 110 INTEGER :: ios ! Local integer output status for namelist read 111 INTEGER :: jtype ! Counter for obs types 112 INTEGER :: jvar ! Counter for variables 113 INTEGER :: jfile ! Counter for files 114 115 CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: & 116 & cn_profbfiles, & ! T/S profile input filenames 117 & cn_sstfbfiles, & ! Sea surface temperature input filenames 118 & cn_slafbfiles, & ! Sea level anomaly input filenames 119 & cn_sicfbfiles, & ! Seaice concentration input filenames 120 & cn_velfbfiles, & ! Velocity profile input filenames 121 & cn_sstbias_files ! SST bias input filenames 122 CHARACTER(LEN=128) :: & 123 & cn_altbiasfile ! Altimeter bias input filename 124 CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: & 125 & clproffiles, & ! Profile filenames 126 & clsurffiles ! Surface filenames 127 128 LOGICAL :: ln_t3d ! Logical switch for temperature profiles 129 LOGICAL :: ln_s3d ! Logical switch for salinity profiles 130 LOGICAL :: ln_sla ! Logical switch for sea level anomalies 131 LOGICAL :: ln_sst ! Logical switch for sea surface temperature 132 LOGICAL :: ln_sic ! Logical switch for sea ice concentration 133 LOGICAL :: ln_vel3d ! Logical switch for velocity (u,v) obs 134 LOGICAL :: ln_nea ! Logical switch to remove obs near land 135 LOGICAL :: ln_altbias ! Logical switch for altimeter bias 136 LOGICAL :: ln_sstbias !: Logical switch for bias corection of SST 137 LOGICAL :: ln_ignmis ! Logical switch for ignoring missing files 138 LOGICAL :: ln_s_at_t ! Logical switch to compute model S at T obs 139 LOGICAL :: llvar1 ! Logical for profile variable 1 140 LOGICAL :: llvar2 ! Logical for profile variable 1 141 LOGICAL :: llnightav ! Logical for calculating night-time averages 142 LOGICAL, DIMENSION(jpmaxnfiles) :: lmask ! Used for finding number of sstbias files 143 144 REAL(dp) :: rn_dobsini ! Obs window start date YYYYMMDD.HHMMSS 145 REAL(dp) :: rn_dobsend ! Obs window end date YYYYMMDD.HHMMSS 146 REAL(wp), POINTER, DIMENSION(:,:) :: & 147 & zglam1, & ! Model longitudes for profile variable 1 148 & zglam2 ! Model longitudes for profile variable 2 149 REAL(wp), POINTER, DIMENSION(:,:) :: & 150 & zgphi1, & ! Model latitudes for profile variable 1 151 & zgphi2 ! Model latitudes for profile variable 2 152 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 153 & zmask1, & ! Model land/sea mask associated with variable 1 154 & zmask2 ! Model land/sea mask associated with variable 2 155 156 NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla, & 157 & ln_sst, ln_sic, ln_vel3d, & 158 & ln_altbias, ln_nea, ln_grid_global, & 173 159 & 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 160 & ln_ignmis, ln_s_at_t, ln_sstnight, & 161 & cn_profbfiles, cn_slafbfiles, & 162 & cn_sstfbfiles, cn_sicfbfiles, & 163 & cn_velfbfiles, cn_altbiasfile, & 164 & cn_gridsearchfile, rn_gridsearchres, & 165 & rn_dobsini, rn_dobsend, nn_1dint, nn_2dint, & 166 & nn_msshc, rn_mdtcorr, rn_mdtcutoff, & 167 & nn_profdavtypes, ln_sstbias, cn_sstbias_files 168 169 INTEGER :: jnumsstbias 170 CALL wrk_alloc( jpi, jpj, zglam1 ) 171 CALL wrk_alloc( jpi, jpj, zglam2 ) 172 CALL wrk_alloc( jpi, jpj, zgphi1 ) 173 CALL wrk_alloc( jpi, jpj, zgphi2 ) 174 CALL wrk_alloc( jpi, jpj, jpk, zmask1 ) 175 CALL wrk_alloc( jpi, jpj, jpk, zmask2 ) 205 176 206 177 !----------------------------------------------------------------------- 207 178 ! Read namelist parameters 208 179 !----------------------------------------------------------------------- 209 210 velcurfiles(:) = '' 211 veladcpfiles(:) = '' 212 CALL ini_date( dobsini ) 213 CALL fin_date( dobsend ) 214 215 ! Read Namelist namobs : control observation diagnostics 216 REWIND( numnam_ref ) ! Namelist namobs in reference namelist : Diagnostic: control observation 180 181 !Initalise all values in namelist arrays 182 ALLOCATE(sstbias_type(jpmaxnfiles)) 183 ! Some namelist arrays need initialising 184 cn_profbfiles(:) = '' 185 cn_slafbfiles(:) = '' 186 cn_sstfbfiles(:) = '' 187 cn_sicfbfiles(:) = '' 188 cn_velfbfiles(:) = '' 189 cn_sstbias_files(:) = '' 190 nn_profdavtypes(:) = -1 191 192 CALL ini_date( rn_dobsini ) 193 CALL fin_date( rn_dobsend ) 194 195 ! Read namelist namobs : control observation diagnostics 196 REWIND( numnam_ref ) ! Namelist namobs in reference namelist 217 197 READ ( numnam_ref, namobs, IOSTAT = ios, ERR = 901) 218 198 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp ) 219 199 220 REWIND( numnam_cfg ) ! Namelist namobs in configuration namelist : Diagnostic: control observation200 REWIND( numnam_cfg ) ! Namelist namobs in configuration namelist 221 201 READ ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 ) 222 202 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp ) 223 203 IF(lwm) WRITE ( numond, namobs ) 224 204 225 ! Count number of files for each type 226 IF (ln_ena) THEN 227 lmask(:) = .FALSE. 228 WHERE (enactfiles(:) /= '') lmask(:) = .TRUE. 229 jnumenact = COUNT(lmask) 230 ENDIF 231 IF (ln_cor) THEN 232 lmask(:) = .FALSE. 233 WHERE (coriofiles(:) /= '') lmask(:) = .TRUE. 234 jnumcorio = COUNT(lmask) 235 ENDIF 236 IF (ln_profb) THEN 237 lmask(:) = .FALSE. 238 WHERE (profbfiles(:) /= '') lmask(:) = .TRUE. 239 jnumprofb = COUNT(lmask) 240 ENDIF 241 IF (ln_sladt) THEN 242 lmask(:) = .FALSE. 243 WHERE (slafilesact(:) /= '') lmask(:) = .TRUE. 244 jnumslaact = COUNT(lmask) 245 lmask(:) = .FALSE. 246 WHERE (slafilespas(:) /= '') lmask(:) = .TRUE. 247 jnumslapas = COUNT(lmask) 248 ENDIF 249 IF (ln_slafb) THEN 250 lmask(:) = .FALSE. 251 WHERE (slafbfiles(:) /= '') lmask(:) = .TRUE. 252 jnumslafb = COUNT(lmask) 253 lmask(:) = .FALSE. 254 ENDIF 255 IF (ln_ghrsst) THEN 256 lmask(:) = .FALSE. 257 WHERE (sstfiles(:) /= '') lmask(:) = .TRUE. 258 jnumsst = COUNT(lmask) 205 IF ( .NOT. ln_diaobs ) THEN 206 IF(lwp) WRITE(numout,cform_war) 207 IF(lwp) WRITE(numout,*)' ln_diaobs is set to false so not calling dia_obs' 208 RETURN 209 ENDIF 210 211 !----------------------------------------------------------------------- 212 ! Set up list of observation types to be used 213 ! and the files associated with each type 214 !----------------------------------------------------------------------- 215 216 nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 217 nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic /) ) 218 219 IF (ln_sstbias) THEN 220 lmask(:) = .FALSE. 221 WHERE (cn_sstbias_files(:) /= '') lmask(:) = .TRUE. 222 jnumsstbias = COUNT(lmask) 223 lmask(:) = .FALSE. 259 224 ENDIF 260 IF (ln_sstfb) THEN 261 lmask(:) = .FALSE. 262 WHERE (sstfbfiles(:) /= '') lmask(:) = .TRUE. 263 jnumsstfb = COUNT(lmask) 264 lmask(:) = .FALSE. 265 ENDIF 266 IF (ln_seaice) THEN 267 lmask(:) = .FALSE. 268 WHERE (seaicefiles(:) /= '') lmask(:) = .TRUE. 269 jnumseaice = COUNT(lmask) 270 ENDIF 271 IF (ln_velavcur) THEN 272 lmask(:) = .FALSE. 273 WHERE (velavcurfiles(:) /= '') lmask(:) = .TRUE. 274 jnumvelavcur = COUNT(lmask) 275 ENDIF 276 IF (ln_velhrcur) THEN 277 lmask(:) = .FALSE. 278 WHERE (velhrcurfiles(:) /= '') lmask(:) = .TRUE. 279 jnumvelhrcur = COUNT(lmask) 280 ENDIF 281 IF (ln_velavadcp) THEN 282 lmask(:) = .FALSE. 283 WHERE (velavadcpfiles(:) /= '') lmask(:) = .TRUE. 284 jnumvelavadcp = COUNT(lmask) 285 ENDIF 286 IF (ln_velhradcp) THEN 287 lmask(:) = .FALSE. 288 WHERE (velhradcpfiles(:) /= '') lmask(:) = .TRUE. 289 jnumvelhradcp = COUNT(lmask) 290 ENDIF 291 IF (ln_velfb) THEN 292 lmask(:) = .FALSE. 293 WHERE (velfbfiles(:) /= '') lmask(:) = .TRUE. 294 jnumvelfb = COUNT(lmask) 295 lmask(:) = .FALSE. 296 ENDIF 297 298 ! Control print 225 226 IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 227 IF(lwp) WRITE(numout,cform_war) 228 IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 229 & ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', & 230 & ' are set to .FALSE. so turning off calls to dia_obs' 231 nwarn = nwarn + 1 232 ln_diaobs = .FALSE. 233 RETURN 234 ENDIF 235 236 IF ( nproftypes > 0 ) THEN 237 238 ALLOCATE( cobstypesprof(nproftypes) ) 239 ALLOCATE( ifilesprof(nproftypes) ) 240 ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 241 242 jtype = 0 243 IF (ln_t3d .OR. ln_s3d) THEN 244 jtype = jtype + 1 245 clproffiles(jtype,:) = cn_profbfiles(:) 246 cobstypesprof(jtype) = 'prof ' 247 ifilesprof(jtype) = 0 248 DO jfile = 1, jpmaxnfiles 249 IF ( trim(clproffiles(jtype,jfile)) /= '' ) & 250 ifilesprof(jtype) = ifilesprof(jtype) + 1 251 END DO 252 ENDIF 253 IF (ln_vel3d) THEN 254 jtype = jtype + 1 255 clproffiles(jtype,:) = cn_velfbfiles(:) 256 cobstypesprof(jtype) = 'vel ' 257 ifilesprof(jtype) = 0 258 DO jfile = 1, jpmaxnfiles 259 IF ( trim(clproffiles(jtype,jfile)) /= '' ) & 260 ifilesprof(jtype) = ifilesprof(jtype) + 1 261 END DO 262 ENDIF 263 264 ENDIF 265 266 IF ( nsurftypes > 0 ) THEN 267 268 ALLOCATE( cobstypessurf(nsurftypes) ) 269 ALLOCATE( ifilessurf(nsurftypes) ) 270 ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 271 272 jtype = 0 273 IF (ln_sla) THEN 274 jtype = jtype + 1 275 clsurffiles(jtype,:) = cn_slafbfiles(:) 276 cobstypessurf(jtype) = 'sla ' 277 ifilessurf(jtype) = 0 278 DO jfile = 1, jpmaxnfiles 279 IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 280 ifilessurf(jtype) = ifilessurf(jtype) + 1 281 END DO 282 ENDIF 283 IF (ln_sst) THEN 284 jtype = jtype + 1 285 clsurffiles(jtype,:) = cn_sstfbfiles(:) 286 cobstypessurf(jtype) = 'sst ' 287 ifilessurf(jtype) = 0 288 DO jfile = 1, jpmaxnfiles 289 IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 290 ifilessurf(jtype) = ifilessurf(jtype) + 1 291 END DO 292 ENDIF 293 #if defined key_lim2 || defined key_lim3 294 IF (ln_sic) THEN 295 jtype = jtype + 1 296 clsurffiles(jtype,:) = cn_sicfbfiles(:) 297 cobstypessurf(jtype) = 'sic ' 298 ifilessurf(jtype) = 0 299 DO jfile = 1, jpmaxnfiles 300 IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 301 ifilessurf(jtype) = ifilessurf(jtype) + 1 302 END DO 303 ENDIF 304 #endif 305 306 ENDIF 307 308 !Write namelist settings to stdout 299 309 IF(lwp) THEN 300 310 WRITE(numout,*) … … 302 312 WRITE(numout,*) '~~~~~~~~~~~~' 303 313 WRITE(numout,*) ' Namelist namobs : set observation diagnostic parameters' 304 WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d 305 WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d 306 WRITE(numout,*) ' Logical switch for ENACT insitu data set ln_ena = ', ln_ena 307 WRITE(numout,*) ' Logical switch for Coriolis insitu data set ln_cor = ', ln_cor 308 WRITE(numout,*) ' Logical switch for feedback insitu data set ln_profb = ', ln_profb 309 WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla 310 WRITE(numout,*) ' Logical switch for AVISO SLA data ln_sladt = ', ln_sladt 311 WRITE(numout,*) ' Logical switch for feedback SLA data ln_slafb = ', ln_slafb 312 WRITE(numout,*) ' Logical switch for SSH observations ln_ssh = ', ln_ssh 313 WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst 314 WRITE(numout,*) ' Logical switch for Reynolds observations ln_reysst = ', ln_reysst 315 WRITE(numout,*) ' Logical switch for GHRSST observations ln_ghrsst = ', ln_ghrsst 316 WRITE(numout,*) ' Logical switch for feedback SST data ln_sstfb = ', ln_sstfb 317 WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight 318 WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss 319 WRITE(numout,*) ' Logical switch for Sea Ice observations ln_seaice = ', ln_seaice 320 WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d 321 WRITE(numout,*) ' Logical switch for velocity daily av. cur. ln_velavcur = ', ln_velavcur 322 WRITE(numout,*) ' Logical switch for velocity high freq. cur. ln_velhrcur = ', ln_velhrcur 323 WRITE(numout,*) ' Logical switch for velocity daily av. ADCP ln_velavadcp = ', ln_velavadcp 324 WRITE(numout,*) ' Logical switch for velocity high freq. ADCP ln_velhradcp = ', ln_velhradcp 325 WRITE(numout,*) ' Logical switch for feedback velocity data ln_velfb = ', ln_velfb 326 WRITE(numout,*) ' Global distribtion of observations ln_grid_global = ',ln_grid_global 327 WRITE(numout,*) & 328 ' Logical switch for obs grid search w/lookup table ln_grid_search_lookup = ',ln_grid_search_lookup 314 WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d 315 WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d 316 WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla 317 WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst 318 WRITE(numout,*) ' Logical switch for Sea Ice observations ln_sic = ', ln_sic 319 WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d 320 WRITE(numout,*) ' Global distribution of observations ln_grid_global = ',ln_grid_global 321 WRITE(numout,*) ' Logical switch for SST bias correction ln_sstbias = ', ln_sstbias 322 WRITE(numout,*) ' Logical switch for obs grid search lookup ln_grid_search_lookup = ',ln_grid_search_lookup 329 323 IF (ln_grid_search_lookup) & 330 WRITE(numout,*) ' Grid search lookup file header grid_search_file = ', grid_search_file 331 IF (ln_ena) THEN 332 DO ji = 1, jnumenact 333 WRITE(numout,'(1X,2A)') ' ENACT input observation file name enactfiles = ', & 334 TRIM(enactfiles(ji)) 324 WRITE(numout,*) ' Grid search lookup file header cn_gridsearchfile = ', cn_gridsearchfile 325 WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS rn_dobsini = ', rn_dobsini 326 WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS rn_dobsend = ', rn_dobsend 327 WRITE(numout,*) ' Type of vertical interpolation method nn_1dint = ', nn_1dint 328 WRITE(numout,*) ' Type of horizontal interpolation method nn_2dint = ', nn_2dint 329 WRITE(numout,*) ' Rejection of observations near land switch ln_nea = ', ln_nea 330 WRITE(numout,*) ' MSSH correction scheme nn_msshc = ', nn_msshc 331 WRITE(numout,*) ' MDT correction rn_mdtcorr = ', rn_mdtcorr 332 WRITE(numout,*) ' MDT cutoff for computed correction rn_mdtcutoff = ', rn_mdtcutoff 333 WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias 334 WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis 335 WRITE(numout,*) ' Daily average types nn_profdavtypes = ', nn_profdavtypes 336 WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight 337 WRITE(numout,*) ' Number of profile obs types: ',nproftypes 338 339 IF ( nproftypes > 0 ) THEN 340 DO jtype = 1, nproftypes 341 DO jfile = 1, ifilesprof(jtype) 342 WRITE(numout,'(1X,2A)') ' '//cobstypesprof(jtype)//' input observation file names = ', & 343 TRIM(clproffiles(jtype,jfile)) 344 END DO 335 345 END DO 336 346 ENDIF 337 IF (ln_cor) THEN 338 DO ji = 1, jnumcorio 339 WRITE(numout,'(1X,2A)') ' Coriolis input observation file name coriofiles = ', & 340 TRIM(coriofiles(ji)) 347 348 WRITE(numout,*)' Number of surface obs types: ',nsurftypes 349 IF ( nsurftypes > 0 ) THEN 350 DO jtype = 1, nsurftypes 351 DO jfile = 1, ifilessurf(jtype) 352 WRITE(numout,'(1X,2A)') ' '//cobstypessurf(jtype)//' input observation file names = ', & 353 TRIM(clsurffiles(jtype,jfile)) 354 END DO 341 355 END DO 342 356 ENDIF 343 IF (ln_profb) THEN 344 DO ji = 1, jnumprofb 345 IF (ln_profb_ena(ji)) THEN 346 WRITE(numout,'(1X,2A)') ' Enact feedback input observation file name profbfiles = ', & 347 TRIM(profbfiles(ji)) 348 ELSE 349 WRITE(numout,'(1X,2A)') ' Feedback input observation file name profbfiles = ', & 350 TRIM(profbfiles(ji)) 351 ENDIF 352 WRITE(numout,'(1X,2A)') ' Enact feedback input time setting switch ln_profb_enatim = ', ln_profb_enatim(ji) 353 END DO 354 ENDIF 355 IF (ln_sladt) THEN 356 DO ji = 1, jnumslaact 357 WRITE(numout,'(1X,2A)') ' Active SLA input observation file name slafilesact = ', & 358 TRIM(slafilesact(ji)) 359 END DO 360 DO ji = 1, jnumslapas 361 WRITE(numout,'(1X,2A)') ' Passive SLA input observation file name slafilespas = ', & 362 TRIM(slafilespas(ji)) 363 END DO 364 ENDIF 365 IF (ln_slafb) THEN 366 DO ji = 1, jnumslafb 367 WRITE(numout,'(1X,2A)') ' Feedback SLA input observation file name slafbfiles = ', & 368 TRIM(slafbfiles(ji)) 369 END DO 370 ENDIF 371 IF (ln_ghrsst) THEN 372 DO ji = 1, jnumsst 373 WRITE(numout,'(1X,2A)') ' GHRSST input observation file name sstfiles = ', & 374 TRIM(sstfiles(ji)) 375 END DO 376 ENDIF 377 IF (ln_sstfb) THEN 378 DO ji = 1, jnumsstfb 379 WRITE(numout,'(1X,2A)') ' Feedback SST input observation file name sstfbfiles = ', & 380 TRIM(sstfbfiles(ji)) 381 END DO 382 ENDIF 383 IF (ln_seaice) THEN 384 DO ji = 1, jnumseaice 385 WRITE(numout,'(1X,2A)') ' Sea Ice input observation file name seaicefiles = ', & 386 TRIM(seaicefiles(ji)) 387 END DO 388 ENDIF 389 IF (ln_velavcur) THEN 390 DO ji = 1, jnumvelavcur 391 WRITE(numout,'(1X,2A)') ' Vel. cur. daily av. input file name velavcurfiles = ', & 392 TRIM(velavcurfiles(ji)) 393 END DO 394 ENDIF 395 IF (ln_velhrcur) THEN 396 DO ji = 1, jnumvelhrcur 397 WRITE(numout,'(1X,2A)') ' Vel. cur. high freq. input file name velhvcurfiles = ', & 398 TRIM(velhrcurfiles(ji)) 399 END DO 400 ENDIF 401 IF (ln_velavadcp) THEN 402 DO ji = 1, jnumvelavadcp 403 WRITE(numout,'(1X,2A)') ' Vel. ADCP daily av. input file name velavadcpfiles = ', & 404 TRIM(velavadcpfiles(ji)) 405 END DO 406 ENDIF 407 IF (ln_velhradcp) THEN 408 DO ji = 1, jnumvelhradcp 409 WRITE(numout,'(1X,2A)') ' Vel. ADCP high freq. input file name velhvadcpfiles = ', & 410 TRIM(velhradcpfiles(ji)) 411 END DO 412 ENDIF 413 IF (ln_velfb) THEN 414 DO ji = 1, jnumvelfb 415 IF (ln_velfb_av(ji)) THEN 416 WRITE(numout,'(1X,2A)') ' Vel. feedback daily av. input file name velfbfiles = ', & 417 TRIM(velfbfiles(ji)) 418 ELSE 419 WRITE(numout,'(1X,2A)') ' Vel. feedback input observation file name velfbfiles = ', & 420 TRIM(velfbfiles(ji)) 421 ENDIF 422 END DO 423 ENDIF 424 WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS dobsini = ', dobsini 425 WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS dobsend = ', dobsend 426 WRITE(numout,*) ' Type of vertical interpolation method n1dint = ', n1dint 427 WRITE(numout,*) ' Type of horizontal interpolation method n2dint = ', n2dint 428 WRITE(numout,*) ' Rejection of observations near land swithch ln_nea = ', ln_nea 429 WRITE(numout,*) ' MSSH correction scheme nmsshc = ', nmsshc 430 WRITE(numout,*) ' MDT correction mdtcorr = ', mdtcorr 431 WRITE(numout,*) ' MDT cutoff for computed correction mdtcutoff = ', mdtcutoff 432 WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias 433 WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis 434 WRITE(numout,*) ' ENACT daily average types = ',endailyavtypes 435 436 ENDIF 437 357 WRITE(numout,*) '~~~~~~~~~~~~' 358 359 ENDIF 360 361 !----------------------------------------------------------------------- 362 ! Obs operator parameter checking and initialisations 363 !----------------------------------------------------------------------- 364 438 365 IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN 439 366 CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' ) … … 441 368 ENDIF 442 369 370 IF ( ln_grid_global ) THEN 371 CALL ctl_warn( 'ln_grid_global=T may cause memory issues when used with a large number of processors' ) 372 ENDIF 373 374 IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN 375 CALL ctl_stop(' Choice of vertical (1D) interpolation method', & 376 & ' is not available') 377 ENDIF 378 379 IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 4 ) ) THEN 380 CALL ctl_stop(' Choice of horizontal (2D) interpolation method', & 381 & ' is not available') 382 ENDIF 383 443 384 CALL obs_typ_init 444 445 CALL mppmap_init 446 447 ! Parameter control 448 #if defined key_diaobs 449 IF ( ( .NOT. ln_t3d ).AND.( .NOT. ln_s3d ).AND.( .NOT. ln_sla ).AND. & 450 & ( .NOT. ln_vel3d ).AND. & 451 & ( .NOT. ln_ssh ).AND.( .NOT. ln_sst ).AND.( .NOT. ln_sss ).AND. & 452 & ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ) ) THEN 453 IF(lwp) WRITE(numout,cform_war) 454 IF(lwp) WRITE(numout,*) ' key_diaobs is activated but logical flags', & 455 & ' ln_t3d, ln_s3d, ln_sla, ln_ssh, ln_sst, ln_sss, ln_seaice, ln_vel3d are all set to .FALSE.' 456 nwarn = nwarn + 1 457 ENDIF 458 #endif 385 IF(ln_grid_global) THEN 386 CALL mppmap_init 387 ENDIF 459 388 460 389 CALL obs_grid_setup( ) 461 IF ( ( n1dint < 0 ).OR.( n1dint > 1 ) ) THEN462 IF(lwp) WRITE(numout,cform_err)463 IF(lwp) WRITE(numout,*) ' Choice of vertical (1D) interpolation method', &464 & ' is not available'465 nstop = nstop + 1466 ENDIF467 IF ( ( n2dint < 0 ).OR.( n2dint > 4 ) ) THEN468 IF(lwp) WRITE(numout,cform_err)469 IF(lwp) WRITE(numout,*) ' Choice of horizontal (2D) interpolation method', &470 & ' is not available'471 nstop = nstop + 1472 ENDIF473 390 474 391 !----------------------------------------------------------------------- 475 392 ! Depending on switches read the various observation types 476 393 !----------------------------------------------------------------------- 477 ! - Temperature/salinity profiles 478 479 IF ( ln_t3d .OR. ln_s3d ) THEN 480 481 ! Set the number of variables for profiles to 2 (T and S) 482 nprofvars = 2 483 ! Set the number of extra variables for profiles to 1 (insitu temp). 484 nprofextr = 1 485 486 ! Count how may insitu data sets we have and allocate data. 487 jprofset = 0 488 IF ( ln_ena ) jprofset = jprofset + 1 489 IF ( ln_cor ) jprofset = jprofset + 1 490 IF ( ln_profb ) jprofset = jprofset + jnumprofb 491 nprofsets = jprofset 492 IF ( nprofsets > 0 ) THEN 493 ALLOCATE(ld_enact(nprofsets)) 494 ALLOCATE(profdata(nprofsets)) 495 ALLOCATE(prodatqc(nprofsets)) 496 ENDIF 497 498 jprofset = 0 499 500 ! ENACT insitu data 501 502 IF ( ln_ena ) THEN 503 504 jprofset = jprofset + 1 505 506 ld_enact(jprofset) = .TRUE. 507 508 CALL obs_rea_pro_dri( 1, profdata(jprofset), & 509 & jnumenact, enactfiles(1:jnumenact), & 510 & nprofvars, nprofextr, & 511 & nitend-nit000+2, & 512 & dobsini, dobsend, ln_t3d, ln_s3d, & 513 & ln_ignmis, ln_s_at_t, .TRUE., .FALSE., & 514 & kdailyavtypes = endailyavtypes ) 515 516 DO jvar = 1, 2 517 518 CALL obs_prof_staend( profdata(jprofset), jvar ) 519 394 395 IF ( nproftypes > 0 ) THEN 396 397 ALLOCATE(profdata(nproftypes)) 398 ALLOCATE(profdataqc(nproftypes)) 399 ALLOCATE(nvarsprof(nproftypes)) 400 ALLOCATE(nextrprof(nproftypes)) 401 402 DO jtype = 1, nproftypes 403 404 nvarsprof(jtype) = 2 405 IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN 406 nextrprof(jtype) = 1 407 llvar1 = ln_t3d 408 llvar2 = ln_s3d 409 zglam1 = glamt 410 zgphi1 = gphit 411 zmask1 = tmask 412 zglam2 = glamt 413 zgphi2 = gphit 414 zmask2 = tmask 415 ENDIF 416 IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN 417 nextrprof(jtype) = 2 418 llvar1 = ln_vel3d 419 llvar2 = ln_vel3d 420 zglam1 = glamu 421 zgphi1 = gphiu 422 zmask1 = umask 423 zglam2 = glamv 424 zgphi2 = gphiv 425 zmask2 = vmask 426 ENDIF 427 428 !Read in profile or profile obs types 429 CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype), & 430 & clproffiles(jtype,1:ifilesprof(jtype)), & 431 & nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, & 432 & rn_dobsini, rn_dobsend, llvar1, llvar2, & 433 & ln_ignmis, ln_s_at_t, .FALSE., & 434 & kdailyavtypes = nn_profdavtypes ) 435 436 DO jvar = 1, nvarsprof(jtype) 437 CALL obs_prof_staend( profdata(jtype), jvar ) 520 438 END DO 521 439 522 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 523 & ln_t3d, ln_s3d, ln_nea, & 524 & kdailyavtypes=endailyavtypes ) 525 526 ENDIF 527 528 ! Coriolis insitu data 529 530 IF ( ln_cor ) THEN 531 532 jprofset = jprofset + 1 533 534 ld_enact(jprofset) = .FALSE. 535 536 CALL obs_rea_pro_dri( 2, profdata(jprofset), & 537 & jnumcorio, coriofiles(1:jnumcorio), & 538 & nprofvars, nprofextr, & 539 & nitend-nit000+2, & 540 & dobsini, dobsend, ln_t3d, ln_s3d, & 541 & ln_ignmis, ln_s_at_t, .FALSE., .FALSE. ) 542 543 DO jvar = 1, 2 544 545 CALL obs_prof_staend( profdata(jprofset), jvar ) 546 547 END DO 548 549 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 550 & ln_t3d, ln_s3d, ln_nea ) 551 552 ENDIF 553 554 ! Feedback insitu data 555 556 IF ( ln_profb ) THEN 557 558 DO jset = 1, jnumprofb 559 560 jprofset = jprofset + 1 561 ld_enact (jprofset) = ln_profb_ena(jset) 562 563 CALL obs_rea_pro_dri( 0, profdata(jprofset), & 564 & 1, profbfiles(jset:jset), & 565 & nprofvars, nprofextr, & 566 & nitend-nit000+2, & 567 & dobsini, dobsend, ln_t3d, ln_s3d, & 568 & ln_ignmis, ln_s_at_t, & 569 & ld_enact(jprofset).AND.& 570 & ln_profb_enatim(jset), & 571 & .FALSE., kdailyavtypes = endailyavtypes ) 572 573 DO jvar = 1, 2 574 575 CALL obs_prof_staend( profdata(jprofset), jvar ) 576 577 END DO 578 579 IF ( ld_enact(jprofset) ) THEN 580 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 581 & ln_t3d, ln_s3d, ln_nea, & 582 & kdailyavtypes = endailyavtypes ) 583 ELSE 584 CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), & 585 & ln_t3d, ln_s3d, ln_nea ) 586 ENDIF 587 588 END DO 589 590 ENDIF 591 592 ENDIF 593 594 ! - Sea level anomalies 595 IF ( ln_sla ) THEN 596 ! Set the number of variables for sla to 1 597 nslavars = 1 598 599 ! Set the number of extra variables for sla to 2 600 nslaextr = 2 440 CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), & 441 & llvar1, llvar2, & 442 & jpi, jpj, jpk, & 443 & zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2, & 444 & ln_nea, kdailyavtypes = nn_profdavtypes ) 445 446 END DO 447 448 DEALLOCATE( ifilesprof, clproffiles ) 449 450 ENDIF 451 452 IF ( nsurftypes > 0 ) THEN 453 454 ALLOCATE(surfdata(nsurftypes)) 455 ALLOCATE(surfdataqc(nsurftypes)) 456 ALLOCATE(nvarssurf(nsurftypes)) 457 ALLOCATE(nextrsurf(nsurftypes)) 458 459 DO jtype = 1, nsurftypes 460 461 nvarssurf(jtype) = 1 462 nextrsurf(jtype) = 0 463 llnightav = .FALSE. 464 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2 465 IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav = ln_sstnight 466 467 !Read in surface obs types 468 CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), & 469 & clsurffiles(jtype,1:ifilessurf(jtype)), & 470 & nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 471 & rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav ) 601 472 602 ! Set the number of sla data sets to 2603 nslasets = 0604 IF ( ln_sladt ) THEN605 nslasets = nslasets + 2606 ENDIF607 IF ( ln_slafb ) THEN608 nslasets = nslasets + jnumslafb609 ENDIF610 473 611 ALLOCATE(sladata(nslasets)) 612 ALLOCATE(sladatqc(nslasets)) 613 sladata(:)%nsurf=0 614 sladatqc(:)%nsurf=0 615 616 nslasets = 0 617 618 ! AVISO SLA data 619 620 IF ( ln_sladt ) THEN 621 622 ! Active SLA observations 623 624 nslasets = nslasets + 1 625 626 CALL obs_rea_sla( 1, sladata(nslasets), jnumslaact, & 627 & slafilesact(1:jnumslaact), & 628 & nslavars, nslaextr, nitend-nit000+2, & 629 & dobsini, dobsend, ln_ignmis, .FALSE. ) 630 CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), & 631 & ln_sla, ln_nea ) 632 633 ! Passive SLA observations 634 635 nslasets = nslasets + 1 636 637 CALL obs_rea_sla( 1, sladata(nslasets), jnumslapas, & 638 & slafilespas(1:jnumslapas), & 639 & nslavars, nslaextr, nitend-nit000+2, & 640 & dobsini, dobsend, ln_ignmis, .FALSE. ) 641 642 CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), & 643 & ln_sla, ln_nea ) 644 645 ENDIF 646 647 ! Feedback SLA data 648 649 IF ( ln_slafb ) THEN 650 651 DO jset = 1, jnumslafb 652 653 nslasets = nslasets + 1 654 655 CALL obs_rea_sla( 0, sladata(nslasets), 1, & 656 & slafbfiles(jset:jset), & 657 & nslavars, nslaextr, nitend-nit000+2, & 658 & dobsini, dobsend, ln_ignmis, .FALSE. ) 659 CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), & 660 & ln_sla, ln_nea ) 661 662 END DO 663 664 ENDIF 665 666 CALL obs_rea_mdt( nslasets, sladatqc, n2dint ) 667 668 ! read in altimeter bias 669 670 IF ( ln_altbias ) THEN 671 CALL obs_rea_altbias ( nslasets, sladatqc, n2dint, bias_file ) 672 ENDIF 673 674 ENDIF 675 676 ! - Sea surface height 677 IF ( ln_ssh ) THEN 678 IF(lwp) WRITE(numout,*) ' SSH currently not available' 679 ENDIF 680 681 ! - Sea surface temperature 682 IF ( ln_sst ) THEN 683 684 ! Set the number of variables for sst to 1 685 nsstvars = 1 686 687 ! Set the number of extra variables for sst to 0 688 nsstextr = 0 689 690 nsstsets = 0 691 692 IF (ln_reysst) nsstsets = nsstsets + 1 693 IF (ln_ghrsst) nsstsets = nsstsets + 1 694 IF ( ln_sstfb ) THEN 695 nsstsets = nsstsets + jnumsstfb 696 ENDIF 697 698 ALLOCATE(sstdata(nsstsets)) 699 ALLOCATE(sstdatqc(nsstsets)) 700 ALLOCATE(ld_sstnight(nsstsets)) 701 sstdata(:)%nsurf=0 702 sstdatqc(:)%nsurf=0 703 ld_sstnight(:)=.false. 704 705 nsstsets = 0 706 707 IF (ln_reysst) THEN 708 709 nsstsets = nsstsets + 1 710 711 ld_sstnight(nsstsets) = ln_sstnight 712 713 CALL obs_rea_sst_rey( reysstname, reysstfmt, sstdata(nsstsets), & 714 & nsstvars, nsstextr, & 715 & nitend-nit000+2, dobsini, dobsend ) 716 CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, & 717 & ln_nea ) 718 719 ENDIF 720 721 IF (ln_ghrsst) THEN 722 723 nsstsets = nsstsets + 1 724 725 ld_sstnight(nsstsets) = ln_sstnight 726 727 CALL obs_rea_sst( 1, sstdata(nsstsets), jnumsst, & 728 & sstfiles(1:jnumsst), & 729 & nsstvars, nsstextr, nitend-nit000+2, & 730 & dobsini, dobsend, ln_ignmis, .FALSE. ) 731 CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, & 732 & ln_nea ) 733 734 ENDIF 735 736 ! Feedback SST data 737 738 IF ( ln_sstfb ) THEN 739 740 DO jset = 1, jnumsstfb 741 742 nsstsets = nsstsets + 1 743 744 ld_sstnight(nsstsets) = ln_sstnight 745 746 CALL obs_rea_sst( 0, sstdata(nsstsets), 1, & 747 & sstfbfiles(jset:jset), & 748 & nsstvars, nsstextr, nitend-nit000+2, & 749 & dobsini, dobsend, ln_ignmis, .FALSE. ) 750 CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), & 751 & ln_sst, ln_nea ) 752 753 END DO 754 755 ENDIF 756 757 ENDIF 758 759 ! - Sea surface salinity 760 IF ( ln_sss ) THEN 761 IF(lwp) WRITE(numout,*) ' SSS currently not available' 762 ENDIF 763 764 ! - Sea Ice Concentration 765 766 IF ( ln_seaice ) THEN 767 768 ! Set the number of variables for seaice to 1 769 nseaicevars = 1 770 771 ! Set the number of extra variables for seaice to 0 772 nseaiceextr = 0 773 774 ! Set the number of data sets to 1 775 nseaicesets = 1 776 777 ALLOCATE(seaicedata(nseaicesets)) 778 ALLOCATE(seaicedatqc(nseaicesets)) 779 seaicedata(:)%nsurf=0 780 seaicedatqc(:)%nsurf=0 781 782 CALL obs_rea_seaice( 1, seaicedata(nseaicesets), jnumseaice, & 783 & seaicefiles(1:jnumseaice), & 784 & nseaicevars, nseaiceextr, nitend-nit000+2, & 785 & dobsini, dobsend, ln_ignmis, .FALSE. ) 786 787 CALL obs_pre_seaice( seaicedata(nseaicesets), seaicedatqc(nseaicesets), & 788 & ln_seaice, ln_nea ) 789 790 ENDIF 791 792 IF (ln_vel3d) THEN 793 794 ! Set the number of variables for profiles to 2 (U and V) 795 nvelovars = 2 796 797 ! Set the number of extra variables for profiles to 2 to store 798 ! rotation parameters 799 nveloextr = 2 800 801 jveloset = 0 802 803 IF ( ln_velavcur ) jveloset = jveloset + 1 804 IF ( ln_velhrcur ) jveloset = jveloset + 1 805 IF ( ln_velavadcp ) jveloset = jveloset + 1 806 IF ( ln_velhradcp ) jveloset = jveloset + 1 807 IF (ln_velfb) jveloset = jveloset + jnumvelfb 808 809 nvelosets = jveloset 810 IF ( nvelosets > 0 ) THEN 811 ALLOCATE( velodata(nvelosets) ) 812 ALLOCATE( veldatqc(nvelosets) ) 813 ALLOCATE( ld_velav(nvelosets) ) 814 ENDIF 815 816 jveloset = 0 817 818 ! Daily averaged data 819 820 IF ( ln_velavcur ) THEN 821 822 jveloset = jveloset + 1 823 824 ld_velav(jveloset) = .TRUE. 825 826 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavcur, & 827 & velavcurfiles(1:jnumvelavcur), & 828 & nvelovars, nveloextr, & 829 & nitend-nit000+2, & 830 & dobsini, dobsend, ln_ignmis, & 831 & ld_velav(jveloset), & 832 & .FALSE. ) 833 834 DO jvar = 1, 2 835 CALL obs_prof_staend( velodata(jveloset), jvar ) 836 END DO 837 838 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 839 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 840 841 ENDIF 842 843 ! High frequency data 844 845 IF ( ln_velhrcur ) THEN 846 847 jveloset = jveloset + 1 848 849 ld_velav(jveloset) = .FALSE. 850 851 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhrcur, & 852 & velhrcurfiles(1:jnumvelhrcur), & 853 & nvelovars, nveloextr, & 854 & nitend-nit000+2, & 855 & dobsini, dobsend, ln_ignmis, & 856 & ld_velav(jveloset), & 857 & .FALSE. ) 858 859 DO jvar = 1, 2 860 CALL obs_prof_staend( velodata(jveloset), jvar ) 861 END DO 862 863 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 864 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 865 866 ENDIF 867 868 ! Daily averaged data 869 870 IF ( ln_velavadcp ) THEN 871 872 jveloset = jveloset + 1 873 874 ld_velav(jveloset) = .TRUE. 875 876 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavadcp, & 877 & velavadcpfiles(1:jnumvelavadcp), & 878 & nvelovars, nveloextr, & 879 & nitend-nit000+2, & 880 & dobsini, dobsend, ln_ignmis, & 881 & ld_velav(jveloset), & 882 & .FALSE. ) 883 884 DO jvar = 1, 2 885 CALL obs_prof_staend( velodata(jveloset), jvar ) 886 END DO 887 888 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 889 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 890 891 ENDIF 892 893 ! High frequency data 894 895 IF ( ln_velhradcp ) THEN 896 897 jveloset = jveloset + 1 898 899 ld_velav(jveloset) = .FALSE. 900 901 CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhradcp, & 902 & velhradcpfiles(1:jnumvelhradcp), & 903 & nvelovars, nveloextr, & 904 & nitend-nit000+2, & 905 & dobsini, dobsend, ln_ignmis, & 906 & ld_velav(jveloset), & 907 & .FALSE. ) 908 909 DO jvar = 1, 2 910 CALL obs_prof_staend( velodata(jveloset), jvar ) 911 END DO 912 913 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 914 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 915 916 ENDIF 917 918 IF ( ln_velfb ) THEN 919 920 DO jset = 1, jnumvelfb 921 922 jveloset = jveloset + 1 923 924 ld_velav(jveloset) = ln_velfb_av(jset) 925 926 CALL obs_rea_vel_dri( 0, velodata(jveloset), 1, & 927 & velfbfiles(jset:jset), & 928 & nvelovars, nveloextr, & 929 & nitend-nit000+2, & 930 & dobsini, dobsend, ln_ignmis, & 931 & ld_velav(jveloset), & 932 & .FALSE. ) 933 934 DO jvar = 1, 2 935 CALL obs_prof_staend( velodata(jveloset), jvar ) 936 END DO 937 938 CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 939 & ln_vel3d, ln_nea, ld_velav(jveloset) ) 940 941 942 END DO 943 944 ENDIF 945 946 ENDIF 947 474 CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea ) 475 476 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 477 CALL obs_rea_mdt( surfdataqc(jtype), nn_2dint ) 478 IF ( ln_altbias ) CALL obs_rea_altbias ( surfdataqc(jtype), nn_2dint, cn_altbiasfile ) 479 ENDIF 480 IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN 481 !Read in bias field and correct SST. 482 IF ( jnumsstbias == 0 ) CALL ctl_stop("ln_sstbias set,"// & 483 " but no bias"// & 484 " files to read in") 485 CALL obs_app_sstbias( surfdataqc(jtype), nn_2dint, & 486 jnumsstbias, cn_sstbias_files(1:jnumsstbias) ) 487 ENDIF 488 END DO 489 490 DEALLOCATE( ifilessurf, clsurffiles ) 491 492 ENDIF 493 494 CALL wrk_dealloc( jpi, jpj, zglam1 ) 495 CALL wrk_dealloc( jpi, jpj, zglam2 ) 496 CALL wrk_dealloc( jpi, jpj, zgphi1 ) 497 CALL wrk_dealloc( jpi, jpj, zgphi2 ) 498 CALL wrk_dealloc( jpi, jpj, jpk, zmask1 ) 499 CALL wrk_dealloc( jpi, jpj, jpk, zmask2 ) 500 948 501 END SUBROUTINE dia_obs_init 949 502 … … 955 508 !! 956 509 !! ** Method : Call the observation operators on each time step to 957 !! compute the model equivalent of the following date: 958 !! - T profiles 959 !! - S profiles 960 !! - Sea surface height (referenced to a mean) 961 !! - Sea surface temperature 962 !! - Sea surface salinity 963 !! - Velocity component (U,V) profiles 964 !! 965 !! ** Action : 510 !! compute the model equivalent of the following data: 511 !! - Profile data, currently T/S or U/V 512 !! - Surface data, currently SST, SLA or sea-ice concentration. 513 !! 514 !! ** Action : 966 515 !! 967 516 !! History : … … 972 521 !! ! 07-04 (G. Smith) Generalized surface operators 973 522 !! ! 08-10 (M. Valdivieso) obs operator for velocity profiles 523 !! ! 14-08 (J. While) observation operator for profiles in 524 !! generalised vertical coordinates 525 !! ! 15-08 (M. Martin) Combined surface/profile routines. 974 526 !!---------------------------------------------------------------------- 975 527 !! * Modules used 976 528 USE dom_oce, ONLY : & ! Ocean space and time domain variables 977 & rdt, & 978 & gdept_1d, & 979 & tmask, umask, vmask 529 & gdept_n, & 530 & gdept_1d 980 531 USE phycst, ONLY : & ! Physical constants 981 532 & rday 982 533 USE oce, ONLY : & ! Ocean dynamics and tracers variables 983 534 & tsn, & 984 & un, vn, & 985 & sshn 535 & un, vn, & 536 & sshn 537 USE phycst, ONLY : & ! Physical constants 538 & rday 986 539 #if defined key_lim3 987 USE ice, ONLY : & ! LIMIce model variables540 USE ice, ONLY : & ! LIM3 Ice model variables 988 541 & frld 989 542 #endif 990 543 #if defined key_lim2 991 USE ice_2, ONLY : & ! LIMIce model variables544 USE ice_2, ONLY : & ! LIM2 Ice model variables 992 545 & frld 993 546 #endif … … 995 548 996 549 !! * Arguments 997 INTEGER, INTENT(IN) :: kstp 550 INTEGER, INTENT(IN) :: kstp ! Current timestep 998 551 !! * Local declarations 999 INTEGER :: idaystp ! Number of timesteps per day 1000 INTEGER :: jprofset ! Profile data set loop variable 1001 INTEGER :: jslaset ! SLA data set loop variable 1002 INTEGER :: jsstset ! SST data set loop variable 1003 INTEGER :: jseaiceset ! sea ice data set loop variable 1004 INTEGER :: jveloset ! velocity profile data loop variable 1005 INTEGER :: jvar ! Variable number 552 INTEGER :: idaystp ! Number of timesteps per day 553 INTEGER :: jtype ! Data loop variable 554 INTEGER :: jvar ! Variable number 555 INTEGER :: ji, jj ! Loop counters 556 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 557 & zprofvar1, & ! Model values for 1st variable in a prof ob 558 & zprofvar2 ! Model values for 2nd variable in a prof ob 559 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 560 & zprofmask1, & ! Mask associated with zprofvar1 561 & zprofmask2 ! Mask associated with zprofvar2 562 REAL(wp), POINTER, DIMENSION(:,:) :: & 563 & zsurfvar ! Model values equivalent to surface ob. 564 REAL(wp), POINTER, DIMENSION(:,:) :: & 565 & zglam1, & ! Model longitudes for prof variable 1 566 & zglam2, & ! Model longitudes for prof variable 2 567 & zgphi1, & ! Model latitudes for prof variable 1 568 & zgphi2 ! Model latitudes for prof variable 2 1006 569 #if ! defined key_lim2 && ! defined key_lim3 1007 REAL(wp), POINTER, DIMENSION(:,:) :: frld 570 REAL(wp), POINTER, DIMENSION(:,:) :: frld 1008 571 #endif 1009 CHARACTER(LEN=20) :: datestr=" ",timestr=" " 1010 572 LOGICAL :: llnightav ! Logical for calculating night-time average 573 574 !Allocate local work arrays 575 CALL wrk_alloc( jpi, jpj, jpk, zprofvar1 ) 576 CALL wrk_alloc( jpi, jpj, jpk, zprofvar2 ) 577 CALL wrk_alloc( jpi, jpj, jpk, zprofmask1 ) 578 CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 ) 579 CALL wrk_alloc( jpi, jpj, zsurfvar ) 580 CALL wrk_alloc( jpi, jpj, zglam1 ) 581 CALL wrk_alloc( jpi, jpj, zglam2 ) 582 CALL wrk_alloc( jpi, jpj, zgphi1 ) 583 CALL wrk_alloc( jpi, jpj, zgphi2 ) 1011 584 #if ! defined key_lim2 && ! defined key_lim3 1012 585 CALL wrk_alloc(jpi,jpj,frld) … … 1028 601 #endif 1029 602 !----------------------------------------------------------------------- 1030 ! Depending on switches call various observation operators 1031 !----------------------------------------------------------------------- 1032 1033 ! - Temperature/salinity profiles 1034 IF ( ln_t3d .OR. ln_s3d ) THEN 1035 DO jprofset = 1, nprofsets 1036 IF ( ld_enact(jprofset) ) THEN 1037 CALL obs_pro_opt( prodatqc(jprofset), & 1038 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1039 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1040 & gdept_1d, tmask, n1dint, n2dint, & 1041 & kdailyavtypes = endailyavtypes ) 603 ! Call the profile and surface observation operators 604 !----------------------------------------------------------------------- 605 606 IF ( nproftypes > 0 ) THEN 607 608 DO jtype = 1, nproftypes 609 610 SELECT CASE ( TRIM(cobstypesprof(jtype)) ) 611 CASE('prof') 612 zprofvar1(:,:,:) = tsn(:,:,:,jp_tem) 613 zprofvar2(:,:,:) = tsn(:,:,:,jp_sal) 614 zprofmask1(:,:,:) = tmask(:,:,:) 615 zprofmask2(:,:,:) = tmask(:,:,:) 616 zglam1(:,:) = glamt(:,:) 617 zglam2(:,:) = glamt(:,:) 618 zgphi1(:,:) = gphit(:,:) 619 zgphi2(:,:) = gphit(:,:) 620 CASE('vel') 621 zprofvar1(:,:,:) = un(:,:,:) 622 zprofvar2(:,:,:) = vn(:,:,:) 623 zprofmask1(:,:,:) = umask(:,:,:) 624 zprofmask2(:,:,:) = vmask(:,:,:) 625 zglam1(:,:) = glamu(:,:) 626 zglam2(:,:) = glamv(:,:) 627 zgphi1(:,:) = gphiu(:,:) 628 zgphi2(:,:) = gphiv(:,:) 629 END SELECT 630 631 IF( ln_zco .OR. ln_zps ) THEN 632 CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk, & 633 & nit000, idaystp, & 634 & zprofvar1, zprofvar2, & 635 & gdept_1d, zprofmask1, zprofmask2, & 636 & zglam1, zglam2, zgphi1, zgphi2, & 637 & nn_1dint, nn_2dint, & 638 & kdailyavtypes = nn_profdavtypes ) 639 ELSE IF(TRIM(cobstypesprof(jtype)) == 'prof') THEN 640 !TG - THIS NEEDS MODIFICATION TO MATCH SIMPLIFICATION 641 CALL obs_pro_sco_opt( profdataqc(jtype), & 642 & kstp, jpi, jpj, jpk, nit000, idaystp, & 643 & zprofvar1, zprofvar2, & 644 & gdept_n(:,:,:), gdepw_n(:,:,:), & 645 & tmask, nn_1dint, nn_2dint, & 646 & kdailyavtypes = nn_profdavtypes ) 1042 647 ELSE 1043 CALL obs_pro_opt( prodatqc(jprofset), & 1044 & kstp, jpi, jpj, jpk, nit000, idaystp, & 1045 & tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), & 1046 & gdept_1d, tmask, n1dint, n2dint ) 648 CALL ctl_stop('DIA_OBS: Generalised vertical interpolation not'// & 649 'yet working for velocity data (turn off velocity observations') 1047 650 ENDIF 651 1048 652 END DO 1049 ENDIF 1050 1051 ! - Sea surface anomaly 1052 IF ( ln_sla ) THEN 1053 DO jslaset = 1, nslasets 1054 CALL obs_sla_opt( sladatqc(jslaset), & 1055 & kstp, jpi, jpj, nit000, sshn, & 1056 & tmask(:,:,1), n2dint ) 1057 END DO 1058 ENDIF 1059 1060 ! - Sea surface temperature 1061 IF ( ln_sst ) THEN 1062 DO jsstset = 1, nsstsets 1063 CALL obs_sst_opt( sstdatqc(jsstset), & 1064 & kstp, jpi, jpj, nit000, idaystp, & 1065 & tsn(:,:,1,jp_tem), tmask(:,:,1), & 1066 & n2dint, ld_sstnight(jsstset) ) 653 654 ENDIF 655 656 IF ( nsurftypes > 0 ) THEN 657 658 DO jtype = 1, nsurftypes 659 660 SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 661 CASE('sst') 662 zsurfvar(:,:) = tsn(:,:,1,jp_tem) 663 llnightav = ln_sstnight 664 CASE('sla') 665 zsurfvar(:,:) = sshn(:,:) 666 llnightav = .FALSE. 667 #if defined key_lim2 || defined key_lim3 668 CASE('sic') 669 IF ( kstp == 0 ) THEN 670 IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN 671 CALL ctl_warn( 'Sea-ice not initialised on zeroth '// & 672 & 'time-step but some obs are valid then.' ) 673 WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & 674 & ' sea-ice obs will be missed' 675 ENDIF 676 surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + & 677 & surfdataqc(jtype)%nsstp(1) 678 CYCLE 679 ELSE 680 zsurfvar(:,:) = 1._wp - frld(:,:) 681 ENDIF 682 683 llnightav = .FALSE. 684 #endif 685 END SELECT 686 687 CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj, & 688 & nit000, idaystp, zsurfvar, tmask(:,:,1), & 689 & nn_2dint, llnightav ) 690 1067 691 END DO 1068 ENDIF 1069 1070 ! - Sea surface salinity 1071 IF ( ln_sss ) THEN 1072 IF(lwp) WRITE(numout,*) ' SSS currently not available' 1073 ENDIF 1074 1075 #if defined key_lim2 || defined key_lim3 1076 IF ( ln_seaice ) THEN 1077 DO jseaiceset = 1, nseaicesets 1078 CALL obs_seaice_opt( seaicedatqc(jseaiceset), & 1079 & kstp, jpi, jpj, nit000, 1.-frld, & 1080 & tmask(:,:,1), n2dint ) 1081 END DO 1082 ENDIF 692 693 ENDIF 694 695 CALL wrk_dealloc( jpi, jpj, jpk, zprofvar1 ) 696 CALL wrk_dealloc( jpi, jpj, jpk, zprofvar2 ) 697 CALL wrk_dealloc( jpi, jpj, jpk, zprofmask1 ) 698 CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 ) 699 CALL wrk_dealloc( jpi, jpj, zsurfvar ) 700 CALL wrk_dealloc( jpi, jpj, zglam1 ) 701 CALL wrk_dealloc( jpi, jpj, zglam2 ) 702 CALL wrk_dealloc( jpi, jpj, zgphi1 ) 703 CALL wrk_dealloc( jpi, jpj, zgphi2 ) 704 #if ! defined key_lim2 && ! defined key_lim3 705 CALL wrk_dealloc(jpi,jpj,frld) 1083 706 #endif 1084 707 1085 ! - Velocity profiles1086 IF ( ln_vel3d ) THEN1087 DO jveloset = 1, nvelosets1088 ! zonal component of velocity1089 CALL obs_vel_opt( veldatqc(jveloset), kstp, jpi, jpj, jpk, &1090 & nit000, idaystp, un, vn, gdept_1d, umask, vmask, &1091 n1dint, n2dint, ld_velav(jveloset) )1092 END DO1093 ENDIF1094 1095 #if ! defined key_lim2 && ! defined key_lim31096 CALL wrk_dealloc(jpi,jpj,frld)1097 #endif1098 1099 708 END SUBROUTINE dia_obs 1100 1101 SUBROUTINE dia_obs_wri 709 710 SUBROUTINE dia_obs_wri 1102 711 !!---------------------------------------------------------------------- 1103 712 !! *** ROUTINE dia_obs_wri *** … … 1107 716 !! ** Method : Call observation diagnostic output routines 1108 717 !! 1109 !! ** Action : 718 !! ** Action : 1110 719 !! 1111 720 !! History : … … 1115 724 !! ! 07-03 (K. Mogensen) General handling of profiles 1116 725 !! ! 08-09 (M. Valdivieso) Velocity component (U,V) profiles 1117 !!---------------------------------------------------------------------- 726 !! ! 15-08 (M. Martin) Combined writing for prof and surf types 727 !!---------------------------------------------------------------------- 728 !! * Modules used 729 USE obs_rot_vel ! Rotation of velocities 730 1118 731 IMPLICIT NONE 1119 732 1120 733 !! * Local declarations 1121 1122 INTEGER :: jprofset ! Profile data set loop variable 1123 INTEGER :: jveloset ! Velocity data set loop variable 1124 INTEGER :: jslaset ! SLA data set loop variable 1125 INTEGER :: jsstset ! SST data set loop variable 1126 INTEGER :: jseaiceset ! Sea Ice data set loop variable 1127 INTEGER :: jset 1128 INTEGER :: jfbini 1129 CHARACTER(LEN=20) :: datestr=" ",timestr=" " 1130 CHARACTER(LEN=10) :: cdtmp 734 INTEGER :: jtype ! Data set loop variable 735 INTEGER :: jo, jvar, jk 736 REAL(wp), DIMENSION(:), ALLOCATABLE :: & 737 & zu, & 738 & zv 739 1131 740 !----------------------------------------------------------------------- 1132 741 ! Depending on switches call various observation output routines 1133 742 !----------------------------------------------------------------------- 1134 743 1135 ! - Temperature/salinity profiles 1136 1137 IF( ln_t3d .OR. ln_s3d ) THEN 1138 1139 ! Copy data from prodatqc to profdata structures 1140 DO jprofset = 1, nprofsets 1141 1142 CALL obs_prof_decompress( prodatqc(jprofset), & 1143 & profdata(jprofset), .TRUE., numout ) 744 IF ( nproftypes > 0 ) THEN 745 746 DO jtype = 1, nproftypes 747 748 IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN 749 750 ! For velocity data, rotate the model velocities to N/S, E/W 751 ! using the compressed data structure. 752 ALLOCATE( & 753 & zu(profdataqc(jtype)%nvprot(1)), & 754 & zv(profdataqc(jtype)%nvprot(2)) & 755 & ) 756 757 CALL obs_rotvel( profdataqc(jtype), nn_2dint, zu, zv ) 758 759 DO jo = 1, profdataqc(jtype)%nprof 760 DO jvar = 1, 2 761 DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar) 762 763 IF ( jvar == 1 ) THEN 764 profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk) 765 ELSE 766 profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk) 767 ENDIF 768 769 END DO 770 END DO 771 END DO 772 773 DEALLOCATE( zu ) 774 DEALLOCATE( zv ) 775 776 END IF 777 778 CALL obs_prof_decompress( profdataqc(jtype), & 779 & profdata(jtype), .TRUE., numout ) 780 781 CALL obs_wri_prof( profdata(jtype) ) 1144 782 1145 783 END DO 1146 784 1147 ! Write the profiles. 1148 1149 jprofset = 0 1150 1151 ! ENACT insitu data 1152 1153 IF ( ln_ena ) THEN 1154 1155 jprofset = jprofset + 1 1156 1157 CALL obs_wri_p3d( 'enact', profdata(jprofset) ) 1158 1159 ENDIF 1160 1161 ! Coriolis insitu data 1162 1163 IF ( ln_cor ) THEN 1164 1165 jprofset = jprofset + 1 1166 1167 CALL obs_wri_p3d( 'corio', profdata(jprofset) ) 1168 1169 ENDIF 1170 1171 ! Feedback insitu data 1172 1173 IF ( ln_profb ) THEN 1174 1175 jfbini = jprofset + 1 1176 1177 DO jprofset = jfbini, nprofsets 1178 1179 jset = jprofset - jfbini + 1 1180 WRITE(cdtmp,'(A,I2.2)')'profb_',jset 1181 CALL obs_wri_p3d( cdtmp, profdata(jprofset) ) 1182 1183 END DO 1184 1185 ENDIF 1186 1187 ENDIF 1188 1189 ! - Sea surface anomaly 1190 IF ( ln_sla ) THEN 1191 1192 ! Copy data from sladatqc to sladata structures 1193 DO jslaset = 1, nslasets 1194 1195 CALL obs_surf_decompress( sladatqc(jslaset), & 1196 & sladata(jslaset), .TRUE., numout ) 785 ENDIF 786 787 IF ( nsurftypes > 0 ) THEN 788 789 DO jtype = 1, nsurftypes 790 791 CALL obs_surf_decompress( surfdataqc(jtype), & 792 & surfdata(jtype), .TRUE., numout ) 793 794 CALL obs_wri_surf( surfdata(jtype) ) 1197 795 1198 796 END DO 1199 797 1200 jslaset = 01201 1202 ! Write the AVISO SLA data1203 1204 IF ( ln_sladt ) THEN1205 1206 jslaset = 11207 CALL obs_wri_sla( 'aviso_act', sladata(jslaset) )1208 jslaset = 21209 CALL obs_wri_sla( 'aviso_pas', sladata(jslaset) )1210 1211 ENDIF1212 1213 IF ( ln_slafb ) THEN1214 1215 jfbini = jslaset + 11216 1217 DO jslaset = jfbini, nslasets1218 1219 jset = jslaset - jfbini + 11220 WRITE(cdtmp,'(A,I2.2)')'slafb_',jset1221 CALL obs_wri_sla( cdtmp, sladata(jslaset) )1222 1223 END DO1224 1225 ENDIF1226 1227 ENDIF1228 1229 ! - Sea surface temperature1230 IF ( ln_sst ) THEN1231 1232 ! Copy data from sstdatqc to sstdata structures1233 DO jsstset = 1, nsstsets1234 1235 CALL obs_surf_decompress( sstdatqc(jsstset), &1236 & sstdata(jsstset), .TRUE., numout )1237 1238 END DO1239 1240 jsstset = 01241 1242 ! Write the AVISO SST data1243 1244 IF ( ln_reysst ) THEN1245 1246 jsstset = jsstset + 11247 CALL obs_wri_sst( 'reynolds', sstdata(jsstset) )1248 1249 ENDIF1250 1251 IF ( ln_ghrsst ) THEN1252 1253 jsstset = jsstset + 11254 CALL obs_wri_sst( 'ghr', sstdata(jsstset) )1255 1256 ENDIF1257 1258 IF ( ln_sstfb ) THEN1259 1260 jfbini = jsstset + 11261 1262 DO jsstset = jfbini, nsstsets1263 1264 jset = jsstset - jfbini + 11265 WRITE(cdtmp,'(A,I2.2)')'sstfb_',jset1266 CALL obs_wri_sst( cdtmp, sstdata(jsstset) )1267 1268 END DO1269 1270 ENDIF1271 1272 ENDIF1273 1274 ! - Sea surface salinity1275 IF ( ln_sss ) THEN1276 IF(lwp) WRITE(numout,*) ' SSS currently not available'1277 ENDIF1278 1279 ! - Sea Ice Concentration1280 IF ( ln_seaice ) THEN1281 1282 ! Copy data from seaicedatqc to seaicedata structures1283 DO jseaiceset = 1, nseaicesets1284 1285 CALL obs_surf_decompress( seaicedatqc(jseaiceset), &1286 & seaicedata(jseaiceset), .TRUE., numout )1287 1288 END DO1289 1290 ! Write the Sea Ice data1291 DO jseaiceset = 1, nseaicesets1292 1293 WRITE(cdtmp,'(A,I2.2)')'seaicefb_',jseaiceset1294 CALL obs_wri_seaice( cdtmp, seaicedata(jseaiceset) )1295 1296 END DO1297 1298 ENDIF1299 1300 ! Velocity data1301 IF( ln_vel3d ) THEN1302 1303 ! Copy data from veldatqc to velodata structures1304 DO jveloset = 1, nvelosets1305 1306 CALL obs_prof_decompress( veldatqc(jveloset), &1307 & velodata(jveloset), .TRUE., numout )1308 1309 END DO1310 1311 ! Write the profiles.1312 1313 jveloset = 01314 1315 ! Daily averaged data1316 1317 IF ( ln_velavcur ) THEN1318 1319 jveloset = jveloset + 11320 1321 CALL obs_wri_vel( 'velavcurr', velodata(jveloset), n2dint )1322 1323 ENDIF1324 1325 ! High frequency data1326 1327 IF ( ln_velhrcur ) THEN1328 1329 jveloset = jveloset + 11330 1331 CALL obs_wri_vel( 'velhrcurr', velodata(jveloset), n2dint )1332 1333 ENDIF1334 1335 ! Daily averaged data1336 1337 IF ( ln_velavadcp ) THEN1338 1339 jveloset = jveloset + 11340 1341 CALL obs_wri_vel( 'velavadcp', velodata(jveloset), n2dint )1342 1343 ENDIF1344 1345 ! High frequency data1346 1347 IF ( ln_velhradcp ) THEN1348 1349 jveloset = jveloset + 11350 1351 CALL obs_wri_vel( 'velhradcp', velodata(jveloset), n2dint )1352 1353 ENDIF1354 1355 ! Feedback velocity data1356 1357 IF ( ln_velfb ) THEN1358 1359 jfbini = jveloset + 11360 1361 DO jveloset = jfbini, nvelosets1362 1363 jset = jveloset - jfbini + 11364 WRITE(cdtmp,'(A,I2.2)')'velfb_',jset1365 CALL obs_wri_vel( cdtmp, velodata(jveloset), n2dint )1366 1367 END DO1368 1369 ENDIF1370 1371 798 ENDIF 1372 799 … … 1386 813 !! 1387 814 !!---------------------------------------------------------------------- 1388 ! !obs_grid deallocation815 ! obs_grid deallocation 1389 816 CALL obs_grid_deallocate 1390 817 1391 !! diaobs deallocation 1392 IF ( nprofsets > 0 ) THEN 1393 DEALLOCATE(ld_enact, & 1394 & profdata, & 1395 & prodatqc) 1396 END IF 1397 IF ( ln_sla ) THEN 1398 DEALLOCATE(sladata, & 1399 & sladatqc) 1400 END IF 1401 IF ( ln_seaice ) THEN 1402 DEALLOCATE(sladata, & 1403 & sladatqc) 1404 END IF 1405 IF ( ln_sst ) THEN 1406 DEALLOCATE(sstdata, & 1407 & sstdatqc) 1408 END IF 1409 IF ( ln_vel3d ) THEN 1410 DEALLOCATE(ld_velav, & 1411 & velodata, & 1412 & veldatqc) 1413 END IF 818 ! diaobs deallocation 819 IF ( nproftypes > 0 ) & 820 & DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof ) 821 822 IF ( nsurftypes > 0 ) & 823 & DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf ) 824 1414 825 END SUBROUTINE dia_obs_dealloc 1415 826 1416 SUBROUTINE ini_date( ddobsini)1417 !!---------------------------------------------------------------------- 1418 !! *** ROUTINE ini_date ***827 SUBROUTINE calc_date( kstp, ddobs ) 828 !!---------------------------------------------------------------------- 829 !! *** ROUTINE calc_date *** 1419 830 !! 1420 !! ** Purpose : Get initial data in double precision YYYYMMDD.HHMMSS format 1421 !! 1422 !! ** Method : Get initial data in double precision YYYYMMDD.HHMMSS format 1423 !! 1424 !! ** Action : Get initial data in double precision YYYYMMDD.HHMMSS format 831 !! ** Purpose : Get date in double precision YYYYMMDD.HHMMSS format 832 !! 833 !! ** Method : Get date in double precision YYYYMMDD.HHMMSS format 834 !! 835 !! ** Action : Get date in double precision YYYYMMDD.HHMMSS format 836 !! 837 !! ** Action : Get initial date in double precision YYYYMMDD.HHMMSS format 1425 838 !! 1426 839 !! History : … … 1430 843 !! ! 06-10 (G. Smith) Calculates initial date the same as method for final date 1431 844 !! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2 845 !! ! 2014-09 (D. Lea) New generic routine now deals with arbitrary initial time of day 1432 846 !!---------------------------------------------------------------------- 1433 847 USE phycst, ONLY : & ! Physical constants 1434 848 & rday 1435 ! USE daymod, ONLY : & ! Time variables1436 ! & nmonth_len1437 849 USE dom_oce, ONLY : & ! Ocean space and time domain variables 1438 850 & rdt … … 1441 853 1442 854 !! * Arguments 1443 REAL(KIND=dp), INTENT(OUT) :: ddobsini ! Initial date in YYYYMMDD.HHMMSS 855 REAL(KIND=dp), INTENT(OUT) :: ddobs ! Date in YYYYMMDD.HHMMSS 856 INTEGER :: kstp 1444 857 1445 858 !! * Local declarations … … 1449 862 INTEGER :: ihou 1450 863 INTEGER :: imin 1451 INTEGER :: imday 1452 REAL( KIND=wp) :: zdayfrc! Fraction of day864 INTEGER :: imday ! Number of days in month. 865 REAL(wp) :: zdayfrc ! Fraction of day 1453 866 1454 867 INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year … … 1456 869 !!---------------------------------------------------------------------- 1457 870 !! Initial date initialization (year, month, day, hour, minute) 1458 !! (This assumes that the initial date is for 00z))1459 871 !!---------------------------------------------------------------------- 1460 872 iyea = ndate0 / 10000 1461 873 imon = ( ndate0 - iyea * 10000 ) / 100 1462 874 iday = ndate0 - iyea * 10000 - imon * 100 1463 ihou = 01464 imin = 0875 ihou = nn_time0 / 100 876 imin = ( nn_time0 - ihou * 100 ) 1465 877 1466 878 !!---------------------------------------------------------------------- 1467 879 !! Compute number of days + number of hours + min since initial time 1468 880 !!---------------------------------------------------------------------- 1469 iday = iday + ( nit000 -1 ) * rdt / rday 1470 zdayfrc = ( nit000 -1 ) * rdt / rday 881 zdayfrc = kstp * rdt / rday 1471 882 zdayfrc = zdayfrc - aint(zdayfrc) 1472 ihou = int( zdayfrc * 24 ) 1473 imin = int( (zdayfrc * 24 - ihou) * 60 ) 1474 1475 !!----------------------------------------------------------------------- 1476 !! Convert number of days (iday) into a real date 1477 !!---------------------------------------------------------------------- 883 imin = imin + int( zdayfrc * 24 * 60 ) 884 DO WHILE (imin >= 60) 885 imin=imin-60 886 ihou=ihou+1 887 END DO 888 DO WHILE (ihou >= 24) 889 ihou=ihou-24 890 iday=iday+1 891 END DO 892 iday = iday + kstp * rdt / rday 893 894 !----------------------------------------------------------------------- 895 ! Convert number of days (iday) into a real date 896 !---------------------------------------------------------------------- 1478 897 1479 898 CALL calc_month_len( iyea, imonth_len ) 1480 899 1481 900 DO WHILE ( iday > imonth_len(imon) ) 1482 901 iday = iday - imonth_len(imon) … … 1489 908 END DO 1490 909 1491 !!---------------------------------------------------------------------- 1492 !! Convert it into YYYYMMDD.HHMMSS format. 1493 !!---------------------------------------------------------------------- 1494 ddobsini = iyea * 10000_dp + imon * 100_dp + & 1495 & iday + ihou * 0.01_dp + imin * 0.0001_dp 1496 1497 1498 END SUBROUTINE ini_date 1499 1500 SUBROUTINE fin_date( ddobsfin ) 1501 !!---------------------------------------------------------------------- 1502 !! *** ROUTINE fin_date *** 910 !---------------------------------------------------------------------- 911 ! Convert it into YYYYMMDD.HHMMSS format. 912 !---------------------------------------------------------------------- 913 ddobs = iyea * 10000_dp + imon * 100_dp + & 914 & iday + ihou * 0.01_dp + imin * 0.0001_dp 915 916 END SUBROUTINE calc_date 917 918 SUBROUTINE ini_date( ddobsini ) 919 !!---------------------------------------------------------------------- 920 !! *** ROUTINE ini_date *** 1503 921 !! 1504 !! ** Purpose : Get final datain double precision YYYYMMDD.HHMMSS format1505 !! 1506 !! ** Method : Get final data in double precision YYYYMMDD.HHMMSS format1507 !! 1508 !! ** Action : Get final data in double precision YYYYMMDD.HHMMSS format922 !! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format 923 !! 924 !! ** Method : 925 !! 926 !! ** Action : 1509 927 !! 1510 928 !! History : … … 1513 931 !! ! 06-10 (A. Weaver) Cleaning 1514 932 !! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2 1515 !!---------------------------------------------------------------------- 1516 USE phycst, ONLY : & ! Physical constants 1517 & rday 1518 ! USE daymod, ONLY : & ! Time variables 1519 ! & nmonth_len 1520 USE dom_oce, ONLY : & ! Ocean space and time domain variables 1521 & rdt 933 !! ! 2014-09 (D. Lea) Change to call generic routine calc_date 934 !!---------------------------------------------------------------------- 1522 935 1523 936 IMPLICIT NONE 1524 937 1525 938 !! * Arguments 1526 REAL(KIND=dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS 1527 1528 !! * Local declarations 1529 INTEGER :: iyea ! date - (year, month, day, hour, minute) 1530 INTEGER :: imon 1531 INTEGER :: iday 1532 INTEGER :: ihou 1533 INTEGER :: imin 1534 INTEGER :: imday ! Number of days in month. 1535 REAL(KIND=wp) :: zdayfrc ! Fraction of day 1536 1537 INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year 1538 1539 !----------------------------------------------------------------------- 1540 ! Initial date initialization (year, month, day, hour, minute) 1541 ! (This assumes that the initial date is for 00z) 1542 !----------------------------------------------------------------------- 1543 iyea = ndate0 / 10000 1544 imon = ( ndate0 - iyea * 10000 ) / 100 1545 iday = ndate0 - iyea * 10000 - imon * 100 1546 ihou = 0 1547 imin = 0 1548 1549 !----------------------------------------------------------------------- 1550 ! Compute number of days + number of hours + min since initial time 1551 !----------------------------------------------------------------------- 1552 iday = iday + nitend * rdt / rday 1553 zdayfrc = nitend * rdt / rday 1554 zdayfrc = zdayfrc - AINT( zdayfrc ) 1555 ihou = INT( zdayfrc * 24 ) 1556 imin = INT( ( zdayfrc * 24 - ihou ) * 60 ) 1557 1558 !----------------------------------------------------------------------- 1559 ! Convert number of days (iday) into a real date 1560 !---------------------------------------------------------------------- 1561 1562 CALL calc_month_len( iyea, imonth_len ) 1563 1564 DO WHILE ( iday > imonth_len(imon) ) 1565 iday = iday - imonth_len(imon) 1566 imon = imon + 1 1567 IF ( imon > 12 ) THEN 1568 imon = 1 1569 iyea = iyea + 1 1570 CALL calc_month_len( iyea, imonth_len ) ! update month lengths 1571 ENDIF 1572 END DO 1573 1574 !----------------------------------------------------------------------- 1575 ! Convert it into YYYYMMDD.HHMMSS format 1576 !----------------------------------------------------------------------- 1577 ddobsfin = iyea * 10000_dp + imon * 100_dp + iday & 1578 & + ihou * 0.01_dp + imin * 0.0001_dp 1579 1580 END SUBROUTINE fin_date 1581 939 REAL(KIND=dp), INTENT(OUT) :: ddobsini ! Initial date in YYYYMMDD.HHMMSS 940 941 CALL calc_date( nit000 - 1, ddobsini ) 942 943 END SUBROUTINE ini_date 944 945 SUBROUTINE fin_date( ddobsfin ) 946 !!---------------------------------------------------------------------- 947 !! *** ROUTINE fin_date *** 948 !! 949 !! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format 950 !! 951 !! ** Method : 952 !! 953 !! ** Action : 954 !! 955 !! History : 956 !! ! 06-03 (K. Mogensen) Original code 957 !! ! 06-05 (K. Mogensen) Reformatted 958 !! ! 06-10 (A. Weaver) Cleaning 959 !! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2 960 !! ! 2014-09 (D. Lea) Change to call generic routine calc_date 961 !!---------------------------------------------------------------------- 962 963 IMPLICIT NONE 964 965 !! * Arguments 966 REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS 967 968 CALL calc_date( nitend, ddobsfin ) 969 970 END SUBROUTINE fin_date 971 1582 972 END MODULE diaobs
Note: See TracChangeset
for help on using the changeset viewer.
