MODULE diaobs !!====================================================================== !! *** MODULE diaobs *** !! Observation diagnostics: Computation of the misfit between data and !! their model equivalent !!====================================================================== !!---------------------------------------------------------------------- !! dia_obs_init : Reading and prepare observations !! dia_obs : Compute model equivalent to observations !! dia_obs_wri : Write observational diagnostics !! ini_date : Compute the initial date YYYYMMDD.HHMMSS !! fin_date : Compute the final date YYYYMMDD.HHMMSS !!---------------------------------------------------------------------- !! * Modules used USE wrk_nemo ! Memory Allocation USE par_kind ! Precision variables USE in_out_manager ! I/O manager USE par_oce USE dom_oce ! Ocean space and time domain variables USE obs_read_prof ! Reading and allocation of profile obs USE obs_read_surf ! Reading and allocation of surface obs USE obs_readmdt ! Reading and allocation of MDT for SLA. USE obs_prep ! Preparation of obs. (grid search etc). USE obs_oper ! Observation operators USE obs_write ! Writing of observation related diagnostics USE obs_grid ! Grid searching USE obs_read_altbias ! Bias treatment for altimeter USE obs_sstbias ! Bias correction routine for SST USE obs_profiles_def ! Profile data definitions USE obs_surf_def ! Surface data definitions USE obs_types ! Definitions for observation types USE mpp_map ! MPP mapping USE lib_mpp ! For ctl_warn/stop IMPLICIT NONE !! * Routine accessibility PRIVATE PUBLIC dia_obs_init, & ! Initialize and read observations & dia_obs, & ! Compute model equivalent to observations & dia_obs_wri, & ! Write model equivalent to observations & dia_obs_dealloc ! Deallocate dia_obs data !! * Module variables LOGICAL, PUBLIC :: & & lk_diaobs = .TRUE. !: Include this for backwards compatibility at NEMO 3.6. LOGICAL :: ln_diaobs !: Logical switch for the obs operator LOGICAL :: ln_sstnight !: Logical switch for night mean SST obs LOGICAL :: ln_sla_fp_indegs !: T=> SLA obs footprint size specified in degrees, F=> in metres LOGICAL :: ln_sst_fp_indegs !: T=> SST obs footprint size specified in degrees, F=> in metres LOGICAL :: ln_sss_fp_indegs !: T=> SSS obs footprint size specified in degrees, F=> in metres LOGICAL :: ln_sic_fp_indegs !: T=> sea-ice obs footprint size specified in degrees, F=> in metres REAL(wp) :: rn_sla_avglamscl !: E/W diameter of SLA observation footprint (metres) REAL(wp) :: rn_sla_avgphiscl !: N/S diameter of SLA observation footprint (metres) REAL(wp) :: rn_sst_avglamscl !: E/W diameter of SST observation footprint (metres) REAL(wp) :: rn_sst_avgphiscl !: N/S diameter of SST observation footprint (metres) REAL(wp) :: rn_sss_avglamscl !: E/W diameter of SSS observation footprint (metres) REAL(wp) :: rn_sss_avgphiscl !: N/S diameter of SSS observation footprint (metres) REAL(wp) :: rn_sic_avglamscl !: E/W diameter of sea-ice observation footprint (metres) REAL(wp) :: rn_sic_avgphiscl !: N/S diameter of sea-ice observation footprint (metres) INTEGER :: nn_1dint !: Vertical interpolation method INTEGER :: nn_2dint !: Default horizontal interpolation method INTEGER :: nn_2dint_sla !: SLA horizontal interpolation method INTEGER :: nn_2dint_sst !: SST horizontal interpolation method INTEGER :: nn_2dint_sss !: SSS horizontal interpolation method INTEGER :: nn_2dint_sic !: Seaice horizontal interpolation method INTEGER, DIMENSION(imaxavtypes) :: & & nn_profdavtypes !: Profile data types representing a daily average INTEGER :: nproftypes !: Number of profile obs types INTEGER :: nsurftypes !: Number of surface obs types INTEGER, DIMENSION(:), ALLOCATABLE :: & & nvarsprof, & !: Number of profile variables & nvarssurf !: Number of surface variables INTEGER, DIMENSION(:), ALLOCATABLE :: & & nextrprof, & !: Number of profile extra variables & nextrsurf !: Number of surface extra variables INTEGER, DIMENSION(:), ALLOCATABLE :: & & n2dintsurf !: Interpolation option for surface variables REAL(wp), DIMENSION(:), ALLOCATABLE :: & & ravglamscl, & !: E/W diameter of averaging footprint for surface variables & ravgphiscl !: N/S diameter of averaging footprint for surface variables LOGICAL, DIMENSION(:), ALLOCATABLE :: & & lfpindegs, & !: T=> surface obs footprint size specified in degrees, F=> in metres & llnightav !: Logical for calculating night-time averages TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: & & surfdata, & !: Initial surface data & surfdataqc !: Surface data after quality control TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: & & profdata, & !: Initial profile data & profdataqc !: Profile data after quality control CHARACTER(len=6), PUBLIC, DIMENSION(:), ALLOCATABLE :: & & cobstypesprof, & !: Profile obs types & cobstypessurf !: Surface obs types !!---------------------------------------------------------------------- !! NEMO/OPA 3.3 , NEMO Consortium (2010) !! $Id$ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- !! * Substitutions # include "domzgr_substitute.h90" CONTAINS SUBROUTINE dia_obs_init !!---------------------------------------------------------------------- !! *** ROUTINE dia_obs_init *** !! !! ** Purpose : Initialize and read observations !! !! ** Method : Read the namelist and call reading routines !! !! ** Action : Read the namelist and call reading routines !! !! History : !! ! 06-03 (K. Mogensen) Original code !! ! 06-05 (A. Weaver) Reformatted !! ! 06-10 (A. Weaver) Cleaning and add controls !! ! 07-03 (K. Mogensen) General handling of profiles !! ! 14-08 (J.While) Incorporated SST bias correction !! ! 15-02 (M. Martin) Simplification of namelist and code !!---------------------------------------------------------------------- IMPLICIT NONE !! * Local declarations INTEGER, PARAMETER :: & & jpmaxnfiles = 1000 ! Maximum number of files for each obs type INTEGER, DIMENSION(:), ALLOCATABLE :: & & ifilesprof, & ! Number of profile files & ifilessurf ! Number of surface files INTEGER :: ios ! Local integer output status for namelist read INTEGER :: jtype ! Counter for obs types INTEGER :: jvar ! Counter for variables INTEGER :: jfile ! Counter for files INTEGER :: jnumsstbias ! Number of SST bias files to read and apply CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: & & cn_profbfiles, & ! T/S profile input filenames & cn_sstfbfiles, & ! Sea surface temperature input filenames & cn_slafbfiles, & ! Sea level anomaly input filenames & cn_sicfbfiles, & ! Seaice concentration input filenames & cn_velfbfiles, & ! Velocity profile input filenames & cn_sssfbfiles, & ! Sea surface salinity input filenames & cn_logchlfbfiles, & ! Log(Chl) input filenames & cn_spmfbfiles, & ! Sediment input filenames & cn_fco2fbfiles, & ! fco2 input filenames & cn_pco2fbfiles, & ! pco2 input filenames & cn_sstbiasfiles ! SST bias input filenames CHARACTER(LEN=128) :: & & cn_altbiasfile ! Altimeter bias input filename LOGICAL :: ln_t3d ! Logical switch for temperature profiles LOGICAL :: ln_s3d ! Logical switch for salinity profiles LOGICAL :: ln_sla ! Logical switch for sea level anomalies LOGICAL :: ln_sst ! Logical switch for sea surface temperature LOGICAL :: ln_sic ! Logical switch for sea ice concentration LOGICAL :: ln_sss ! Logical switch for sea surface salinity obs LOGICAL :: ln_vel3d ! Logical switch for velocity (u,v) obs LOGICAL :: ln_logchl ! Logical switch for log(Chl) obs LOGICAL :: ln_spm ! Logical switch for sediment obs LOGICAL :: ln_fco2 ! Logical switch for fco2 obs LOGICAL :: ln_pco2 ! Logical switch for pco2 obs LOGICAL :: ln_nea ! Logical switch to remove obs near land LOGICAL :: ln_altbias ! Logical switch for altimeter bias LOGICAL :: ln_sstbias ! Logical switch for bias correction of SST LOGICAL :: ln_ignmis ! Logical switch for ignoring missing files LOGICAL :: ln_s_at_t ! Logical switch to compute model S at T obs LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary REAL(dp) :: rn_dobsini ! Obs window start date YYYYMMDD.HHMMSS REAL(dp) :: rn_dobsend ! Obs window end date YYYYMMDD.HHMMSS CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: & & clproffiles, & ! Profile filenames & clsurffiles ! Surface filenames LOGICAL :: llvar1 ! Logical for profile variable 1 LOGICAL :: llvar2 ! Logical for profile variable 1 REAL(wp), POINTER, DIMENSION(:,:) :: & & zglam1, & ! Model longitudes for profile variable 1 & zglam2 ! Model longitudes for profile variable 2 REAL(wp), POINTER, DIMENSION(:,:) :: & & zgphi1, & ! Model latitudes for profile variable 1 & zgphi2 ! Model latitudes for profile variable 2 REAL(wp), POINTER, DIMENSION(:,:,:) :: & & zmask1, & ! Model land/sea mask associated with variable 1 & zmask2 ! Model land/sea mask associated with variable 2 NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla, & & ln_sst, ln_sic, ln_sss, ln_vel3d, & & ln_logchl, ln_spm, ln_fco2, ln_pco2, & & ln_altbias, ln_sstbias, ln_nea, & & ln_grid_global, ln_grid_search_lookup, & & ln_ignmis, ln_s_at_t, ln_bound_reject, & & ln_sstnight, & & ln_sla_fp_indegs, ln_sst_fp_indegs, & & ln_sss_fp_indegs, ln_sic_fp_indegs, & & cn_profbfiles, cn_slafbfiles, & & cn_sstfbfiles, cn_sicfbfiles, & & cn_velfbfiles, cn_sssfbfiles, & & cn_logchlfbfiles, cn_spmfbfiles, & & cn_fco2fbfiles, cn_pco2fbfiles, & & cn_sstbiasfiles, cn_altbiasfile, & & cn_gridsearchfile, rn_gridsearchres, & & rn_dobsini, rn_dobsend, & & rn_sla_avglamscl, rn_sla_avgphiscl, & & rn_sst_avglamscl, rn_sst_avgphiscl, & & rn_sss_avglamscl, rn_sss_avgphiscl, & & rn_sic_avglamscl, rn_sic_avgphiscl, & & nn_1dint, nn_2dint, & & nn_2dint_sla, nn_2dint_sst, & & nn_2dint_sss, nn_2dint_sic, & & nn_msshc, rn_mdtcorr, rn_mdtcutoff, & & nn_profdavtypes CALL wrk_alloc( jpi, jpj, zglam1 ) CALL wrk_alloc( jpi, jpj, zglam2 ) CALL wrk_alloc( jpi, jpj, zgphi1 ) CALL wrk_alloc( jpi, jpj, zgphi2 ) CALL wrk_alloc( jpi, jpj, jpk, zmask1 ) CALL wrk_alloc( jpi, jpj, jpk, zmask2 ) !----------------------------------------------------------------------- ! Read namelist parameters !----------------------------------------------------------------------- ! Some namelist arrays need initialising cn_profbfiles(:) = '' cn_slafbfiles(:) = '' cn_sstfbfiles(:) = '' cn_sicfbfiles(:) = '' cn_velfbfiles(:) = '' cn_sssfbfiles(:) = '' cn_logchlfbfiles(:) = '' cn_spmfbfiles(:) = '' cn_fco2fbfiles(:) = '' cn_pco2fbfiles(:) = '' cn_sstbiasfiles(:) = '' nn_profdavtypes(:) = -1 CALL ini_date( rn_dobsini ) CALL fin_date( rn_dobsend ) ! Read namelist namobs : control observation diagnostics REWIND( numnam_ref ) ! Namelist namobs in reference namelist READ ( numnam_ref, namobs, IOSTAT = ios, ERR = 901) 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp ) REWIND( numnam_cfg ) ! Namelist namobs in configuration namelist READ ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 ) 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp ) IF(lwm) WRITE ( numond, namobs ) lk_diaobs = .FALSE. #if defined key_diaobs IF ( ln_diaobs ) lk_diaobs = .TRUE. #endif IF ( .NOT. lk_diaobs ) THEN IF(lwp) WRITE(numout,cform_war) IF(lwp) WRITE(numout,*)' ln_diaobs is set to false or key_diaobs is not set, so not calling dia_obs' RETURN ENDIF IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'dia_obs_init : Observation diagnostic initialization' WRITE(numout,*) '~~~~~~~~~~~~' WRITE(numout,*) ' Namelist namobs : set observation diagnostic parameters' WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst WRITE(numout,*) ' Logical switch for Sea Ice observations ln_sic = ', ln_sic WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss WRITE(numout,*) ' Logical switch for log(Chl) observations ln_logchl = ', ln_logchl WRITE(numout,*) ' Logical switch for SPM observations ln_spm = ', ln_spm WRITE(numout,*) ' Logical switch for FCO2 observations ln_fco2 = ', ln_fco2 WRITE(numout,*) ' Logical switch for PCO2 observations ln_pco2 = ', ln_pco2 WRITE(numout,*) ' Global distribution of observations ln_grid_global = ', ln_grid_global WRITE(numout,*) ' Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup IF (ln_grid_search_lookup) & WRITE(numout,*) ' Grid search lookup file header cn_gridsearchfile = ', cn_gridsearchfile WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS rn_dobsini = ', rn_dobsini WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS rn_dobsend = ', rn_dobsend WRITE(numout,*) ' Type of vertical interpolation method nn_1dint = ', nn_1dint WRITE(numout,*) ' Type of horizontal interpolation method nn_2dint = ', nn_2dint WRITE(numout,*) ' Rejection of observations near land switch ln_nea = ', ln_nea WRITE(numout,*) ' Rejection of obs near open bdys ln_bound_reject = ', ln_bound_reject WRITE(numout,*) ' MSSH correction scheme nn_msshc = ', nn_msshc WRITE(numout,*) ' MDT correction rn_mdtcorr = ', rn_mdtcorr WRITE(numout,*) ' MDT cutoff for computed correction rn_mdtcutoff = ', rn_mdtcutoff WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias WRITE(numout,*) ' Logical switch for sst bias ln_sstbias = ', ln_sstbias WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis WRITE(numout,*) ' Daily average types nn_profdavtypes = ', nn_profdavtypes WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight ENDIF !----------------------------------------------------------------------- ! Set up list of observation types to be used ! and the files associated with each type !----------------------------------------------------------------------- nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss, & & ln_logchl, ln_spm, ln_fco2, ln_pco2 /) ) IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN IF(lwp) WRITE(numout,cform_war) IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & & ' are set to .FALSE. so turning off calls to dia_obs' nwarn = nwarn + 1 lk_diaobs = .FALSE. RETURN ENDIF IF(lwp) WRITE(numout,*) ' Number of profile obs types: ',nproftypes IF ( nproftypes > 0 ) THEN ALLOCATE( cobstypesprof(nproftypes) ) ALLOCATE( ifilesprof(nproftypes) ) ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) jtype = 0 IF (ln_t3d .OR. ln_s3d) THEN jtype = jtype + 1 CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof ', & & cn_profbfiles, ifilesprof, cobstypesprof, clproffiles ) ENDIF IF (ln_vel3d) THEN jtype = jtype + 1 CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel ', & & cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles ) ENDIF ENDIF IF(lwp) WRITE(numout,*)' Number of surface obs types: ',nsurftypes IF ( nsurftypes > 0 ) THEN ALLOCATE( cobstypessurf(nsurftypes) ) ALLOCATE( ifilessurf(nsurftypes) ) ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) ALLOCATE(n2dintsurf(nsurftypes)) ALLOCATE(ravglamscl(nsurftypes)) ALLOCATE(ravgphiscl(nsurftypes)) ALLOCATE(lfpindegs(nsurftypes)) ALLOCATE(llnightav(nsurftypes)) jtype = 0 IF (ln_sla) THEN jtype = jtype + 1 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sla ', & & cn_slafbfiles, ifilessurf, cobstypessurf, clsurffiles ) CALL obs_setinterpopts( nsurftypes, jtype, 'sla ', & & nn_2dint, nn_2dint_sla, & & rn_sla_avglamscl, rn_sla_avgphiscl, & & ln_sla_fp_indegs, .FALSE., & & n2dintsurf, ravglamscl, ravgphiscl, & & lfpindegs, llnightav ) ENDIF IF (ln_sst) THEN jtype = jtype + 1 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sst ', & & cn_sstfbfiles, ifilessurf, cobstypessurf, clsurffiles ) CALL obs_setinterpopts( nsurftypes, jtype, 'sst ', & & nn_2dint, nn_2dint_sst, & & rn_sst_avglamscl, rn_sst_avgphiscl, & & ln_sst_fp_indegs, ln_sstnight, & & n2dintsurf, ravglamscl, ravgphiscl, & & lfpindegs, llnightav ) ENDIF #if defined key_lim2 || defined key_lim3 || defined key_cice IF (ln_sic) THEN jtype = jtype + 1 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sic ', & & cn_sicfbfiles, ifilessurf, cobstypessurf, clsurffiles ) CALL obs_setinterpopts( nsurftypes, jtype, 'sic ', & & nn_2dint, nn_2dint_sic, & & rn_sic_avglamscl, rn_sic_avgphiscl, & & ln_sic_fp_indegs, .FALSE., & & n2dintsurf, ravglamscl, ravgphiscl, & & lfpindegs, llnightav ) ENDIF #endif IF (ln_sss) THEN jtype = jtype + 1 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sss ', & & cn_sssfbfiles, ifilessurf, cobstypessurf, clsurffiles ) CALL obs_setinterpopts( nsurftypes, jtype, 'sss ', & & nn_2dint, nn_2dint_sss, & & rn_sss_avglamscl, rn_sss_avgphiscl, & & ln_sss_fp_indegs, .FALSE., & & n2dintsurf, ravglamscl, ravgphiscl, & & lfpindegs, llnightav ) ENDIF IF (ln_logchl) THEN jtype = jtype + 1 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'logchl', & & cn_logchlfbfiles, ifilessurf, cobstypessurf, clsurffiles ) CALL obs_setinterpopts( nsurftypes, jtype, 'logchl', & & nn_2dint, -1, 0., 0., .TRUE., .FALSE., & & n2dintsurf, ravglamscl, ravgphiscl, & & lfpindegs, llnightav ) ENDIF IF (ln_spm) THEN jtype = jtype + 1 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'spm ', & & cn_spmfbfiles, ifilessurf, cobstypessurf, clsurffiles ) CALL obs_setinterpopts( nsurftypes, jtype, 'spm ', & & nn_2dint, -1, 0., 0., .TRUE., .FALSE., & & n2dintsurf, ravglamscl, ravgphiscl, & & lfpindegs, llnightav ) ENDIF IF (ln_fco2) THEN jtype = jtype + 1 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'fco2 ', & & cn_fco2fbfiles, ifilessurf, cobstypessurf, clsurffiles ) CALL obs_setinterpopts( nsurftypes, jtype, 'fco2 ', & & nn_2dint, -1, 0., 0., .TRUE., .FALSE., & & n2dintsurf, ravglamscl, ravgphiscl, & & lfpindegs, llnightav ) ENDIF IF (ln_pco2) THEN jtype = jtype + 1 CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'pco2 ', & & cn_pco2fbfiles, ifilessurf, cobstypessurf, clsurffiles ) CALL obs_setinterpopts( nsurftypes, jtype, 'pco2 ', & & nn_2dint, -1, 0., 0., .TRUE., .FALSE., & & n2dintsurf, ravglamscl, ravgphiscl, & & lfpindegs, llnightav ) ENDIF ENDIF IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' !----------------------------------------------------------------------- ! Obs operator parameter checking and initialisations !----------------------------------------------------------------------- IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' ) RETURN ENDIF IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN CALL ctl_stop(' Choice of vertical (1D) interpolation method', & & ' is not available') ENDIF IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 6 ) ) THEN CALL ctl_stop(' Choice of horizontal (2D) interpolation method', & & ' is not available') ENDIF CALL obs_typ_init CALL mppmap_init CALL obs_grid_setup( ) !----------------------------------------------------------------------- ! Depending on switches read the various observation types !----------------------------------------------------------------------- IF ( nproftypes > 0 ) THEN ALLOCATE(profdata(nproftypes)) ALLOCATE(profdataqc(nproftypes)) ALLOCATE(nvarsprof(nproftypes)) ALLOCATE(nextrprof(nproftypes)) DO jtype = 1, nproftypes nvarsprof(jtype) = 2 IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN nextrprof(jtype) = 1 llvar1 = ln_t3d llvar2 = ln_s3d zglam1 = glamt zgphi1 = gphit zmask1 = tmask zglam2 = glamt zgphi2 = gphit zmask2 = tmask ENDIF IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN nextrprof(jtype) = 2 llvar1 = ln_vel3d llvar2 = ln_vel3d zglam1 = glamu zgphi1 = gphiu zmask1 = umask zglam2 = glamv zgphi2 = gphiv zmask2 = vmask ENDIF !Read in profile or profile obs types CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype), & & clproffiles(jtype,1:ifilesprof(jtype)), & & nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, & & rn_dobsini, rn_dobsend, llvar1, llvar2, & & ln_ignmis, ln_s_at_t, .FALSE., & & kdailyavtypes = nn_profdavtypes ) DO jvar = 1, nvarsprof(jtype) CALL obs_prof_staend( profdata(jtype), jvar ) END DO CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), & & llvar1, llvar2, & & jpi, jpj, jpk, & & zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2, & & ln_nea, ln_bound_reject, & & kdailyavtypes = nn_profdavtypes ) END DO DEALLOCATE( ifilesprof, clproffiles ) ENDIF IF ( nsurftypes > 0 ) THEN ALLOCATE(surfdata(nsurftypes)) ALLOCATE(surfdataqc(nsurftypes)) ALLOCATE(nvarssurf(nsurftypes)) ALLOCATE(nextrsurf(nsurftypes)) DO jtype = 1, nsurftypes nvarssurf(jtype) = 1 nextrsurf(jtype) = 0 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2 !Read in surface obs types CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), & & clsurffiles(jtype,1:ifilessurf(jtype)), & & nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & & rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) ) CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject ) IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) ) IF ( ln_altbias ) & & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile ) ENDIF IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN jnumsstbias = 0 DO jfile = 1, jpmaxnfiles IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) & & jnumsstbias = jnumsstbias + 1 END DO IF ( jnumsstbias == 0 ) THEN CALL ctl_stop("ln_sstbias set but no bias files to read in") ENDIF CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), & & jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) ) ENDIF END DO DEALLOCATE( ifilessurf, clsurffiles ) ENDIF CALL wrk_dealloc( jpi, jpj, zglam1 ) CALL wrk_dealloc( jpi, jpj, zglam2 ) CALL wrk_dealloc( jpi, jpj, zgphi1 ) CALL wrk_dealloc( jpi, jpj, zgphi2 ) CALL wrk_dealloc( jpi, jpj, jpk, zmask1 ) CALL wrk_dealloc( jpi, jpj, jpk, zmask2 ) END SUBROUTINE dia_obs_init SUBROUTINE dia_obs( kstp ) !!---------------------------------------------------------------------- !! *** ROUTINE dia_obs *** !! !! ** Purpose : Call the observation operators on each time step !! !! ** Method : Call the observation operators on each time step to !! compute the model equivalent of the following data: !! - Profile data, currently T/S or U/V !! - Surface data, currently SST, SLA or sea-ice concentration. !! !! ** Action : !! !! History : !! ! 06-03 (K. Mogensen) Original code !! ! 06-05 (K. Mogensen) Reformatted !! ! 06-10 (A. Weaver) Cleaning !! ! 07-03 (K. Mogensen) General handling of profiles !! ! 07-04 (G. Smith) Generalized surface operators !! ! 08-10 (M. Valdivieso) obs operator for velocity profiles !! ! 15-08 (M. Martin) Combined surface/profile routines. !!---------------------------------------------------------------------- !! * Modules used USE phycst, ONLY : & ! Physical constants & rday USE oce, ONLY : & ! Ocean dynamics and tracers variables & tsn, & & un, & & vn, & & sshn #if defined key_lim3 USE ice, ONLY : & ! LIM3 Ice model variables & frld #endif #if defined key_lim2 USE ice_2, ONLY : & ! LIM2 Ice model variables & frld #endif #if defined key_cice USE sbc_oce, ONLY : fr_i ! ice fraction #endif #if defined key_hadocc USE trc, ONLY : & ! HadOCC chlorophyll, fCO2 and pCO2 & HADOCC_CHL, & & HADOCC_FCO2, & & HADOCC_PCO2, & & HADOCC_FILL_FLT #elif defined key_medusa && defined key_foam_medusa USE trc, ONLY : & ! MEDUSA chlorophyll, fCO2 and pCO2 & trn USE par_medusa, ONLY: & & jpchn, & & jpchd #if defined key_roam USE sms_medusa, ONLY: & & f2_pco2w, & & f2_fco2w #endif #elif defined key_fabm USE fabm USE par_fabm #endif #if defined key_spm USE par_spm, ONLY: & ! ERSEM/SPM sediments & jp_spm USE trc, ONLY : & & trn #endif IMPLICIT NONE !! * Arguments INTEGER, INTENT(IN) :: kstp ! Current timestep !! * Local declarations INTEGER :: idaystp ! Number of timesteps per day INTEGER :: jtype ! Data loop variable INTEGER :: jvar ! Variable number INTEGER :: ji, jj ! Loop counters REAL(wp) :: tiny ! small number REAL(wp), POINTER, DIMENSION(:,:,:) :: & & zprofvar1, & ! Model values for 1st variable in a prof ob & zprofvar2 ! Model values for 2nd variable in a prof ob REAL(wp), POINTER, DIMENSION(:,:,:) :: & & zprofmask1, & ! Mask associated with zprofvar1 & zprofmask2 ! Mask associated with zprofvar2 REAL(wp), POINTER, DIMENSION(:,:) :: & & zsurfvar, & ! Model values equivalent to surface ob. & zsurfmask ! Mask associated with surface variable REAL(wp), POINTER, DIMENSION(:,:) :: & & zglam1, & ! Model longitudes for prof variable 1 & zglam2, & ! Model longitudes for prof variable 2 & zgphi1, & ! Model latitudes for prof variable 1 & zgphi2 ! Model latitudes for prof variable 2 !Allocate local work arrays CALL wrk_alloc( jpi, jpj, jpk, zprofvar1 ) CALL wrk_alloc( jpi, jpj, jpk, zprofvar2 ) CALL wrk_alloc( jpi, jpj, jpk, zprofmask1 ) CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 ) CALL wrk_alloc( jpi, jpj, zsurfvar ) CALL wrk_alloc( jpi, jpj, zsurfmask ) CALL wrk_alloc( jpi, jpj, zglam1 ) CALL wrk_alloc( jpi, jpj, zglam2 ) CALL wrk_alloc( jpi, jpj, zgphi1 ) CALL wrk_alloc( jpi, jpj, zgphi2 ) IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'dia_obs : Call the observation operators', kstp WRITE(numout,*) '~~~~~~~' CALL FLUSH(numout) ENDIF idaystp = NINT( rday / rdt ) !----------------------------------------------------------------------- ! Call the profile and surface observation operators !----------------------------------------------------------------------- IF ( nproftypes > 0 ) THEN DO jtype = 1, nproftypes SELECT CASE ( TRIM(cobstypesprof(jtype)) ) CASE('prof') zprofvar1(:,:,:) = tsn(:,:,:,jp_tem) zprofvar2(:,:,:) = tsn(:,:,:,jp_sal) zprofmask1(:,:,:) = tmask(:,:,:) zprofmask2(:,:,:) = tmask(:,:,:) zglam1(:,:) = glamt(:,:) zglam2(:,:) = glamt(:,:) zgphi1(:,:) = gphit(:,:) zgphi2(:,:) = gphit(:,:) CASE('vel') zprofvar1(:,:,:) = un(:,:,:) zprofvar2(:,:,:) = vn(:,:,:) zprofmask1(:,:,:) = umask(:,:,:) zprofmask2(:,:,:) = vmask(:,:,:) zglam1(:,:) = glamu(:,:) zglam2(:,:) = glamv(:,:) zgphi1(:,:) = gphiu(:,:) zgphi2(:,:) = gphiv(:,:) CASE DEFAULT CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) END SELECT CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk, & & nit000, idaystp, & & zprofvar1, zprofvar2, & & fsdept(:,:,:), fsdepw(:,:,:), & & zprofmask1, zprofmask2, & & zglam1, zglam2, zgphi1, zgphi2, & & nn_1dint, nn_2dint, & & kdailyavtypes = nn_profdavtypes ) END DO ENDIF IF ( nsurftypes > 0 ) THEN DO jtype = 1, nsurftypes !Defaults which might be changed zsurfmask(:,:) = tmask(:,:,1) SELECT CASE ( TRIM(cobstypessurf(jtype)) ) CASE('sst') zsurfvar(:,:) = tsn(:,:,1,jp_tem) CASE('sla') zsurfvar(:,:) = sshn(:,:) CASE('sss') zsurfvar(:,:) = tsn(:,:,1,jp_sal) CASE('sic') IF ( kstp == 0 ) THEN IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN CALL ctl_warn( 'Sea-ice not initialised on zeroth '// & & 'time-step but some obs are valid then.' ) WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & & ' sea-ice obs will be missed' ENDIF surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + & & surfdataqc(jtype)%nsstp(1) CYCLE ELSE #if defined key_cice zsurfvar(:,:) = fr_i(:,:) #elif defined key_lim2 || defined key_lim3 zsurfvar(:,:) = 1._wp - frld(:,:) #else CALL ctl_stop( ' Trying to run sea-ice observation operator', & & ' but no sea-ice model appears to have been defined' ) #endif ENDIF CASE('logchl') #if defined key_hadocc zsurfvar(:,:) = HADOCC_CHL(:,:,1) ! (not log) chlorophyll from HadOCC #elif defined key_medusa && defined key_foam_medusa ! Add non-diatom and diatom surface chlorophyll from MEDUSA zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) #elif defined key_fabm chl_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabmdia_chltot) zsurfvar(:,:) = chl_3d(:,:,1) #else CALL ctl_stop( ' Trying to run logchl observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif zsurfmask(:,:) = tmask(:,:,1) ! create a special mask to exclude certain things ! Take the log10 where we can, otherwise exclude tiny = 1.0e-20 WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt ) zsurfvar(:,:) = LOG10(zsurfvar(:,:)) ELSEWHERE zsurfvar(:,:) = obfillflt zsurfmask(:,:) = 0 END WHERE CASE('spm') #if defined key_spm zsurfvar(:,:) = 0.0 DO jn = 1, jp_spm zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn) ! sum SPM sizes END DO #else CALL ctl_stop( ' Trying to run spm observation operator', & & ' but no spm model appears to have been defined' ) #endif CASE('fco2') #if defined key_hadocc zsurfvar(:,:) = HADOCC_FCO2(:,:) ! fCO2 from HadOCC IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. & & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN zsurfvar(:,:) = obfillflt zsurfmask(:,:) = 0 CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', & & ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' ) ENDIF #elif defined key_medusa && defined key_foam_medusa && defined key_roam zsurfvar(:,:) = f2_fco2w(:,:) #elif defined key_fabm ! First, get pCO2 from FABM pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc) zsurfvar(:,:) = pco2_3d(:,:,1) ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of: ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems ! and data reduction routines, Deep-Sea Research II, 56: 512-522. ! and ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas, ! Marine Chemistry, 2: 203-215. ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so ! not explicitly included - atmospheric pressure is not necessarily available so this is ! the best assumption. ! Further, the (1-xCO2)^2 term has been neglected. This is common practice ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes) ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway. zsurfvar(:,:) = zsurfvar(:,:) * EXP((-1636.75 + & & 12.0408 * (tsn(:,:,1,jp_tem)+rt0) - & & 0.0327957 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & & 0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & & 2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0))) / & & (82.0578 * (tsn(:,:,1,jp_tem)+rt0))) #else CALL ctl_stop( ' Trying to run fco2 observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif CASE('pco2') #if defined key_hadocc zsurfvar(:,:) = HADOCC_PCO2(:,:) ! pCO2 from HadOCC IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. & & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN zsurfvar(:,:) = obfillflt zsurfmask(:,:) = 0 CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', & & ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' ) ENDIF #elif defined key_medusa && defined key_foam_medusa && defined key_roam zsurfvar(:,:) = f2_pco2w(:,:) #elif defined key_fabm pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc) zsurfvar(:,:) = pco2_3d(:,:,1) #else CALL ctl_stop( ' Trying to run pCO2 observation operator', & & ' but no biogeochemical model appears to have been defined' ) #endif CASE DEFAULT CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' ) END SELECT CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj, & & nit000, idaystp, zsurfvar, zsurfmask, & & n2dintsurf(jtype), llnightav(jtype), & & ravglamscl(jtype), ravgphiscl(jtype), & & lfpindegs(jtype) ) END DO ENDIF CALL wrk_dealloc( jpi, jpj, jpk, zprofvar1 ) CALL wrk_dealloc( jpi, jpj, jpk, zprofvar2 ) CALL wrk_dealloc( jpi, jpj, jpk, zprofmask1 ) CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 ) CALL wrk_dealloc( jpi, jpj, zsurfvar ) CALL wrk_dealloc( jpi, jpj, zsurfmask ) CALL wrk_dealloc( jpi, jpj, zglam1 ) CALL wrk_dealloc( jpi, jpj, zglam2 ) CALL wrk_dealloc( jpi, jpj, zgphi1 ) CALL wrk_dealloc( jpi, jpj, zgphi2 ) END SUBROUTINE dia_obs SUBROUTINE dia_obs_wri !!---------------------------------------------------------------------- !! *** ROUTINE dia_obs_wri *** !! !! ** Purpose : Call observation diagnostic output routines !! !! ** Method : Call observation diagnostic output routines !! !! ** Action : !! !! History : !! ! 06-03 (K. Mogensen) Original code !! ! 06-05 (K. Mogensen) Reformatted !! ! 06-10 (A. Weaver) Cleaning !! ! 07-03 (K. Mogensen) General handling of profiles !! ! 08-09 (M. Valdivieso) Velocity component (U,V) profiles !! ! 15-08 (M. Martin) Combined writing for prof and surf types !!---------------------------------------------------------------------- !! * Modules used USE obs_rot_vel ! Rotation of velocities IMPLICIT NONE !! * Local declarations INTEGER :: jtype ! Data set loop variable INTEGER :: jo, jvar, jk REAL(wp), DIMENSION(:), ALLOCATABLE :: & & zu, & & zv !----------------------------------------------------------------------- ! Depending on switches call various observation output routines !----------------------------------------------------------------------- IF ( nproftypes > 0 ) THEN DO jtype = 1, nproftypes IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN ! For velocity data, rotate the model velocities to N/S, E/W ! using the compressed data structure. ALLOCATE( & & zu(profdataqc(jtype)%nvprot(1)), & & zv(profdataqc(jtype)%nvprot(2)) & & ) CALL obs_rotvel( profdataqc(jtype), nn_2dint, zu, zv ) DO jo = 1, profdataqc(jtype)%nprof DO jvar = 1, 2 DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar) IF ( jvar == 1 ) THEN profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk) ELSE profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk) ENDIF END DO END DO END DO DEALLOCATE( zu ) DEALLOCATE( zv ) END IF CALL obs_prof_decompress( profdataqc(jtype), & & profdata(jtype), .TRUE., numout ) CALL obs_wri_prof( profdata(jtype) ) END DO ENDIF IF ( nsurftypes > 0 ) THEN DO jtype = 1, nsurftypes CALL obs_surf_decompress( surfdataqc(jtype), & & surfdata(jtype), .TRUE., numout ) CALL obs_wri_surf( surfdata(jtype) ) END DO ENDIF END SUBROUTINE dia_obs_wri SUBROUTINE dia_obs_dealloc IMPLICIT NONE !!---------------------------------------------------------------------- !! *** ROUTINE dia_obs_dealloc *** !! !! ** Purpose : To deallocate data to enable the obs_oper online loop. !! Specifically: dia_obs_init --> dia_obs --> dia_obs_wri !! !! ** Method : Clean up various arrays left behind by the obs_oper. !! !! ** Action : !! !!---------------------------------------------------------------------- ! obs_grid deallocation CALL obs_grid_deallocate ! diaobs deallocation IF ( nproftypes > 0 ) & & DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof ) IF ( nsurftypes > 0 ) & & DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, & & n2dintsurf, ravglamscl, ravgphiscl, lfpindegs, llnightav ) END SUBROUTINE dia_obs_dealloc SUBROUTINE ini_date( ddobsini ) !!---------------------------------------------------------------------- !! *** ROUTINE ini_date *** !! !! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format !! !! ** Method : Get initial date in double precision YYYYMMDD.HHMMSS format !! !! ** Action : Get initial date in double precision YYYYMMDD.HHMMSS format !! !! History : !! ! 06-03 (K. Mogensen) Original code !! ! 06-05 (K. Mogensen) Reformatted !! ! 06-10 (A. Weaver) Cleaning !! ! 06-10 (G. Smith) Calculates initial date the same as method for final date !! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2 !!---------------------------------------------------------------------- USE phycst, ONLY : & ! Physical constants & rday USE dom_oce, ONLY : & ! Ocean space and time domain variables & rdt IMPLICIT NONE !! * Arguments REAL(dp), INTENT(OUT) :: ddobsini ! Initial date in YYYYMMDD.HHMMSS !! * Local declarations INTEGER :: iyea ! date - (year, month, day, hour, minute) INTEGER :: imon INTEGER :: iday INTEGER :: ihou INTEGER :: imin INTEGER :: imday ! Number of days in month. INTEGER, DIMENSION(12) :: & & imonth_len ! Length in days of the months of the current year REAL(wp) :: zdayfrc ! Fraction of day !---------------------------------------------------------------------- ! Initial date initialization (year, month, day, hour, minute) ! (This assumes that the initial date is for 00z)) !---------------------------------------------------------------------- iyea = ndate0 / 10000 imon = ( ndate0 - iyea * 10000 ) / 100 iday = ndate0 - iyea * 10000 - imon * 100 ihou = 0 imin = 0 !---------------------------------------------------------------------- ! Compute number of days + number of hours + min since initial time !---------------------------------------------------------------------- iday = iday + ( nit000 -1 ) * rdt / rday zdayfrc = ( nit000 -1 ) * rdt / rday zdayfrc = zdayfrc - aint(zdayfrc) ihou = int( zdayfrc * 24 ) imin = int( (zdayfrc * 24 - ihou) * 60 ) !----------------------------------------------------------------------- ! Convert number of days (iday) into a real date !---------------------------------------------------------------------- CALL calc_month_len( iyea, imonth_len ) DO WHILE ( iday > imonth_len(imon) ) iday = iday - imonth_len(imon) imon = imon + 1 IF ( imon > 12 ) THEN imon = 1 iyea = iyea + 1 CALL calc_month_len( iyea, imonth_len ) ! update month lengths ENDIF END DO !---------------------------------------------------------------------- ! Convert it into YYYYMMDD.HHMMSS format. !---------------------------------------------------------------------- ddobsini = iyea * 10000_dp + imon * 100_dp + & & iday + ihou * 0.01_dp + imin * 0.0001_dp END SUBROUTINE ini_date SUBROUTINE fin_date( ddobsfin ) !!---------------------------------------------------------------------- !! *** ROUTINE fin_date *** !! !! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format !! !! ** Method : Get final date in double precision YYYYMMDD.HHMMSS format !! !! ** Action : Get final date in double precision YYYYMMDD.HHMMSS format !! !! History : !! ! 06-03 (K. Mogensen) Original code !! ! 06-05 (K. Mogensen) Reformatted !! ! 06-10 (A. Weaver) Cleaning !! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2 !!---------------------------------------------------------------------- USE phycst, ONLY : & ! Physical constants & rday USE dom_oce, ONLY : & ! Ocean space and time domain variables & rdt IMPLICIT NONE !! * Arguments REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS !! * Local declarations INTEGER :: iyea ! date - (year, month, day, hour, minute) INTEGER :: imon INTEGER :: iday INTEGER :: ihou INTEGER :: imin INTEGER :: imday ! Number of days in month. INTEGER, DIMENSION(12) :: & & imonth_len ! Length in days of the months of the current year REAL(wp) :: zdayfrc ! Fraction of day !----------------------------------------------------------------------- ! Initial date initialization (year, month, day, hour, minute) ! (This assumes that the initial date is for 00z) !----------------------------------------------------------------------- iyea = ndate0 / 10000 imon = ( ndate0 - iyea * 10000 ) / 100 iday = ndate0 - iyea * 10000 - imon * 100 ihou = 0 imin = 0 !----------------------------------------------------------------------- ! Compute number of days + number of hours + min since initial time !----------------------------------------------------------------------- iday = iday + nitend * rdt / rday zdayfrc = nitend * rdt / rday zdayfrc = zdayfrc - AINT( zdayfrc ) ihou = INT( zdayfrc * 24 ) imin = INT( ( zdayfrc * 24 - ihou ) * 60 ) !----------------------------------------------------------------------- ! Convert number of days (iday) into a real date !---------------------------------------------------------------------- CALL calc_month_len( iyea, imonth_len ) DO WHILE ( iday > imonth_len(imon) ) iday = iday - imonth_len(imon) imon = imon + 1 IF ( imon > 12 ) THEN imon = 1 iyea = iyea + 1 CALL calc_month_len( iyea, imonth_len ) ! update month lengths ENDIF END DO !----------------------------------------------------------------------- ! Convert it into YYYYMMDD.HHMMSS format !----------------------------------------------------------------------- ddobsfin = iyea * 10000_dp + imon * 100_dp + iday & & + ihou * 0.01_dp + imin * 0.0001_dp END SUBROUTINE fin_date SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, jtype, ctypein, & & cfilestype, ifiles, cobstypes, cfiles ) INTEGER, INTENT(IN) :: ntypes ! Total number of obs types INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type INTEGER, INTENT(IN) :: jtype ! Index of the current type of obs INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & & ifiles ! Out appended number of files for this type CHARACTER(len=6), INTENT(IN) :: ctypein CHARACTER(len=128), DIMENSION(jpmaxnfiles), INTENT(IN) :: & & cfilestype ! In list of files for this obs type CHARACTER(len=6), DIMENSION(ntypes), INTENT(INOUT) :: & & cobstypes ! Out appended list of obs types CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(INOUT) :: & & cfiles ! Out appended list of files for all types !Local variables INTEGER :: jfile cfiles(jtype,:) = cfilestype(:) cobstypes(jtype) = ctypein ifiles(jtype) = 0 DO jfile = 1, jpmaxnfiles IF ( trim(cfiles(jtype,jfile)) /= '' ) & ifiles(jtype) = ifiles(jtype) + 1 END DO IF ( ifiles(jtype) == 0 ) THEN CALL ctl_stop( 'Logical for observation type '//TRIM(ctypein)// & & ' set to true but no files available to read' ) ENDIF IF(lwp) THEN WRITE(numout,*) ' '//cobstypes(jtype)//' input observation file names:' DO jfile = 1, ifiles(jtype) WRITE(numout,*) ' '//TRIM(cfiles(jtype,jfile)) END DO ENDIF END SUBROUTINE obs_settypefiles SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein, & & n2dint_default, n2dint_type, & & ravglamscl_type, ravgphiscl_type, & & lfp_indegs_type, lavnight_type, & & n2dint, ravglamscl, ravgphiscl, & & lfpindegs, lavnight ) INTEGER, INTENT(IN) :: ntypes ! Total number of obs types INTEGER, INTENT(IN) :: jtype ! Index of the current type of obs INTEGER, INTENT(IN) :: n2dint_default ! Default option for interpolation type INTEGER, INTENT(IN) :: n2dint_type ! Option for interpolation type REAL(wp), INTENT(IN) :: & & ravglamscl_type, & !E/W diameter of obs footprint for this type & ravgphiscl_type !N/S diameter of obs footprint for this type LOGICAL, INTENT(IN) :: lfp_indegs_type !T=> footprint in degrees, F=> in metres LOGICAL, INTENT(IN) :: lavnight_type !T=> obs represent night time average CHARACTER(len=6), INTENT(IN) :: ctypein INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & & n2dint REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: & & ravglamscl, ravgphiscl LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: & & lfpindegs, lavnight lavnight(jtype) = lavnight_type IF ( (n2dint_type >= 1) .AND. (n2dint_type <= 6) ) THEN n2dint(jtype) = n2dint_type ELSE n2dint(jtype) = n2dint_default ENDIF ! For averaging observation footprints set options for size of footprint IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN IF ( ravglamscl_type > 0._wp ) THEN ravglamscl(jtype) = ravglamscl_type ELSE CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 'scale (ravglamscl) for observation type '//TRIM(ctypein) ) ENDIF IF ( ravgphiscl_type > 0._wp ) THEN ravgphiscl(jtype) = ravgphiscl_type ELSE CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 'scale (ravgphiscl) for observation type '//TRIM(ctypein) ) ENDIF lfpindegs(jtype) = lfp_indegs_type ENDIF ! Write out info IF(lwp) THEN IF ( n2dint(jtype) <= 4 ) THEN WRITE(numout,*) ' '//TRIM(ctypein)// & & ' model counterparts will be interpolated horizontally' ELSE IF ( n2dint(jtype) <= 6 ) THEN WRITE(numout,*) ' '//TRIM(ctypein)// & & ' model counterparts will be averaged horizontally' WRITE(numout,*) ' '//' with E/W scale: ',ravglamscl(jtype) WRITE(numout,*) ' '//' with N/S scale: ',ravgphiscl(jtype) IF ( lfpindegs(jtype) ) THEN WRITE(numout,*) ' '//' (in degrees)' ELSE WRITE(numout,*) ' '//' (in metres)' ENDIF ENDIF ENDIF END SUBROUTINE obs_setinterpopts END MODULE diaobs