- Timestamp:
- 2022-03-21T13:12:29+01:00 (2 years ago)
- Location:
- branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 1 added
- 14 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
r11203 r15764 21 21 !! dyn_asm_inc : Apply the dynamic (u and v) increments 22 22 !! ssh_asm_inc : Apply the SSH increment 23 !! seaice_asm_inc : Apply the seaice increment 23 !! seaice_asm_inc : Apply the sea ice concentration increment 24 !! sit_asm_inc : Apply the sea ice thickness increment 24 25 !!---------------------------------------------------------------------- 25 26 USE wrk_nemo ! Memory Allocation … … 49 50 PUBLIC dyn_asm_inc !: Apply the dynamic (u and v) increments 50 51 PUBLIC ssh_asm_inc !: Apply the SSH increment 51 PUBLIC seaice_asm_inc !: Apply the seaice increment 52 PUBLIC seaice_asm_inc !: Apply the seaice concentration increment 53 PUBLIC sit_asm_inc !: Apply the seaice thickness increment 54 PUBLIC bgc_asm_inc !: Apply the biogeochemistry increments 52 55 53 56 #if defined key_asminc … … 57 60 #endif 58 61 LOGICAL, PUBLIC :: ln_bkgwri = .FALSE. !: No output of the background state fields 62 LOGICAL, PUBLIC :: ln_avgbkg = .FALSE. !: No output of the mean background state fields 59 63 LOGICAL, PUBLIC :: ln_asmiau = .FALSE. !: No applying forcing with an assimilation increment 60 64 LOGICAL, PUBLIC :: ln_asmdin = .FALSE. !: No direct initialization … … 62 66 LOGICAL, PUBLIC :: ln_dyninc = .FALSE. !: No dynamics (u and v) assimilation increments 63 67 LOGICAL, PUBLIC :: ln_sshinc = .FALSE. !: No sea surface height assimilation increment 64 LOGICAL, PUBLIC :: ln_seaiceinc !: No sea ice concentration increment 68 LOGICAL, PUBLIC :: ln_seaiceinc = .FALSE. !: No sea ice concentration increment 69 LOGICAL, PUBLIC :: ln_sitinc = .FALSE. !: No sea ice thickness increment 70 LOGICAL, PUBLIC :: lk_bgcinc = .FALSE. !: No biogeochemistry increments 65 71 LOGICAL, PUBLIC :: ln_salfix = .FALSE. !: Apply minimum salinity check 66 72 LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing … … 85 91 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ssh_bkg, ssh_bkginc ! Background sea surface height and its increment 86 92 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: seaice_bkginc ! Increment to the background sea ice conc 93 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: sit_bkginc ! Increment to the background sea ice thickness 87 94 88 95 !! * Substitutions … … 127 134 REAL(wp), POINTER, DIMENSION(:,:) :: hdiv ! 2D workspace 128 135 !! 129 NAMELIST/nam_asminc/ ln_bkgwri, 136 NAMELIST/nam_asminc/ ln_bkgwri, ln_avgbkg, ln_balwri, & 130 137 & ln_trainc, ln_dyninc, ln_sshinc, & 131 138 & ln_asmdin, ln_asmiau, & … … 138 145 !----------------------------------------------------------------------- 139 146 ln_seaiceinc = .FALSE. 147 ln_sitinc = .FALSE. 140 148 ln_temnofreeze = .FALSE. 141 149 … … 154 162 WRITE(numout,*) 'asm_inc_init : Assimilation increment initialization :' 155 163 WRITE(numout,*) '~~~~~~~~~~~~' 156 WRITE(numout,*) ' Namelist nam asm: set assimilation increment parameters'164 WRITE(numout,*) ' Namelist nam_asminc : set assimilation increment parameters' 157 165 WRITE(numout,*) ' Logical switch for writing out background state ln_bkgwri = ', ln_bkgwri 166 WRITE(numout,*) ' Logical switch for writing mean background state ln_avgbkg = ', ln_avgbkg 167 WRITE(numout,*) ' Logical switch for writing out balancing increments ln_balwri = ', ln_balwri 158 168 WRITE(numout,*) ' Logical switch for applying tracer increments ln_trainc = ', ln_trainc 159 169 WRITE(numout,*) ' Logical switch for applying velocity increments ln_dyninc = ', ln_dyninc 160 170 WRITE(numout,*) ' Logical switch for applying SSH increments ln_sshinc = ', ln_sshinc 171 WRITE(numout,*) ' Logical switch for applying SIC increments ln_seaiceinc = ', ln_seaiceinc 172 WRITE(numout,*) ' Logical switch for applying SIT increments ln_sitinc = ', ln_sitinc 161 173 WRITE(numout,*) ' Logical switch for Direct Initialization (DI) ln_asmdin = ', ln_asmdin 162 WRITE(numout,*) ' Logical switch for applying sea ice increments ln_seaiceinc = ', ln_seaiceinc163 174 WRITE(numout,*) ' Logical switch for Incremental Analysis Updating (IAU) ln_asmiau = ', ln_asmiau 164 175 WRITE(numout,*) ' Timestep of background in [0,nitend-nit000-1] nitbkg = ', nitbkg … … 216 227 217 228 IF ( ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) & 218 .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) .OR. ( ln_seaiceinc) )) & 219 & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc is set to .true.', & 229 & .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ).OR. & 230 & ( ln_sitinc ).OR.( lk_bgcinc ) )) & 231 & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', & 232 & ' ln_sitinc and ln_(bgc-variable)inc is set to .true.', & 220 233 & ' but ln_asmdin and ln_asmiau are both set to .false. :', & 221 234 & ' Inconsistent options') … … 226 239 227 240 IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) & 228 & ) & 229 & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc are set to .false. :', & 241 & .AND.( .NOT. ln_sitinc ).AND.( .NOT. lk_bgcinc ) ) & 242 & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', & 243 & ' ln_sitinc and ln_(bgc-variable)inc are set to .false. :', & 230 244 & ' The assimilation increments are not applied') 231 245 … … 325 339 !-------------------------------------------------------------------- 326 340 327 ALLOCATE( t_bkginc(jpi,jpj,jpk) ) 328 ALLOCATE( s_bkginc(jpi,jpj,jpk) ) 329 ALLOCATE( u_bkginc(jpi,jpj,jpk) ) 330 ALLOCATE( v_bkginc(jpi,jpj,jpk) ) 331 ALLOCATE( ssh_bkginc(jpi,jpj) ) 332 ALLOCATE( seaice_bkginc(jpi,jpj)) 341 IF ( ln_trainc ) THEN 342 ALLOCATE( t_bkginc(jpi,jpj,jpk) ) 343 ALLOCATE( s_bkginc(jpi,jpj,jpk) ) 344 t_bkginc(:,:,:) = 0.0 345 s_bkginc(:,:,:) = 0.0 346 ENDIF 347 IF ( ln_dyninc ) THEN 348 ALLOCATE( u_bkginc(jpi,jpj,jpk) ) 349 ALLOCATE( v_bkginc(jpi,jpj,jpk) ) 350 u_bkginc(:,:,:) = 0.0 351 v_bkginc(:,:,:) = 0.0 352 ENDIF 353 IF ( ln_sshinc ) THEN 354 ALLOCATE( ssh_bkginc(jpi,jpj) ) 355 ssh_bkginc(:,:) = 0.0 356 ENDIF 357 IF ( ln_seaiceinc ) THEN 358 ALLOCATE( seaice_bkginc(jpi,jpj)) 359 seaice_bkginc(:,:) = 0.0 360 ENDIF 361 IF ( ln_sitinc ) THEN 362 ALLOCATE( sit_bkginc(jpi,jpj)) 363 sit_bkginc(:,:) = 0.0 364 ENDIF 333 365 #if defined key_asminc 334 366 ALLOCATE( ssh_iau(jpi,jpj) ) 335 #endif336 t_bkginc(:,:,:) = 0.0337 s_bkginc(:,:,:) = 0.0338 u_bkginc(:,:,:) = 0.0339 v_bkginc(:,:,:) = 0.0340 ssh_bkginc(:,:) = 0.0341 seaice_bkginc(:,:) = 0.0342 #if defined key_asminc343 367 ssh_iau(:,:) = 0.0 344 368 #endif 345 IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) ) THEN 369 IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) & 370 & .OR.( ln_sitinc ).OR.( lk_bgcinc ) ) THEN 346 371 347 372 !-------------------------------------------------------------------- … … 406 431 ! to allow for differences in masks 407 432 WHERE( ABS( ssh_bkginc(:,:) ) > 1.0e+10 ) ssh_bkginc(:,:) = 0.0 433 ENDIF 434 435 IF ( ln_sitinc ) THEN 436 CALL iom_get( inum, jpdom_autoglo, 'bckinsit', sit_bkginc, 1 ) 437 ! Apply the masks 438 sit_bkginc(:,:) = sit_bkginc(:,:) * tmask(:,:,1) 439 ! Set missing increments to 0.0 rather than 1e+20 440 ! to allow for differences in masks 441 WHERE( ABS( sit_bkginc(:,:) ) > 1.0e+10 ) sit_bkginc(:,:) = 0.0 408 442 ENDIF 409 443 … … 766 800 ! Perhaps the following call should be in step 767 801 IF ( ln_seaiceinc ) CALL seaice_asm_inc ( kt ) ! apply sea ice concentration increment 802 IF ( ln_sitinc ) CALL sit_asm_inc ( kt ) ! apply sea ice thickness increment 768 803 ! 769 804 END SUBROUTINE tra_asm_inc … … 923 958 END SUBROUTINE ssh_asm_inc 924 959 960 SUBROUTINE sit_asm_inc( kt, kindic ) 961 !!---------------------------------------------------------------------- 962 !! *** ROUTINE sit_asm_inc *** 963 !! 964 !! ** Purpose : Apply the sea ice thickness assimilation increment. 965 !! 966 !! ** Method : Direct initialization or Incremental Analysis Updating. 967 !! 968 !! ** Action : 969 !! 970 !!---------------------------------------------------------------------- 971 IMPLICIT NONE 972 ! 973 INTEGER, INTENT(in) :: kt ! Current time step 974 INTEGER, INTENT(in), OPTIONAL :: kindic ! flag for disabling the deallocation 975 ! 976 INTEGER :: it 977 REAL(wp) :: zincwgt ! IAU weight for current time step 978 !!---------------------------------------------------------------------- 979 980 IF ( ln_asmiau ) THEN 981 982 !-------------------------------------------------------------------- 983 ! Incremental Analysis Updating 984 !-------------------------------------------------------------------- 985 986 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 987 988 it = kt - nit000 + 1 989 zincwgt = wgtiau(it) ! IAU weight for the current time step 990 ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments) 991 ! EF: Actually CICE is expecting a tendency so is divided by rdt below 992 993 IF(lwp) THEN 994 WRITE(numout,*) 995 WRITE(numout,*) 'sit_asm_inc : sea ice thick IAU at time step = ', & 996 & kt,' with IAU weight = ', wgtiau(it) 997 WRITE(numout,*) '~~~~~~~~~~~~' 998 ENDIF 999 1000 #if defined key_cice && defined key_asminc 1001 ! Sea-ice thickness : CICE case. Pass ice thickness increment tendency into CICE 1002 ndsit_da(:,:) = sit_bkginc(:,:) * zincwgt / rdt 1003 #endif 1004 1005 IF ( kt == nitiaufin_r ) THEN 1006 DEALLOCATE( sit_bkginc ) 1007 ENDIF 1008 1009 ELSE 1010 1011 #if defined key_cice && defined key_asminc 1012 ! Sea-ice thickness : CICE case. Zero ice increment tendency into CICE 1013 ndsit_da(:,:) = 0.0_wp 1014 #endif 1015 1016 ENDIF 1017 1018 ELSEIF ( ln_asmdin ) THEN 1019 1020 !-------------------------------------------------------------------- 1021 ! Direct Initialization 1022 !-------------------------------------------------------------------- 1023 1024 IF ( kt == nitdin_r ) THEN 1025 1026 neuler = 0 ! Force Euler forward step 1027 1028 #if defined key_cice && defined key_asminc 1029 ! Sea-ice thickness : CICE case. Pass ice thickness increment tendency into CICE 1030 ndsit_da(:,:) = sit_bkginc(:,:) / rdt 1031 #endif 1032 IF ( .NOT. PRESENT(kindic) ) THEN 1033 DEALLOCATE( sit_bkginc ) 1034 END IF 1035 1036 ELSE 1037 1038 #if defined key_cice && defined key_asminc 1039 ! Sea-ice thicnkness : CICE case. Zero ice thickness increment tendency into CICE 1040 ndsit_da(:,:) = 0.0_wp 1041 #endif 1042 1043 ENDIF 1044 1045 ENDIF 1046 1047 END SUBROUTINE sit_asm_inc 925 1048 926 1049 SUBROUTINE seaice_asm_inc( kt, kindic ) … … 928 1051 !! *** ROUTINE seaice_asm_inc *** 929 1052 !! 930 !! ** Purpose : Apply the sea ice assimilation increment.1053 !! ** Purpose : Apply the sea ice concentration assimilation increment. 931 1054 !! 932 1055 !! ** Method : Direct initialization or Incremental Analysis Updating. -
branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r11203 r15764 45 45 USE agrif_opa_interp ! agrif 46 46 #endif 47 #if defined key_asminc48 USE asminc ! Assimilation increment49 #endif50 47 51 48 IMPLICIT NONE … … 459 456 & + fwfisf(:,:) + fwfisf_b(:,:) ) 460 457 ENDIF 461 #if defined key_asminc 462 ! ! Include the IAU weighted SSH increment 463 IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 464 zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:) 465 ENDIF 466 #endif 467 ! !* Fill boundary data arrays for AGRIF 468 ! ! ------------------------------------ 458 ! !* Fill boundary data arrays with AGRIF 459 ! ! ------------------------------------- 469 460 #if defined key_agrif 470 461 IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt ) -
branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
r11239 r15764 22 22 USE obs_read_surf ! Reading and allocation of surface obs 23 23 USE obs_readmdt ! Reading and allocation of MDT for SLA. 24 USE obs_readsnowdepth ! Get model snow depth for conversion of freeboard to ice thickness 24 25 USE obs_prep ! Preparation of obs. (grid search etc). 25 26 USE obs_oper ! Observation operators … … 33 34 USE mpp_map ! MPP mapping 34 35 USE lib_mpp ! For ctl_warn/stop 36 USE tradmp ! For climatological temperature & salinity 35 37 36 38 IMPLICIT NONE … … 44 46 45 47 !! * Module variables 46 LOGICAL, PUBLIC :: & 47 & lk_diaobs = .TRUE. !: Include this for backwards compatibility at NEMO 3.6. 48 LOGICAL :: ln_diaobs !: Logical switch for the obs operator 49 LOGICAL :: ln_sstnight !: Logical switch for night mean SST obs 50 LOGICAL :: ln_default_fp_indegs !: T=> Default obs footprint size specified in degrees, F=> in metres 51 LOGICAL :: ln_sla_fp_indegs !: T=> SLA obs footprint size specified in degrees, F=> in metres 52 LOGICAL :: ln_sst_fp_indegs !: T=> SST obs footprint size specified in degrees, F=> in metres 53 LOGICAL :: ln_sss_fp_indegs !: T=> SSS obs footprint size specified in degrees, F=> in metres 54 LOGICAL :: ln_sic_fp_indegs !: T=> sea-ice obs footprint size specified in degrees, F=> in metres 55 56 REAL(wp) :: rn_default_avglamscl !: Default E/W diameter of observation footprint 57 REAL(wp) :: rn_default_avgphiscl !: Default N/S diameter of observation footprint 58 REAL(wp) :: rn_sla_avglamscl !: E/W diameter of SLA observation footprint 59 REAL(wp) :: rn_sla_avgphiscl !: N/S diameter of SLA observation footprint 60 REAL(wp) :: rn_sst_avglamscl !: E/W diameter of SST observation footprint 61 REAL(wp) :: rn_sst_avgphiscl !: N/S diameter of SST observation footprint 62 REAL(wp) :: rn_sss_avglamscl !: E/W diameter of SSS observation footprint 63 REAL(wp) :: rn_sss_avgphiscl !: N/S diameter of SSS observation footprint 64 REAL(wp) :: rn_sic_avglamscl !: E/W diameter of sea-ice observation footprint 65 REAL(wp) :: rn_sic_avgphiscl !: N/S diameter of sea-ice observation footprint 66 67 INTEGER :: nn_1dint !: Vertical interpolation method 68 INTEGER :: nn_2dint_default !: Default horizontal interpolation method 69 INTEGER :: nn_2dint_sla !: SLA horizontal interpolation method (-1 = default) 70 INTEGER :: nn_2dint_sst !: SST horizontal interpolation method (-1 = default) 71 INTEGER :: nn_2dint_sss !: SSS horizontal interpolation method (-1 = default) 72 INTEGER :: nn_2dint_sic !: Seaice horizontal interpolation method (-1 = default) 73 48 LOGICAL, PUBLIC :: & 49 & lk_diaobs = .TRUE. !: Include this for backwards compatibility at NEMO 3.6. 50 LOGICAL :: ln_diaobs !: Logical switch for the obs operator 51 LOGICAL :: ln_sstnight !: Logical switch for night mean SST obs 52 LOGICAL :: ln_default_fp_indegs !: T=> Default obs footprint size specified in degrees, F=> in metres 53 LOGICAL :: ln_sla_fp_indegs !: T=> SLA obs footprint size specified in degrees, F=> in metres 54 LOGICAL :: ln_sst_fp_indegs !: T=> SST obs footprint size specified in degrees, F=> in metres 55 LOGICAL :: ln_sss_fp_indegs !: T=> SSS obs footprint size specified in degrees, F=> in metres 56 LOGICAL :: ln_ssv_fp_indegs !: T=> SSV obs footprint size specified in degrees, F=> in metres 57 LOGICAL :: ln_sic_fp_indegs !: T=> SIC obs footprint size specified in degrees, F=> in metres 58 LOGICAL :: ln_sit_fp_indegs !: T=> SIT obs footprint size specified in degrees, F=> in metres 59 LOGICAL :: ln_fbd_fp_indegs !: T=> SIT obs footprint size specified in degrees, F=> in metres 60 LOGICAL :: ln_output_clim !: Logical switch for interpolating and writing T/S climatology 61 LOGICAL :: ln_time_mean_sla_bkg !: Logical switch for applying time mean of SLA background to remove tidal signal 62 63 REAL(wp) :: rn_default_avglamscl !: Default E/W diameter of observation footprint 64 REAL(wp) :: rn_default_avgphiscl !: Default N/S diameter of observation footprint 65 REAL(wp) :: rn_sla_avglamscl !: E/W diameter of SLA observation footprint 66 REAL(wp) :: rn_sla_avgphiscl !: N/S diameter of SLA observation footprint 67 REAL(wp) :: rn_sst_avglamscl !: E/W diameter of SST observation footprint 68 REAL(wp) :: rn_sst_avgphiscl !: N/S diameter of SST observation footprint 69 REAL(wp) :: rn_sss_avglamscl !: E/W diameter of SSS observation footprint 70 REAL(wp) :: rn_sss_avgphiscl !: N/S diameter of SSS observation footprint 71 REAL(wp) :: rn_ssv_avglamscl !: E/W diameter of SSV observation footprint 72 REAL(wp) :: rn_ssv_avgphiscl !: N/S diameter of SSV observation footprint 73 REAL(wp) :: rn_sic_avglamscl !: E/W diameter of SIC observation footprint 74 REAL(wp) :: rn_sic_avgphiscl !: N/S diameter of SIC observation footprint 75 REAL(wp) :: rn_sit_avglamscl !: E/W diameter of SIT observation footprint 76 REAL(wp) :: rn_sit_avgphiscl !: N/S diameter of SIT observation footprint 77 REAL(wp) :: rn_fbd_avglamscl !: E/W diameter of FBD observation footprint 78 REAL(wp) :: rn_fbd_avgphiscl !: N/S diameter of FBD observation footprint 79 REAL(wp), PUBLIC :: & 80 & MeanPeriodHours = 24. + (5./6.) !: Meaning period for surface data. 81 82 83 INTEGER :: nn_1dint !: Vertical interpolation method 84 INTEGER :: nn_2dint_default !: Default horizontal interpolation method 85 INTEGER :: nn_2dint_sla !: SLA horizontal interpolation method (-1 = default) 86 INTEGER :: nn_2dint_sst !: SST horizontal interpolation method (-1 = default) 87 INTEGER :: nn_2dint_sss !: SSS horizontal interpolation method (-1 = default) 88 INTEGER :: nn_2dint_ssv !: SSV horizontal interpolation method (-1 = default) 89 INTEGER :: nn_2dint_sic !: SIC horizontal interpolation method (-1 = default) 90 INTEGER :: nn_2dint_sit !: SIT horizontal interpolation method (-1 = default) 91 INTEGER :: nn_2dint_fbd !: FBD horizontal interpolation method (-1 = default) 92 74 93 INTEGER, DIMENSION(imaxavtypes) :: & 75 94 & nn_profdavtypes !: Profile data types representing a daily average … … 130 149 !! ! 15-02 (M. Martin) Simplification of namelist and code 131 150 !!---------------------------------------------------------------------- 151 #if defined key_cice 152 USE sbc_oce, ONLY : & ! CICE variables 153 & thick_s ! snow depth for freeboard conversion 154 #endif 132 155 133 156 IMPLICIT NONE 134 157 135 158 !! * Local declarations 136 INTEGER, PARAMETER :: & 137 & jpmaxnfiles = 1000 ! Maximum number of files for each obs type 138 INTEGER, DIMENSION(:), ALLOCATABLE :: & 139 & ifilesprof, & ! Number of profile files 140 & ifilessurf ! Number of surface files 141 INTEGER :: ios ! Local integer output status for namelist read 142 INTEGER :: jtype ! Counter for obs types 143 INTEGER :: jvar ! Counter for variables 144 INTEGER :: jfile ! Counter for files 145 INTEGER :: jnumsstbias ! Number of SST bias files to read and apply 146 INTEGER :: n2dint_type ! Local version of nn_2dint* 147 148 CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: & 149 & cn_profbfiles, & ! T/S profile input filenames 150 & cn_sstfbfiles, & ! Sea surface temperature input filenames 151 & cn_slafbfiles, & ! Sea level anomaly input filenames 152 & cn_sicfbfiles, & ! Seaice concentration input filenames 153 & cn_velfbfiles, & ! Velocity profile input filenames 154 & cn_sssfbfiles, & ! Sea surface salinity input filenames 155 & cn_slchltotfbfiles, & ! Surface total log10(chlorophyll) input filenames 156 & cn_slchldiafbfiles, & ! Surface diatom log10(chlorophyll) input filenames 157 & cn_slchlnonfbfiles, & ! Surface non-diatom log10(chlorophyll) input filenames 158 & cn_slchldinfbfiles, & ! Surface dinoflagellate log10(chlorophyll) input filenames 159 & cn_slchlmicfbfiles, & ! Surface microphytoplankton log10(chlorophyll) input filenames 160 & cn_slchlnanfbfiles, & ! Surface nanophytoplankton log10(chlorophyll) input filenames 161 & cn_slchlpicfbfiles, & ! Surface picophytoplankton log10(chlorophyll) input filenames 162 & cn_schltotfbfiles, & ! Surface total chlorophyll input filenames 163 & cn_slphytotfbfiles, & ! Surface total log10(phytoplankton carbon) input filenames 164 & cn_slphydiafbfiles, & ! Surface diatom log10(phytoplankton carbon) input filenames 165 & cn_slphynonfbfiles, & ! Surface non-diatom log10(phytoplankton carbon) input filenames 166 & cn_sspmfbfiles, & ! Surface suspended particulate matter input filenames 167 & cn_sfco2fbfiles, & ! Surface fugacity of carbon dioxide input filenames 168 & cn_spco2fbfiles, & ! Surface partial pressure of carbon dioxide input filenames 169 & cn_plchltotfbfiles, & ! Profile total log10(chlorophyll) input filenames 170 & cn_pchltotfbfiles, & ! Profile total chlorophyll input filenames 171 & cn_pno3fbfiles, & ! Profile nitrate input filenames 172 & cn_psi4fbfiles, & ! Profile silicate input filenames 173 & cn_ppo4fbfiles, & ! Profile phosphate input filenames 174 & cn_pdicfbfiles, & ! Profile dissolved inorganic carbon input filenames 175 & cn_palkfbfiles, & ! Profile alkalinity input filenames 176 & cn_pphfbfiles, & ! Profile pH input filenames 177 & cn_po2fbfiles, & ! Profile dissolved oxygen input filenames 178 & cn_sstbiasfiles ! SST bias input filenames 179 180 CHARACTER(LEN=128) :: & 181 & cn_altbiasfile ! Altimeter bias input filename 182 183 184 LOGICAL :: ln_t3d ! Logical switch for temperature profiles 185 LOGICAL :: ln_s3d ! Logical switch for salinity profiles 186 LOGICAL :: ln_sla ! Logical switch for sea level anomalies 187 LOGICAL :: ln_sst ! Logical switch for sea surface temperature 188 LOGICAL :: ln_sic ! Logical switch for sea ice concentration 189 LOGICAL :: ln_sss ! Logical switch for sea surface salinity obs 190 LOGICAL :: ln_vel3d ! Logical switch for velocity (u,v) obs 191 LOGICAL :: ln_slchltot ! Logical switch for surface total log10(chlorophyll) obs 192 LOGICAL :: ln_slchldia ! Logical switch for surface diatom log10(chlorophyll) obs 193 LOGICAL :: ln_slchlnon ! Logical switch for surface non-diatom log10(chlorophyll) obs 194 LOGICAL :: ln_slchldin ! Logical switch for surface dinoflagellate log10(chlorophyll) obs 195 LOGICAL :: ln_slchlmic ! Logical switch for surface microphytoplankton log10(chlorophyll) obs 196 LOGICAL :: ln_slchlnan ! Logical switch for surface nanophytoplankton log10(chlorophyll) obs 197 LOGICAL :: ln_slchlpic ! Logical switch for surface picophytoplankton log10(chlorophyll) obs 198 LOGICAL :: ln_schltot ! Logical switch for surface total chlorophyll obs 199 LOGICAL :: ln_slphytot ! Logical switch for surface total log10(phytoplankton carbon) obs 200 LOGICAL :: ln_slphydia ! Logical switch for surface diatom log10(phytoplankton carbon) obs 201 LOGICAL :: ln_slphynon ! Logical switch for surface non-diatom log10(phytoplankton carbon) obs 202 LOGICAL :: ln_sspm ! Logical switch for surface suspended particulate matter obs 203 LOGICAL :: ln_sfco2 ! Logical switch for surface fugacity of carbon dioxide obs 204 LOGICAL :: ln_spco2 ! Logical switch for surface partial pressure of carbon dioxide obs 205 LOGICAL :: ln_plchltot ! Logical switch for profile total log10(chlorophyll) obs 206 LOGICAL :: ln_pchltot ! Logical switch for profile total chlorophyll obs 207 LOGICAL :: ln_pno3 ! Logical switch for profile nitrate obs 208 LOGICAL :: ln_psi4 ! Logical switch for profile silicate obs 209 LOGICAL :: ln_ppo4 ! Logical switch for profile phosphate obs 210 LOGICAL :: ln_pdic ! Logical switch for profile dissolved inorganic carbon obs 211 LOGICAL :: ln_palk ! Logical switch for profile alkalinity obs 212 LOGICAL :: ln_pph ! Logical switch for profile pH obs 213 LOGICAL :: ln_po2 ! Logical switch for profile dissolved oxygen obs 214 LOGICAL :: ln_nea ! Logical switch to remove obs near land 215 LOGICAL :: ln_altbias ! Logical switch for altimeter bias 216 LOGICAL :: ln_sstbias ! Logical switch for bias correction of SST 217 LOGICAL :: ln_ignmis ! Logical switch for ignoring missing files 218 LOGICAL :: ln_s_at_t ! Logical switch to compute model S at T obs 219 LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary 220 221 REAL(dp) :: rn_dobsini ! Obs window start date YYYYMMDD.HHMMSS 222 REAL(dp) :: rn_dobsend ! Obs window end date YYYYMMDD.HHMMSS 223 224 REAL(wp) :: ztype_avglamscl ! Local version of rn_*_avglamscl 225 REAL(wp) :: ztype_avgphiscl ! Local version of rn_*_avgphiscl 226 227 CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: & 228 & clproffiles, & ! Profile filenames 229 & clsurffiles ! Surface filenames 230 231 LOGICAL, DIMENSION(:), ALLOCATABLE :: llvar ! Logical for profile variable read 232 LOGICAL :: ltype_fp_indegs ! Local version of ln_*_fp_indegs 233 LOGICAL :: ltype_night ! Local version of ln_sstnight (false for other variables) 234 235 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 236 & zglam ! Model longitudes for profile variables 237 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 238 & zgphi ! Model latitudes for profile variables 239 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 240 & zmask ! Model land/sea mask associated with variables 241 242 243 NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla, & 244 & ln_sst, ln_sic, ln_sss, ln_vel3d, & 245 & ln_slchltot, ln_slchldia, ln_slchlnon, & 246 & ln_slchldin, ln_slchlmic, ln_slchlnan, & 247 & ln_slchlpic, ln_schltot, & 248 & ln_slphytot, ln_slphydia, ln_slphynon, & 249 & ln_sspm, ln_sfco2, ln_spco2, & 250 & ln_plchltot, ln_pchltot, ln_pno3, & 251 & ln_psi4, ln_ppo4, ln_pdic, & 252 & ln_palk, ln_pph, ln_po2, & 253 & ln_altbias, ln_sstbias, ln_nea, & 254 & ln_grid_global, ln_grid_search_lookup, & 255 & ln_ignmis, ln_s_at_t, ln_bound_reject, & 256 & ln_sstnight, ln_default_fp_indegs, & 257 & ln_sla_fp_indegs, ln_sst_fp_indegs, & 258 & ln_sss_fp_indegs, ln_sic_fp_indegs, & 259 & cn_profbfiles, cn_slafbfiles, & 260 & cn_sstfbfiles, cn_sicfbfiles, & 261 & cn_velfbfiles, cn_sssfbfiles, & 262 & cn_slchltotfbfiles, cn_slchldiafbfiles, & 263 & cn_slchlnonfbfiles, cn_slchldinfbfiles, & 264 & cn_slchlmicfbfiles, cn_slchlnanfbfiles, & 265 & cn_slchlpicfbfiles, cn_schltotfbfiles, & 266 & cn_slphytotfbfiles, cn_slphydiafbfiles, & 267 & cn_slphynonfbfiles, cn_sspmfbfiles, & 268 & cn_sfco2fbfiles, cn_spco2fbfiles, & 269 & cn_plchltotfbfiles, cn_pchltotfbfiles, & 270 & cn_pno3fbfiles, cn_psi4fbfiles, cn_ppo4fbfiles, & 271 & cn_pdicfbfiles, cn_palkfbfiles, cn_pphfbfiles, & 272 & cn_po2fbfiles, & 273 & cn_sstbiasfiles, cn_altbiasfile, & 274 & cn_gridsearchfile, rn_gridsearchres, & 275 & rn_dobsini, rn_dobsend, & 276 & rn_default_avglamscl, rn_default_avgphiscl, & 277 & rn_sla_avglamscl, rn_sla_avgphiscl, & 278 & rn_sst_avglamscl, rn_sst_avgphiscl, & 279 & rn_sss_avglamscl, rn_sss_avgphiscl, & 280 & rn_sic_avglamscl, rn_sic_avgphiscl, & 281 & nn_1dint, nn_2dint_default, & 282 & nn_2dint_sla, nn_2dint_sst, & 283 & nn_2dint_sss, nn_2dint_sic, & 284 & nn_msshc, rn_mdtcorr, rn_mdtcutoff, & 159 INTEGER, PARAMETER :: & 160 & jpmaxnfiles = 1000 ! Maximum number of files for each obs type 161 INTEGER, DIMENSION(:), ALLOCATABLE :: & 162 & ifilesprof, & ! Number of profile files 163 & ifilessurf ! Number of surface files 164 INTEGER :: ios ! Local integer output status for namelist read 165 INTEGER :: jtype ! Counter for obs types 166 INTEGER :: jvar ! Counter for variables 167 INTEGER :: jfile ! Counter for files 168 INTEGER :: jnumsstbias ! Number of SST bias files to read and apply 169 INTEGER :: n2dint_type ! Local version of nn_2dint* 170 171 CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: & 172 & cn_profbfiles, & ! T/S profile input filenames 173 & cn_sstfbfiles, & ! Sea surface temperature input filenames 174 & cn_slafbfiles, & ! Sea level anomaly input filenames 175 & cn_sicfbfiles, & ! Seaice concentration input filenames 176 & cn_sitfbfiles, & ! Seaice thickness input filenames 177 & cn_fbdfbfiles, & ! Seaice freeboard input filenames 178 & cn_velfbfiles, & ! Velocity profile input filenames 179 & cn_sssfbfiles, & ! Sea surface salinity input filenames 180 & cn_ssvfbfiles, & ! Sea surface velocity input filenames 181 & cn_slchltotfbfiles, & ! Surface total log10(chlorophyll) input filenames 182 & cn_slchldiafbfiles, & ! Surface diatom log10(chlorophyll) input filenames 183 & cn_slchlnonfbfiles, & ! Surface non-diatom log10(chlorophyll) input filenames 184 & cn_slchldinfbfiles, & ! Surface dinoflagellate log10(chlorophyll) input filenames 185 & cn_slchlmicfbfiles, & ! Surface microphytoplankton log10(chlorophyll) input filenames 186 & cn_slchlnanfbfiles, & ! Surface nanophytoplankton log10(chlorophyll) input filenames 187 & cn_slchlpicfbfiles, & ! Surface picophytoplankton log10(chlorophyll) input filenames 188 & cn_schltotfbfiles, & ! Surface total chlorophyll input filenames 189 & cn_slphytotfbfiles, & ! Surface total log10(phytoplankton carbon) input filenames 190 & cn_slphydiafbfiles, & ! Surface diatom log10(phytoplankton carbon) input filenames 191 & cn_slphynonfbfiles, & ! Surface non-diatom log10(phytoplankton carbon) input filenames 192 & cn_sspmfbfiles, & ! Surface suspended particulate matter input filenames 193 & cn_skd490fbfiles, & ! Surface Kd490 input filenames 194 & cn_sfco2fbfiles, & ! Surface fugacity of carbon dioxide input filenames 195 & cn_spco2fbfiles, & ! Surface partial pressure of carbon dioxide input filenames 196 & cn_plchltotfbfiles, & ! Profile total log10(chlorophyll) input filenames 197 & cn_pchltotfbfiles, & ! Profile total chlorophyll input filenames 198 & cn_pno3fbfiles, & ! Profile nitrate input filenames 199 & cn_psi4fbfiles, & ! Profile silicate input filenames 200 & cn_ppo4fbfiles, & ! Profile phosphate input filenames 201 & cn_pdicfbfiles, & ! Profile dissolved inorganic carbon input filenames 202 & cn_palkfbfiles, & ! Profile alkalinity input filenames 203 & cn_pphfbfiles, & ! Profile pH input filenames 204 & cn_po2fbfiles, & ! Profile dissolved oxygen input filenames 205 & cn_sstbiasfiles ! SST bias input filenames 206 207 CHARACTER(LEN=128) :: & 208 & cn_altbiasfile ! Altimeter bias input filename 209 210 211 LOGICAL :: ln_seaicetypes = .FALSE. ! Logical switch indicating data type is sea ice 212 LOGICAL :: ln_t3d ! Logical switch for temperature profiles 213 LOGICAL :: ln_s3d ! Logical switch for salinity profiles 214 LOGICAL :: ln_sla ! Logical switch for sea level anomalies 215 LOGICAL :: ln_sst ! Logical switch for sea surface temperature 216 LOGICAL :: ln_sic ! Logical switch for sea ice concentration 217 LOGICAL :: ln_sit ! Logical switch for sea ice thickness 218 LOGICAL :: ln_fbd ! Logical switch for sea ice freeboard 219 LOGICAL :: ln_sss ! Logical switch for sea surface salinity obs 220 LOGICAL :: ln_ssv ! Logical switch for sea surface velocity obs 221 LOGICAL :: ln_vel3d ! Logical switch for velocity (u,v) obs 222 LOGICAL :: ln_slchltot ! Logical switch for surface total log10(chlorophyll) obs 223 LOGICAL :: ln_slchldia ! Logical switch for surface diatom log10(chlorophyll) obs 224 LOGICAL :: ln_slchlnon ! Logical switch for surface non-diatom log10(chlorophyll) obs 225 LOGICAL :: ln_slchldin ! Logical switch for surface dinoflagellate log10(chlorophyll) obs 226 LOGICAL :: ln_slchlmic ! Logical switch for surface microphytoplankton log10(chlorophyll) obs 227 LOGICAL :: ln_slchlnan ! Logical switch for surface nanophytoplankton log10(chlorophyll) obs 228 LOGICAL :: ln_slchlpic ! Logical switch for surface picophytoplankton log10(chlorophyll) obs 229 LOGICAL :: ln_schltot ! Logical switch for surface total chlorophyll obs 230 LOGICAL :: ln_slphytot ! Logical switch for surface total log10(phytoplankton carbon) obs 231 LOGICAL :: ln_slphydia ! Logical switch for surface diatom log10(phytoplankton carbon) obs 232 LOGICAL :: ln_slphynon ! Logical switch for surface non-diatom log10(phytoplankton carbon) obs 233 LOGICAL :: ln_sspm ! Logical switch for surface suspended particulate matter obs 234 LOGICAL :: ln_skd490 ! Logical switch for surface Kd490 235 LOGICAL :: ln_sfco2 ! Logical switch for surface fugacity of carbon dioxide obs 236 LOGICAL :: ln_spco2 ! Logical switch for surface partial pressure of carbon dioxide obs 237 LOGICAL :: ln_plchltot ! Logical switch for profile total log10(chlorophyll) obs 238 LOGICAL :: ln_pchltot ! Logical switch for profile total chlorophyll obs 239 LOGICAL :: ln_pno3 ! Logical switch for profile nitrate obs 240 LOGICAL :: ln_psi4 ! Logical switch for profile silicate obs 241 LOGICAL :: ln_ppo4 ! Logical switch for profile phosphate obs 242 LOGICAL :: ln_pdic ! Logical switch for profile dissolved inorganic carbon obs 243 LOGICAL :: ln_palk ! Logical switch for profile alkalinity obs 244 LOGICAL :: ln_pph ! Logical switch for profile pH obs 245 LOGICAL :: ln_po2 ! Logical switch for profile dissolved oxygen obs 246 LOGICAL :: ln_nea ! Logical switch to remove obs near land 247 LOGICAL :: ln_altbias ! Logical switch for altimeter bias 248 LOGICAL :: ln_sstbias ! Logical switch for bias correction of SST 249 LOGICAL :: ln_ignmis ! Logical switch for ignoring missing files 250 LOGICAL :: ln_s_at_t ! Logical switch to compute model S at T obs 251 LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary 252 253 REAL(dp) :: rn_dobsini ! Obs window start date YYYYMMDD.HHMMSS 254 REAL(dp) :: rn_dobsend ! Obs window end date YYYYMMDD.HHMMSS 255 256 REAL(wp) :: ztype_avglamscl ! Local version of rn_*_avglamscl 257 REAL(wp) :: ztype_avgphiscl ! Local version of rn_*_avgphiscl 258 259 CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: & 260 & clproffiles, & ! Profile filenames 261 & clsurffiles ! Surface filenames 262 CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars ! Expected variable names 263 264 LOGICAL, DIMENSION(:), ALLOCATABLE :: llvar ! Logical for profile variable read 265 LOGICAL :: ltype_fp_indegs ! Local version of ln_*_fp_indegs 266 LOGICAL :: ltype_night ! Local version of ln_sstnight (false for other variables) 267 LOGICAL :: ltype_clim ! Local version of ln_output_clim 268 269 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 270 & zglam ! Model longitudes for profile variables 271 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 272 & zgphi ! Model latitudes for profile variables 273 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 274 & zmask ! Model land/sea mask associated with variables 275 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 276 & zmask_surf ! Surface model land/sea mask associated with variables 277 278 279 NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla, & 280 & ln_sst, ln_sic, ln_sit, ln_fbd, & 281 & ln_sss, ln_ssv, ln_vel3d, & 282 & ln_slchltot, ln_slchldia, ln_slchlnon, & 283 & ln_slchldin, ln_slchlmic, ln_slchlnan, & 284 & ln_slchlpic, ln_schltot, & 285 & ln_slphytot, ln_slphydia, ln_slphynon, & 286 & ln_sspm, ln_sfco2, ln_spco2, & 287 & ln_skd490, & 288 & ln_plchltot, ln_pchltot, ln_pno3, & 289 & ln_psi4, ln_ppo4, ln_pdic, & 290 & ln_palk, ln_pph, ln_po2, & 291 & ln_altbias, ln_sstbias, ln_nea, & 292 & ln_grid_global, ln_grid_search_lookup, & 293 & ln_ignmis, ln_s_at_t, ln_bound_reject, & 294 & ln_sstnight, ln_output_clim, & 295 & ln_time_mean_sla_bkg, ln_default_fp_indegs, & 296 & ln_sla_fp_indegs, ln_sst_fp_indegs, & 297 & ln_sss_fp_indegs, ln_ssv_fp_indegs, & 298 & ln_sic_fp_indegs, ln_sit_fp_indegs, & 299 & ln_fbd_fp_indegs, & 300 & cn_profbfiles, cn_slafbfiles, & 301 & cn_sstfbfiles, cn_sicfbfiles, & 302 & cn_sitfbfiles, cn_fbdfbfiles, & 303 & cn_velfbfiles, cn_sssfbfiles, cn_ssvfbfiles, & 304 & cn_slchltotfbfiles, cn_slchldiafbfiles, & 305 & cn_slchlnonfbfiles, cn_slchldinfbfiles, & 306 & cn_slchlmicfbfiles, cn_slchlnanfbfiles, & 307 & cn_slchlpicfbfiles, cn_schltotfbfiles, & 308 & cn_slphytotfbfiles, cn_slphydiafbfiles, & 309 & cn_slphynonfbfiles, cn_sspmfbfiles, & 310 & cn_skd490fbfiles, & 311 & cn_sfco2fbfiles, cn_spco2fbfiles, & 312 & cn_plchltotfbfiles, cn_pchltotfbfiles, & 313 & cn_pno3fbfiles, cn_psi4fbfiles, cn_ppo4fbfiles, & 314 & cn_pdicfbfiles, cn_palkfbfiles, cn_pphfbfiles, & 315 & cn_po2fbfiles, & 316 & cn_sstbiasfiles, cn_altbiasfile, & 317 & cn_gridsearchfile, rn_gridsearchres, & 318 & rn_dobsini, rn_dobsend, & 319 & rn_default_avglamscl, rn_default_avgphiscl, & 320 & rn_sla_avglamscl, rn_sla_avgphiscl, & 321 & rn_sst_avglamscl, rn_sst_avgphiscl, & 322 & rn_sss_avglamscl, rn_sss_avgphiscl, & 323 & rn_ssv_avglamscl, rn_ssv_avgphiscl, & 324 & rn_sic_avglamscl, rn_sic_avgphiscl, & 325 & rn_sit_avglamscl, rn_sit_avgphiscl, & 326 & rn_fbd_avglamscl, rn_fbd_avgphiscl, & 327 & nn_1dint, nn_2dint_default, & 328 & nn_2dint_sla, nn_2dint_sst, & 329 & nn_2dint_sss, nn_2dint_ssv, & 330 & nn_2dint_sic, nn_2dint_sit, & 331 & nn_2dint_fbd, & 332 & nn_msshc, rn_mdtcorr, rn_mdtcutoff, & 285 333 & nn_profdavtypes 286 334 … … 289 337 !----------------------------------------------------------------------- 290 338 291 ! Some namelist arrays need initialising 292 cn_profbfiles(:) = '' 293 cn_slafbfiles(:) = '' 294 cn_sstfbfiles(:) = '' 295 cn_sicfbfiles(:) = '' 296 cn_velfbfiles(:) = '' 297 cn_sssfbfiles(:) = '' 298 cn_slchltotfbfiles(:) = '' 299 cn_slchldiafbfiles(:) = '' 300 cn_slchlnonfbfiles(:) = '' 301 cn_slchldinfbfiles(:) = '' 302 cn_slchlmicfbfiles(:) = '' 303 cn_slchlnanfbfiles(:) = '' 304 cn_slchlpicfbfiles(:) = '' 305 cn_schltotfbfiles(:) = '' 306 cn_slphytotfbfiles(:) = '' 307 cn_slphydiafbfiles(:) = '' 308 cn_slphynonfbfiles(:) = '' 309 cn_sspmfbfiles(:) = '' 310 cn_sfco2fbfiles(:) = '' 311 cn_spco2fbfiles(:) = '' 312 cn_plchltotfbfiles(:) = '' 313 cn_pchltotfbfiles(:) = '' 314 cn_pno3fbfiles(:) = '' 315 cn_psi4fbfiles(:) = '' 316 cn_ppo4fbfiles(:) = '' 317 cn_pdicfbfiles(:) = '' 318 cn_palkfbfiles(:) = '' 319 cn_pphfbfiles(:) = '' 320 cn_po2fbfiles(:) = '' 321 cn_sstbiasfiles(:) = '' 322 nn_profdavtypes(:) = -1 323 324 CALL ini_date( rn_dobsini ) 325 CALL fin_date( rn_dobsend ) 326 327 ! Read namelist namobs : control observation diagnostics 339 cn_profbfiles(:) = '' 340 cn_slafbfiles(:) = '' 341 cn_sstfbfiles(:) = '' 342 cn_sicfbfiles(:) = '' 343 cn_sitfbfiles(:) = '' 344 cn_fbdfbfiles(:) = '' 345 cn_velfbfiles(:) = '' 346 cn_sssfbfiles(:) = '' 347 cn_ssvfbfiles(:) = '' 348 cn_slchltotfbfiles(:) = '' 349 cn_slchldiafbfiles(:) = '' 350 cn_slchlnonfbfiles(:) = '' 351 cn_slchldinfbfiles(:) = '' 352 cn_slchlmicfbfiles(:) = '' 353 cn_slchlnanfbfiles(:) = '' 354 cn_slchlpicfbfiles(:) = '' 355 cn_schltotfbfiles(:) = '' 356 cn_slphytotfbfiles(:) = '' 357 cn_slphydiafbfiles(:) = '' 358 cn_slphynonfbfiles(:) = '' 359 cn_sspmfbfiles(:) = '' 360 cn_skd490fbfiles(:) = '' 361 cn_sfco2fbfiles(:) = '' 362 cn_spco2fbfiles(:) = '' 363 cn_plchltotfbfiles(:) = '' 364 cn_pchltotfbfiles(:) = '' 365 cn_pno3fbfiles(:) = '' 366 cn_psi4fbfiles(:) = '' 367 cn_ppo4fbfiles(:) = '' 368 cn_pdicfbfiles(:) = '' 369 cn_palkfbfiles(:) = '' 370 cn_pphfbfiles(:) = '' 371 cn_po2fbfiles(:) = '' 372 cn_sstbiasfiles(:) = '' 373 nn_profdavtypes(:) = -1 374 375 CALL ini_date( rn_dobsini ) 376 CALL fin_date( rn_dobsend ) 377 378 ! Read namelist namobs : control observation diagnostics 328 379 REWIND( numnam_ref ) ! Namelist namobs in reference namelist 329 380 READ ( numnam_ref, namobs, IOSTAT = ios, ERR = 901) … … 351 402 WRITE(numout,*) '~~~~~~~~~~~~' 352 403 WRITE(numout,*) ' Namelist namobs : set observation diagnostic parameters' 353 WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d 354 WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d 355 WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla 356 WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst 357 WRITE(numout,*) ' Logical switch for Sea Ice observations ln_sic = ', ln_sic 358 WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d 359 WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss 360 WRITE(numout,*) ' Logical switch for surface total logchl obs ln_slchltot = ', ln_slchltot 361 WRITE(numout,*) ' Logical switch for surface diatom logchl obs ln_slchldia = ', ln_slchldia 362 WRITE(numout,*) ' Logical switch for surface non-diatom logchl obs ln_slchlnon = ', ln_slchlnon 363 WRITE(numout,*) ' Logical switch for surface dino logchl obs ln_slchldin = ', ln_slchldin 364 WRITE(numout,*) ' Logical switch for surface micro logchl obs ln_slchlmic = ', ln_slchlmic 365 WRITE(numout,*) ' Logical switch for surface nano logchl obs ln_slchlnan = ', ln_slchlnan 366 WRITE(numout,*) ' Logical switch for surface pico logchl obs ln_slchlpic = ', ln_slchlpic 367 WRITE(numout,*) ' Logical switch for surface total chl obs ln_schltot = ', ln_schltot 368 WRITE(numout,*) ' Logical switch for surface total log(phyC) obs ln_slphytot = ', ln_slphytot 369 WRITE(numout,*) ' Logical switch for surface diatom log(phyC) obs ln_slphydia = ', ln_slphydia 370 WRITE(numout,*) ' Logical switch for surface non-diatom log(phyC) obs ln_slphynon = ', ln_slphynon 371 WRITE(numout,*) ' Logical switch for surface SPM observations ln_sspm = ', ln_sspm 372 WRITE(numout,*) ' Logical switch for surface fCO2 observations ln_sfco2 = ', ln_sfco2 373 WRITE(numout,*) ' Logical switch for surface pCO2 observations ln_spco2 = ', ln_spco2 374 WRITE(numout,*) ' Logical switch for profile total logchl obs ln_plchltot = ', ln_plchltot 375 WRITE(numout,*) ' Logical switch for profile total chl obs ln_pchltot = ', ln_pchltot 376 WRITE(numout,*) ' Logical switch for profile nitrate obs ln_pno3 = ', ln_pno3 377 WRITE(numout,*) ' Logical switch for profile silicate obs ln_psi4 = ', ln_psi4 378 WRITE(numout,*) ' Logical switch for profile phosphate obs ln_ppo4 = ', ln_ppo4 379 WRITE(numout,*) ' Logical switch for profile DIC obs ln_pdic = ', ln_pdic 380 WRITE(numout,*) ' Logical switch for profile alkalinity obs ln_palk = ', ln_palk 381 WRITE(numout,*) ' Logical switch for profile pH obs ln_pph = ', ln_pph 382 WRITE(numout,*) ' Logical switch for profile oxygen obs ln_po2 = ', ln_po2 383 WRITE(numout,*) ' Global distribution of observations ln_grid_global = ', ln_grid_global 404 WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d 405 WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d 406 WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla 407 WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst 408 WRITE(numout,*) ' Logical switch for SIC observations ln_sic = ', ln_sic 409 WRITE(numout,*) ' Logical switch for SIT observations ln_sit = ', ln_sit 410 WRITE(numout,*) ' Logical switch for FBD observations ln_fbd = ', ln_fbd 411 WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d 412 WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss 413 WRITE(numout,*) ' Logical switch for SSV observations ln_ssv = ', ln_ssv 414 WRITE(numout,*) ' Logical switch for surface total logchl obs ln_slchltot = ', ln_slchltot 415 WRITE(numout,*) ' Logical switch for surface diatom logchl obs ln_slchldia = ', ln_slchldia 416 WRITE(numout,*) ' Logical switch for surface non-diatom logchl obs ln_slchlnon = ', ln_slchlnon 417 WRITE(numout,*) ' Logical switch for surface dino logchl obs ln_slchldin = ', ln_slchldin 418 WRITE(numout,*) ' Logical switch for surface micro logchl obs ln_slchlmic = ', ln_slchlmic 419 WRITE(numout,*) ' Logical switch for surface nano logchl obs ln_slchlnan = ', ln_slchlnan 420 WRITE(numout,*) ' Logical switch for surface pico logchl obs ln_slchlpic = ', ln_slchlpic 421 WRITE(numout,*) ' Logical switch for surface total chl obs ln_schltot = ', ln_schltot 422 WRITE(numout,*) ' Logical switch for surface total log(phyC) obs ln_slphytot = ', ln_slphytot 423 WRITE(numout,*) ' Logical switch for surface diatom log(phyC) obs ln_slphydia = ', ln_slphydia 424 WRITE(numout,*) ' Logical switch for surface non-diatom log(phyC) obs ln_slphynon = ', ln_slphynon 425 WRITE(numout,*) ' Logical switch for surface SPM observations ln_sspm = ', ln_sspm 426 WRITE(numout,*) ' Logical switch for surface Kd490 observations ln_skd490 = ', ln_skd490 427 WRITE(numout,*) ' Logical switch for surface fCO2 observations ln_sfco2 = ', ln_sfco2 428 WRITE(numout,*) ' Logical switch for surface pCO2 observations ln_spco2 = ', ln_spco2 429 WRITE(numout,*) ' Logical switch for profile total logchl obs ln_plchltot = ', ln_plchltot 430 WRITE(numout,*) ' Logical switch for profile total chl obs ln_pchltot = ', ln_pchltot 431 WRITE(numout,*) ' Logical switch for profile nitrate obs ln_pno3 = ', ln_pno3 432 WRITE(numout,*) ' Logical switch for profile silicate obs ln_psi4 = ', ln_psi4 433 WRITE(numout,*) ' Logical switch for profile phosphate obs ln_ppo4 = ', ln_ppo4 434 WRITE(numout,*) ' Logical switch for profile DIC obs ln_pdic = ', ln_pdic 435 WRITE(numout,*) ' Logical switch for profile alkalinity obs ln_palk = ', ln_palk 436 WRITE(numout,*) ' Logical switch for profile pH obs ln_pph = ', ln_pph 437 WRITE(numout,*) ' Logical switch for profile oxygen obs ln_po2 = ', ln_po2 438 WRITE(numout,*) ' Global distribution of observations ln_grid_global = ', ln_grid_global 384 439 WRITE(numout,*) ' Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup 385 440 IF (ln_grid_search_lookup) & 386 441 WRITE(numout,*) ' Grid search lookup file header cn_gridsearchfile = ', cn_gridsearchfile 387 WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS rn_dobsini = ', rn_dobsini 388 WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS rn_dobsend = ', rn_dobsend 389 WRITE(numout,*) ' Type of vertical interpolation method nn_1dint = ', nn_1dint 390 WRITE(numout,*) ' Default horizontal interpolation method nn_2dint_default = ', nn_2dint_default 391 WRITE(numout,*) ' Type of horizontal interpolation method for SLA nn_2dint_sla = ', nn_2dint_sla 392 WRITE(numout,*) ' Type of horizontal interpolation method for SST nn_2dint_sst = ', nn_2dint_sst 393 WRITE(numout,*) ' Type of horizontal interpolation method for SSS nn_2dint_sss = ', nn_2dint_sss 394 WRITE(numout,*) ' Type of horizontal interpolation method for SIC nn_2dint_sic = ', nn_2dint_sic 395 WRITE(numout,*) ' Default E/W diameter of obs footprint rn_default_avglamscl = ', rn_default_avglamscl 396 WRITE(numout,*) ' Default N/S diameter of obs footprint rn_default_avgphiscl = ', rn_default_avgphiscl 397 WRITE(numout,*) ' Default obs footprint in deg [T] or m [F] ln_default_fp_indegs = ', ln_default_fp_indegs 398 WRITE(numout,*) ' SLA E/W diameter of obs footprint rn_sla_avglamscl = ', rn_sla_avglamscl 399 WRITE(numout,*) ' SLA N/S diameter of obs footprint rn_sla_avgphiscl = ', rn_sla_avgphiscl 400 WRITE(numout,*) ' SLA obs footprint in deg [T] or m [F] ln_sla_fp_indegs = ', ln_sla_fp_indegs 401 WRITE(numout,*) ' SST E/W diameter of obs footprint rn_sst_avglamscl = ', rn_sst_avglamscl 402 WRITE(numout,*) ' SST N/S diameter of obs footprint rn_sst_avgphiscl = ', rn_sst_avgphiscl 403 WRITE(numout,*) ' SST obs footprint in deg [T] or m [F] ln_sst_fp_indegs = ', ln_sst_fp_indegs 404 WRITE(numout,*) ' SIC E/W diameter of obs footprint rn_sic_avglamscl = ', rn_sic_avglamscl 405 WRITE(numout,*) ' SIC N/S diameter of obs footprint rn_sic_avgphiscl = ', rn_sic_avgphiscl 406 WRITE(numout,*) ' SIC obs footprint in deg [T] or m [F] ln_sic_fp_indegs = ', ln_sic_fp_indegs 407 WRITE(numout,*) ' Rejection of observations near land switch ln_nea = ', ln_nea 408 WRITE(numout,*) ' Rejection of obs near open bdys ln_bound_reject = ', ln_bound_reject 409 WRITE(numout,*) ' MSSH correction scheme nn_msshc = ', nn_msshc 410 WRITE(numout,*) ' MDT correction rn_mdtcorr = ', rn_mdtcorr 411 WRITE(numout,*) ' MDT cutoff for computed correction rn_mdtcutoff = ', rn_mdtcutoff 412 WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias 413 WRITE(numout,*) ' Logical switch for sst bias ln_sstbias = ', ln_sstbias 414 WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis 415 WRITE(numout,*) ' Daily average types nn_profdavtypes = ', nn_profdavtypes 416 WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight 442 WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS rn_dobsini = ', rn_dobsini 443 WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS rn_dobsend = ', rn_dobsend 444 WRITE(numout,*) ' Type of vertical interpolation method nn_1dint = ', nn_1dint 445 WRITE(numout,*) ' Default horizontal interpolation method nn_2dint_default = ', nn_2dint_default 446 WRITE(numout,*) ' Type of horizontal interpolation method for SLA nn_2dint_sla = ', nn_2dint_sla 447 WRITE(numout,*) ' Type of horizontal interpolation method for SST nn_2dint_sst = ', nn_2dint_sst 448 WRITE(numout,*) ' Type of horizontal interpolation method for SSS nn_2dint_sss = ', nn_2dint_sss 449 WRITE(numout,*) ' Type of horizontal interpolation method for SSV nn_2dint_ssv = ', nn_2dint_ssv 450 WRITE(numout,*) ' Type of horizontal interpolation method for SIC nn_2dint_sic = ', nn_2dint_sic 451 WRITE(numout,*) ' Type of horizontal interpolation method for SIT nn_2dint_sit = ', nn_2dint_sit 452 WRITE(numout,*) ' Type of horizontal interpolation method for FBD nn_2dint_fbd = ', nn_2dint_fbd 453 WRITE(numout,*) ' Default E/W diameter of obs footprint rn_default_avglamscl = ', rn_default_avglamscl 454 WRITE(numout,*) ' Default N/S diameter of obs footprint rn_default_avgphiscl = ', rn_default_avgphiscl 455 WRITE(numout,*) ' Default obs footprint in deg [T] or m [F] ln_default_fp_indegs = ', ln_default_fp_indegs 456 WRITE(numout,*) ' SLA E/W diameter of obs footprint rn_sla_avglamscl = ', rn_sla_avglamscl 457 WRITE(numout,*) ' SLA N/S diameter of obs footprint rn_sla_avgphiscl = ', rn_sla_avgphiscl 458 WRITE(numout,*) ' SLA obs footprint in deg [T] or m [F] ln_sla_fp_indegs = ', ln_sla_fp_indegs 459 WRITE(numout,*) ' SST E/W diameter of obs footprint rn_sst_avglamscl = ', rn_sst_avglamscl 460 WRITE(numout,*) ' SST N/S diameter of obs footprint rn_sst_avgphiscl = ', rn_sst_avgphiscl 461 WRITE(numout,*) ' SST obs footprint in deg [T] or m [F] ln_sst_fp_indegs = ', ln_sst_fp_indegs 462 WRITE(numout,*) ' SIC E/W diameter of obs footprint rn_sic_avglamscl = ', rn_sic_avglamscl 463 WRITE(numout,*) ' SIC N/S diameter of obs footprint rn_sic_avgphiscl = ', rn_sic_avgphiscl 464 WRITE(numout,*) ' SIC obs footprint in deg [T] or m [F] ln_sic_fp_indegs = ', ln_sic_fp_indegs 465 WRITE(numout,*) ' SIT E/W diameter of obs footprint rn_sit_avglamscl = ', rn_sit_avglamscl 466 WRITE(numout,*) ' SIT N/S diameter of obs footprint rn_sit_avgphiscl = ', rn_sit_avgphiscl 467 WRITE(numout,*) ' SIT obs footprint in deg [T] or m [F] ln_sit_fp_indegs = ', ln_sit_fp_indegs 468 WRITE(numout,*) ' FBD E/W diameter of obs footprint rn_fbd_avglamscl = ', rn_fbd_avglamscl 469 WRITE(numout,*) ' FBD N/S diameter of obs footprint rn_fbd_avgphiscl = ', rn_fbd_avgphiscl 470 WRITE(numout,*) ' FBD obs footprint in deg [T] or m [F] ln_fbd_fp_indegs = ', ln_fbd_fp_indegs 471 WRITE(numout,*) ' Rejection of observations near land switch ln_nea = ', ln_nea 472 WRITE(numout,*) ' Rejection of obs near open bdys ln_bound_reject = ', ln_bound_reject 473 WRITE(numout,*) ' MSSH correction scheme nn_msshc = ', nn_msshc 474 WRITE(numout,*) ' MDT correction rn_mdtcorr = ', rn_mdtcorr 475 WRITE(numout,*) ' MDT cutoff for computed correction rn_mdtcutoff = ', rn_mdtcutoff 476 WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias 477 WRITE(numout,*) ' Logical switch for sst bias ln_sstbias = ', ln_sstbias 478 WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis 479 WRITE(numout,*) ' Daily average types nn_profdavtypes = ', nn_profdavtypes 480 WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight 481 WRITE(numout,*) ' Logical switch for writing climat. at obs points ln_output_clim = ', ln_output_clim 482 WRITE(numout,*) ' Logical switch for time-mean of SLA ln_time_mean_sla_bkg = ', ln_time_mean_sla_bkg 417 483 ENDIF 418 484 !----------------------------------------------------------------------- … … 421 487 !----------------------------------------------------------------------- 422 488 423 nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d, ln_plchltot, & 424 & ln_pchltot, ln_pno3, ln_psi4, ln_ppo4, & 425 & ln_pdic, ln_palk, ln_pph, ln_po2 /) ) 426 nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss, & 427 & ln_slchltot, ln_slchldia, ln_slchlnon, ln_slchldin, & 428 & ln_slchlmic, ln_slchlnan, ln_slchlpic, ln_schltot, & 429 & ln_slphytot, ln_slphydia, ln_slphynon, ln_sspm, & 430 & ln_sfco2, ln_spco2 /) ) 489 nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d, ln_plchltot, & 490 & ln_pchltot, ln_pno3, ln_psi4, ln_ppo4, & 491 & ln_pdic, ln_palk, ln_pph, ln_po2 /) ) 492 nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sit, ln_fbd, & 493 & ln_sss, ln_ssv, & 494 & ln_slchltot, ln_slchldia, ln_slchlnon, ln_slchldin, & 495 & ln_slchlmic, ln_slchlnan, ln_slchlpic, ln_schltot, & 496 & ln_slphytot, ln_slphydia, ln_slphynon, ln_sspm, & 497 & ln_skd490, ln_sfco2, ln_spco2 /) ) 431 498 432 499 IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN … … 439 506 ENDIF 440 507 508 IF ( ln_output_clim .AND. ( .NOT. ln_tradmp ) ) THEN 509 IF(lwp) WRITE(numout,cform_war) 510 IF(lwp) WRITE(numout,*) ' ln_output_clim is true, but ln_tradmp is false', & 511 & ' so climatological T/S not available and will not be output' 512 nwarn = nwarn + 1 513 ln_output_clim = .FALSE. 514 ENDIF 515 441 516 IF(lwp) WRITE(numout,*) ' Number of profile obs types: ',nproftypes 442 517 IF ( nproftypes > 0 ) THEN … … 535 610 clsurffiles(jtype,:) = cn_sicfbfiles 536 611 ENDIF 612 IF (ln_sit) THEN 613 jtype = jtype + 1 614 cobstypessurf(jtype) = 'sit' 615 clsurffiles(jtype,:) = cn_sitfbfiles 616 ENDIF 617 IF (ln_fbd) THEN 618 jtype = jtype + 1 619 cobstypessurf(jtype) = 'fbd' 620 clsurffiles(jtype,:) = cn_fbdfbfiles 621 ENDIF 537 622 IF (ln_sss) THEN 538 623 jtype = jtype + 1 … … 540 625 clsurffiles(jtype,:) = cn_sssfbfiles 541 626 ENDIF 627 IF (ln_ssv) THEN 628 jtype = jtype + 1 629 cobstypessurf(jtype) = 'ssv' 630 clsurffiles(jtype,:) = cn_ssvfbfiles 631 ENDIF 542 632 IF (ln_slchltot) THEN 543 633 jtype = jtype + 1 … … 599 689 cobstypessurf(jtype) = 'sspm' 600 690 clsurffiles(jtype,:) = cn_sspmfbfiles 691 ENDIF 692 IF (ln_skd490) THEN 693 jtype = jtype + 1 694 cobstypessurf(jtype) = 'skd490' 695 clsurffiles(jtype,:) = cn_skd490fbfiles 601 696 ENDIF 602 697 IF (ln_sfco2) THEN … … 645 740 ltype_fp_indegs = ln_sic_fp_indegs 646 741 ltype_night = .FALSE. 742 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sit' ) THEN 743 IF ( nn_2dint_sit == -1 ) THEN 744 n2dint_type = nn_2dint_default 745 ELSE 746 n2dint_type = nn_2dint_sit 747 ENDIF 748 ztype_avglamscl = rn_sit_avglamscl 749 ztype_avgphiscl = rn_sit_avgphiscl 750 ltype_fp_indegs = ln_sit_fp_indegs 751 ltype_night = .FALSE. 752 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'fbd' ) THEN 753 IF ( nn_2dint_fbd == -1 ) THEN 754 n2dint_type = nn_2dint_default 755 ELSE 756 n2dint_type = nn_2dint_fbd 757 ENDIF 758 ztype_avglamscl = rn_fbd_avglamscl 759 ztype_avgphiscl = rn_fbd_avgphiscl 760 ltype_fp_indegs = ln_fbd_fp_indegs 761 ltype_night = .FALSE. 647 762 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN 648 763 IF ( nn_2dint_sss == -1 ) THEN … … 654 769 ztype_avgphiscl = rn_sss_avgphiscl 655 770 ltype_fp_indegs = ln_sss_fp_indegs 771 ltype_night = .FALSE. 772 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'ssv' ) THEN 773 IF ( nn_2dint_ssv == -1 ) THEN 774 n2dint_type = nn_2dint_default 775 ELSE 776 n2dint_type = nn_2dint_ssv 777 ENDIF 778 ztype_avglamscl = rn_ssv_avglamscl 779 ztype_avgphiscl = rn_ssv_avgphiscl 780 ltype_fp_indegs = ln_ssv_fp_indegs 656 781 ltype_night = .FALSE. 657 782 ELSE … … 715 840 DO jtype = 1, nproftypes 716 841 842 ltype_clim = .FALSE. 843 717 844 IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN 718 845 nvarsprof(jtype) = 2 719 846 nextrprof(jtype) = 1 847 IF ( ln_output_clim ) ltype_clim = .TRUE. 720 848 ALLOCATE(llvar(nvarsprof(jtype))) 849 ALLOCATE(clvars(nvarsprof(jtype))) 721 850 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zglam ) 722 851 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zgphi ) … … 734 863 nextrprof(jtype) = 2 735 864 ALLOCATE(llvar(nvarsprof(jtype))) 865 ALLOCATE(clvars(nvarsprof(jtype))) 736 866 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zglam ) 737 867 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zgphi ) … … 739 869 llvar(1) = ln_vel3d 740 870 llvar(2) = ln_vel3d 871 clvars(1) = 'UVEL' 872 clvars(2) = 'VVEL' 741 873 zglam(:,:,1) = glamu(:,:) 742 874 zglam(:,:,2) = glamv(:,:) … … 749 881 nextrprof(jtype) = 0 750 882 ALLOCATE(llvar(nvarsprof(jtype))) 883 ALLOCATE(clvars(nvarsprof(jtype))) 751 884 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zglam ) 752 885 CALL wrk_alloc( jpi, jpj, nvarsprof(jtype), zgphi ) … … 756 889 zgphi(:,:,1) = gphit(:,:) 757 890 zmask(:,:,:,1) = tmask(:,:,:) 891 IF ( TRIM(cobstypesprof(jtype)) == 'plchltot' ) THEN 892 clvars(1) = 'PLCHLTOT' 893 ELSE IF ( TRIM(cobstypesprof(jtype)) == 'pchltot' ) THEN 894 clvars(1) = 'PCHLTOT' 895 ELSE IF ( TRIM(cobstypesprof(jtype)) == 'pno3' ) THEN 896 clvars(1) = 'PNO3' 897 ELSE IF ( TRIM(cobstypesprof(jtype)) == 'psi4' ) THEN 898 clvars(1) = 'PSI4' 899 ELSE IF ( TRIM(cobstypesprof(jtype)) == 'ppo4' ) THEN 900 clvars(1) = 'PPO4' 901 ELSE IF ( TRIM(cobstypesprof(jtype)) == 'pdic' ) THEN 902 clvars(1) = 'PDIC' 903 ELSE IF ( TRIM(cobstypesprof(jtype)) == 'palk' ) THEN 904 clvars(1) = 'PALK' 905 ELSE IF ( TRIM(cobstypesprof(jtype)) == 'pph' ) THEN 906 clvars(1) = 'PPH' 907 ELSE IF ( TRIM(cobstypesprof(jtype)) == 'po2' ) THEN 908 clvars(1) = 'PO2' 909 ENDIF 758 910 ENDIF 759 911 … … 763 915 & nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, & 764 916 & rn_dobsini, rn_dobsend, llvar, & 765 & ln_ignmis, ln_s_at_t, .FALSE., &917 & ln_ignmis, ln_s_at_t, .FALSE., ltype_clim, clvars, & 766 918 & kdailyavtypes = nn_profdavtypes ) 767 919 … … 797 949 DO jtype = 1, nsurftypes 798 950 799 nvarssurf(jtype) = 1 800 nextrsurf(jtype) = 0 801 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2 951 ltype_clim = .FALSE. 952 IF ( ln_output_clim .AND. & 953 & ( ( TRIM(cobstypessurf(jtype)) == 'sst' ) .OR. & 954 & ( TRIM(cobstypessurf(jtype)) == 'sss' ) ) ) & 955 & ltype_clim = .TRUE. 956 957 IF ( (TRIM(cobstypessurf(jtype)) == 'sla') .OR. & 958 & (TRIM(cobstypessurf(jtype)) == 'sit') .OR. & 959 & (TRIM(cobstypessurf(jtype)) == 'fbd') ) THEN 960 nvarssurf(jtype) = 1 961 nextrsurf(jtype) = 2 962 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'ssv' ) THEN 963 nvarssurf(jtype) = 2 964 nextrsurf(jtype) = 0 965 ELSE 966 nvarssurf(jtype) = 1 967 nextrsurf(jtype) = 0 968 ENDIF 969 970 ALLOCATE( clvars( nvarssurf(jtype) ) ) 971 CALL wrk_alloc( jpi, jpj, nvarssurf(jtype), zglam ) 972 CALL wrk_alloc( jpi, jpj, nvarssurf(jtype), zgphi ) 973 CALL wrk_alloc( jpi, jpj, nvarssurf(jtype), zmask_surf ) 974 975 IF ( TRIM(cobstypessurf(jtype)) == 'ssv' ) THEN 976 zglam(:,:,1) = glamu(:,:) 977 zglam(:,:,2) = glamv(:,:) 978 zgphi(:,:,1) = gphiu(:,:) 979 zgphi(:,:,2) = gphiv(:,:) 980 zmask_surf(:,:,1) = umask(:,:,1) 981 zmask_surf(:,:,2) = vmask(:,:,1) 982 ELSE 983 DO jvar = 1, nvarssurf(jtype) 984 zglam(:,:,jvar) = glamt(:,:) 985 zgphi(:,:,jvar) = gphit(:,:) 986 zmask_surf(:,:,jvar) = tmask(:,:,1) 987 END DO 988 ENDIF 989 990 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 991 clvars(1) = 'SLA' 992 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN 993 clvars(1) = 'SST' 994 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN 995 clvars(1) = 'ICECONC' 996 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'fbd' ) THEN 997 clvars(1) = 'FBD' 998 ln_seaicetypes = .TRUE. 999 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sit' ) THEN 1000 clvars(1) = 'SIT' 1001 ln_seaicetypes = .TRUE. 1002 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN 1003 clvars(1) = 'SSS' 1004 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'ssv' ) THEN 1005 clvars(1) = 'UVEL' 1006 clvars(2) = 'VVEL' 1007 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchltot' ) THEN 1008 clvars(1) = 'SLCHLTOT' 1009 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchldia' ) THEN 1010 clvars(1) = 'SLCHLDIA' 1011 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchlnon' ) THEN 1012 clvars(1) = 'SLCHLNON' 1013 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchldin' ) THEN 1014 clvars(1) = 'SLCHLDIN' 1015 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchlmic' ) THEN 1016 clvars(1) = 'SLCHLMIC' 1017 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchlnan' ) THEN 1018 clvars(1) = 'SLCHLNAN' 1019 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchlpic' ) THEN 1020 clvars(1) = 'SLCHLPIC' 1021 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'schltot' ) THEN 1022 clvars(1) = 'SCHLTOT' 1023 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slphytot' ) THEN 1024 clvars(1) = 'SLPHYTOT' 1025 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slphydia' ) THEN 1026 clvars(1) = 'SLPHYDIA' 1027 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slphynon' ) THEN 1028 clvars(1) = 'SLPHYNON' 1029 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sspm' ) THEN 1030 clvars(1) = 'SSPM' 1031 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'skd490' ) THEN 1032 clvars(1) = 'SKD490' 1033 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sfco2' ) THEN 1034 clvars(1) = 'SFCO2' 1035 ELSE IF ( TRIM(cobstypessurf(jtype)) == 'spco2' ) THEN 1036 clvars(1) = 'SPCO2' 1037 ENDIF 802 1038 803 1039 !Read in surface obs types 804 CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), & 805 & clsurffiles(jtype,1:ifilessurf(jtype)), & 806 & nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 807 & rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) ) 808 809 CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject ) 810 811 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 812 CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) ) 813 IF ( ln_altbias ) & 814 & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile ) 815 ENDIF 816 817 IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN 818 jnumsstbias = 0 819 DO jfile = 1, jpmaxnfiles 820 IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) & 821 & jnumsstbias = jnumsstbias + 1 822 END DO 823 IF ( jnumsstbias == 0 ) THEN 824 CALL ctl_stop("ln_sstbias set but no bias files to read in") 1040 CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), & 1041 & clsurffiles(jtype,1:ifilessurf(jtype)), & 1042 & nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 1043 & rn_dobsini, rn_dobsend, MeanPeriodHours, ln_ignmis, .FALSE., & 1044 & llnightav(jtype), ltype_clim, ln_time_mean_sla_bkg, clvars ) 1045 1046 CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), & 1047 & jpi, jpj, & 1048 & zmask_surf, zglam, zgphi, & 1049 & ln_nea, ln_bound_reject, ln_seaicetypes ) 1050 1051 IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 1052 CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) ) 1053 IF ( ln_altbias ) & 1054 & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile ) 1055 ENDIF 1056 1057 #if defined key_cice 1058 IF ( TRIM(cobstypessurf(jtype)) == 'fbd' ) THEN 1059 CALL obs_rea_snowdepth( surfdataqc(jtype), n2dintsurf(jtype), thick_s(:,:) ) 1060 ENDIF 1061 #endif 1062 1063 IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN 1064 jnumsstbias = 0 1065 DO jfile = 1, jpmaxnfiles 1066 IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) & 1067 & jnumsstbias = jnumsstbias + 1 1068 END DO 1069 IF ( jnumsstbias == 0 ) THEN 1070 CALL ctl_stop("ln_sstbias set but no bias files to read in") 1071 825 1072 ENDIF 826 1073 827 CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), & 828 & jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) ) 829 830 ENDIF 831 832 END DO 833 1074 CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), & 1075 & jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) ) 1076 1077 ENDIF 1078 1079 DEALLOCATE( clvars ) 1080 CALL wrk_dealloc( jpi, jpj, nvarssurf(jtype), zglam ) 1081 CALL wrk_dealloc( jpi, jpj, nvarssurf(jtype), zgphi ) 1082 CALL wrk_dealloc( jpi, jpj, nvarssurf(jtype), zmask_surf ) 1083 1084 END DO 1085 834 1086 DEALLOCATE( ifilessurf, clsurffiles ) 835 1087 … … 879 1131 & frld 880 1132 #endif 881 #if defined key_cice 882 USE sbc_oce, ONLY : fr_i ! ice fraction883 #endif 884 #if defined key_top 885 USE trc, ONLY : & ! Biogeochemical state variables 886 & trn 887 #endif 888 #if defined key_hadocc 889 USE par_hadocc ! HadOCC parameters 890 USE trc, ONLY : & 891 & HADOCC_CHL, &892 & HADOCC_FCO2, &893 & HADOCC_ PCO2, &894 & HADOCC_F ILL_FLT895 USE had_bgc_const, ONLY: c2n_p896 #elif defined key_medusa 897 USE par_medusa ! MEDUSA parameters898 USE sms_medusa, ONLY: & 899 & xthetapn, &900 & xthetapd901 #if defined key_roam 902 USE sms_medusa, ONLY: &903 & f2_pco2w, & 904 & f2_fco2w, &905 & f 3_pH906 #endif 907 #elif defined key_fabm 908 USE par_fabm ! FABM parameters 909 USE fabm, ONLY: & 910 & fabm_get_interior_diagnostic_data911 #endif 912 #if defined key_spm 913 USE par_spm, ONLY: & ! Sediment parameters 914 & jp_spm 1133 #if defined key_cice 1134 USE sbc_oce, ONLY : & ! CICE variables 1135 & fr_i, & ! ice fraction 1136 & thick_i ! ice thickness 1137 #endif 1138 #if defined key_top 1139 USE trc, ONLY : & ! Biogeochemical state variables 1140 & trn 1141 #endif 1142 #if defined key_hadocc 1143 USE par_hadocc ! HadOCC parameters 1144 USE trc, ONLY : & 1145 & HADOCC_CHL, & 1146 & HADOCC_FCO2, & 1147 & HADOCC_PCO2, & 1148 & HADOCC_FILL_FLT 1149 USE had_bgc_const, ONLY: c2n_p 1150 #elif defined key_medusa 1151 USE par_medusa ! MEDUSA parameters 1152 USE sms_medusa, ONLY: & 1153 & xthetapn, & 1154 & xthetapd 1155 #if defined key_roam 1156 USE sms_medusa, ONLY: & 1157 & f2_pco2w, & 1158 & f2_fco2w, & 1159 & f3_pH 1160 #endif 1161 #elif defined key_fabm 1162 USE par_fabm ! FABM parameters 1163 #endif 1164 #if defined key_spm 1165 USE par_spm, ONLY: & ! Sediment parameters 1166 & jp_spm 915 1167 #endif 916 1168 … … 920 1172 INTEGER, INTENT(IN) :: kstp ! Current timestep 921 1173 !! * Local declarations 922 INTEGER :: idaystp ! Number of timesteps per day 923 INTEGER :: jtype ! Data loop variable 924 INTEGER :: jvar ! Variable number 925 INTEGER :: ji, jj, jk ! Loop counters 926 REAL(wp) :: tiny ! small number 927 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 928 & zprofvar ! Model values for variables in a prof ob 929 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 930 & zprofmask ! Mask associated with zprofvar 931 REAL(wp), POINTER, DIMENSION(:,:) :: & 932 & zsurfvar, & ! Model values equivalent to surface ob. 933 & zsurfmask ! Mask associated with surface variable 934 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 935 & zglam, & ! Model longitudes for prof variables 936 & zgphi ! Model latitudes for prof variables 937 LOGICAL :: llog10 ! Perform log10 transform of variable 938 #if defined key_fabm 939 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 940 & pco2_3d ! 3D pCO2 from FABM 1174 INTEGER :: idaystp ! Number of timesteps per day 1175 INTEGER :: imeanstp ! Number of timesteps for sla averaging 1176 INTEGER :: jtype ! Data loop variable 1177 INTEGER :: jvar ! Variable number 1178 INTEGER :: ji, jj, jk ! Loop counters 1179 REAL(wp) :: tiny ! small number 1180 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 1181 & zprofvar, & ! Model values for variables in a prof ob 1182 & zprofclim ! Climatology values for variables in a prof ob 1183 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 1184 & zprofmask ! Mask associated with zprofvar 1185 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 1186 & zsurfvar, & ! Model values equivalent to surface ob. 1187 & zsurfclim, & ! Climatology values for variables in a surface ob. 1188 & zsurfmask ! Mask associated with surface variable 1189 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 1190 & zglam, & ! Model longitudes 1191 & zgphi ! Model latitudes 1192 LOGICAL :: llog10 ! Perform log10 transform of variable 1193 #if defined key_fabm 1194 REAL(wp), POINTER, DIMENSION(:,:,:) :: & 1195 & fabm_3d ! 3D variable from FABM 941 1196 #endif 942 1197 … … 954 1209 !----------------------------------------------------------------------- 955 1210 956 IF ( nproftypes > 0 ) THEN 957 958 DO jtype = 1, nproftypes 959 960 ! Allocate local work arrays 961 CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar ) 962 CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask ) 963 CALL wrk_alloc( jpi, jpj, profdataqc(jtype)%nvar, zglam ) 964 CALL wrk_alloc( jpi, jpj, profdataqc(jtype)%nvar, zgphi ) 965 966 ! Defaults which might change 967 DO jvar = 1, profdataqc(jtype)%nvar 968 zprofmask(:,:,:,jvar) = tmask(:,:,:) 969 zglam(:,:,jvar) = glamt(:,:) 970 zgphi(:,:,jvar) = gphit(:,:) 971 END DO 972 973 SELECT CASE ( TRIM(cobstypesprof(jtype)) ) 974 975 CASE('prof') 976 zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem) 977 zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal) 978 979 CASE('vel') 980 zprofvar(:,:,:,1) = un(:,:,:) 981 zprofvar(:,:,:,2) = vn(:,:,:) 982 zprofmask(:,:,:,1) = umask(:,:,:) 983 zprofmask(:,:,:,2) = vmask(:,:,:) 984 zglam(:,:,1) = glamu(:,:) 985 zglam(:,:,2) = glamv(:,:) 986 zgphi(:,:,1) = gphiu(:,:) 987 zgphi(:,:,2) = gphiv(:,:) 988 989 CASE('plchltot') 990 #if defined key_hadocc 991 ! Chlorophyll from HadOCC 992 zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) 993 #elif defined key_medusa 994 ! Add non-diatom and diatom chlorophyll from MEDUSA 995 zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 996 #elif defined key_fabm 997 ! Add all chlorophyll groups from ERSEM 998 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + & 999 & trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4) 1000 #else 1001 CALL ctl_stop( ' Trying to run plchltot observation operator', & 1002 & ' but no biogeochemical model appears to have been defined' ) 1003 #endif 1004 ! Take the log10 where we can, otherwise exclude 1005 tiny = 1.0e-20 1006 WHERE(zprofvar(:,:,:,:) > tiny .AND. zprofvar(:,:,:,:) /= obfillflt ) 1007 zprofvar(:,:,:,:) = LOG10(zprofvar(:,:,:,:)) 1008 ELSEWHERE 1009 zprofvar(:,:,:,:) = obfillflt 1010 zprofmask(:,:,:,:) = 0 1011 END WHERE 1012 ! Mask out model below any excluded values, 1013 ! to avoid interpolation issues 1014 DO jvar = 1, profdataqc(jtype)%nvar 1015 DO jj = 1, jpj 1016 DO ji = 1, jpi 1017 depth_loop: DO jk = 1, jpk 1018 IF ( zprofmask(ji,jj,jk,jvar) == 0 ) THEN 1019 zprofmask(ji,jj,jk:jpk,jvar) = 0 1020 EXIT depth_loop 1021 ENDIF 1022 END DO depth_loop 1023 END DO 1024 END DO 1025 END DO 1026 1027 CASE('pchltot') 1028 #if defined key_hadocc 1029 ! Chlorophyll from HadOCC 1030 zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) 1031 #elif defined key_medusa 1032 ! Add non-diatom and diatom chlorophyll from MEDUSA 1033 zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 1034 #elif defined key_fabm 1035 ! Add all chlorophyll groups from ERSEM 1036 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + & 1037 & trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4) 1038 #else 1039 CALL ctl_stop( ' Trying to run pchltot observation operator', & 1040 & ' but no biogeochemical model appears to have been defined' ) 1041 #endif 1042 1043 CASE('pno3') 1044 #if defined key_hadocc 1045 ! Dissolved inorganic nitrogen from HadOCC 1046 zprofvar(:,:,:,1) = trn(:,:,:,jp_had_nut) 1047 #elif defined key_medusa 1048 ! Dissolved inorganic nitrogen from MEDUSA 1049 zprofvar(:,:,:,1) = trn(:,:,:,jpdin) 1050 #elif defined key_fabm 1051 ! Nitrate from ERSEM 1052 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n) 1053 #else 1054 CALL ctl_stop( ' Trying to run pno3 observation operator', & 1055 & ' but no biogeochemical model appears to have been defined' ) 1056 #endif 1057 1058 CASE('psi4') 1059 #if defined key_hadocc 1060 CALL ctl_stop( ' Trying to run psi4 observation operator', & 1061 & ' but HadOCC does not simulate silicate' ) 1062 #elif defined key_medusa 1063 ! Silicate from MEDUSA 1064 zprofvar(:,:,:,1) = trn(:,:,:,jpsil) 1065 #elif defined key_fabm 1066 ! Silicate from ERSEM 1067 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n5s) 1068 #else 1069 CALL ctl_stop( ' Trying to run psi4 observation operator', & 1070 & ' but no biogeochemical model appears to have been defined' ) 1071 #endif 1072 1073 CASE('ppo4') 1074 #if defined key_hadocc 1075 CALL ctl_stop( ' Trying to run ppo4 observation operator', & 1076 & ' but HadOCC does not simulate phosphate' ) 1077 #elif defined key_medusa 1078 CALL ctl_stop( ' Trying to run ppo4 observation operator', & 1079 & ' but MEDUSA does not simulate phosphate' ) 1080 #elif defined key_fabm 1081 ! Phosphate from ERSEM 1082 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n1p) 1083 #else 1084 CALL ctl_stop( ' Trying to run ppo4 observation operator', & 1085 & ' but no biogeochemical model appears to have been defined' ) 1086 #endif 1087 1088 CASE('pdic') 1089 #if defined key_hadocc 1090 ! Dissolved inorganic carbon from HadOCC 1091 zprofvar(:,:,:,1) = trn(:,:,:,jp_had_dic) 1092 #elif defined key_medusa 1093 ! Dissolved inorganic carbon from MEDUSA 1094 zprofvar(:,:,:,1) = trn(:,:,:,jpdic) 1095 #elif defined key_fabm 1096 ! Dissolved inorganic carbon from ERSEM 1097 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3c) 1098 #else 1099 CALL ctl_stop( ' Trying to run pdic observation operator', & 1100 & ' but no biogeochemical model appears to have been defined' ) 1101 #endif 1102 1103 CASE('palk') 1104 #if defined key_hadocc 1105 ! Alkalinity from HadOCC 1106 zprofvar(:,:,:,1) = trn(:,:,:,jp_had_alk) 1107 #elif defined key_medusa 1108 ! Alkalinity from MEDUSA 1109 zprofvar(:,:,:,1) = trn(:,:,:,jpalk) 1110 #elif defined key_fabm 1111 ! Alkalinity from ERSEM 1112 zprofvar(:,:,:,1) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta) 1113 #else 1114 CALL ctl_stop( ' Trying to run palk observation operator', & 1115 & ' but no biogeochemical model appears to have been defined' ) 1116 #endif 1117 1118 CASE('pph') 1119 #if defined key_hadocc 1120 CALL ctl_stop( ' Trying to run pph observation operator', & 1121 & ' but HadOCC has no pH diagnostic defined' ) 1122 #elif defined key_medusa && defined key_roam 1123 ! pH from MEDUSA 1124 zprofvar(:,:,:,1) = f3_pH(:,:,:) 1125 #elif defined key_fabm 1126 ! pH from ERSEM 1127 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3ph) 1128 #else 1129 CALL ctl_stop( ' Trying to run pph observation operator', & 1130 & ' but no biogeochemical model appears to have been defined' ) 1131 #endif 1132 1133 CASE('po2') 1134 #if defined key_hadocc 1135 CALL ctl_stop( ' Trying to run po2 observation operator', & 1136 & ' but HadOCC does not simulate oxygen' ) 1137 #elif defined key_medusa 1138 ! Oxygen from MEDUSA 1139 zprofvar(:,:,:,1) = trn(:,:,:,jpoxy) 1140 #elif defined key_fabm 1141 ! Oxygen from ERSEM 1142 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o2o) 1143 #else 1144 CALL ctl_stop( ' Trying to run po2 observation operator', & 1145 & ' but no biogeochemical model appears to have been defined' ) 1146 #endif 1147 1148 CASE DEFAULT 1149 CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) 1150 1151 END SELECT 1152 1153 DO jvar = 1, profdataqc(jtype)%nvar 1154 CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk, & 1155 & nit000, idaystp, jvar, & 1156 & zprofvar(:,:,:,jvar), & 1157 & fsdept(:,:,:), fsdepw(:,:,:), & 1158 & zprofmask(:,:,:,jvar), & 1159 & zglam(:,:,jvar), zgphi(:,:,jvar), & 1160 & nn_1dint, nn_2dint_default, & 1161 & kdailyavtypes = nn_profdavtypes ) 1162 END DO 1163 1164 CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar ) 1165 CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask ) 1166 CALL wrk_dealloc( jpi, jpj, profdataqc(jtype)%nvar, zglam ) 1167 CALL wrk_dealloc( jpi, jpj, profdataqc(jtype)%nvar, zgphi ) 1168 1169 END DO 1170 1171 ENDIF 1172 1173 IF ( nsurftypes > 0 ) THEN 1174 1175 !Allocate local work arrays 1176 CALL wrk_alloc( jpi, jpj, zsurfvar ) 1177 CALL wrk_alloc( jpi, jpj, zsurfmask ) 1178 #if defined key_fabm 1179 CALL wrk_alloc( jpi, jpj, jpk, pco2_3d ) 1180 #endif 1181 1182 DO jtype = 1, nsurftypes 1183 1184 !Defaults which might be changed 1185 zsurfmask(:,:) = tmask(:,:,1) 1186 llog10 = .FALSE. 1187 1188 SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 1189 CASE('sst') 1190 zsurfvar(:,:) = tsn(:,:,1,jp_tem) 1191 CASE('sla') 1192 zsurfvar(:,:) = sshn(:,:) 1193 CASE('sss') 1194 zsurfvar(:,:) = tsn(:,:,1,jp_sal) 1195 CASE('sic') 1196 IF ( kstp == 0 ) THEN 1197 IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN 1198 CALL ctl_warn( 'Sea-ice not initialised on zeroth '// & 1199 & 'time-step but some obs are valid then.' ) 1200 WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & 1201 & ' sea-ice obs will be missed' 1202 ENDIF 1203 surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + & 1204 & surfdataqc(jtype)%nsstp(1) 1205 CYCLE 1206 ELSE 1207 #if defined key_cice 1208 zsurfvar(:,:) = fr_i(:,:) 1209 #elif defined key_lim2 || defined key_lim3 1210 zsurfvar(:,:) = 1._wp - frld(:,:) 1211 #else 1212 CALL ctl_stop( ' Trying to run sea-ice observation operator', & 1213 & ' but no sea-ice model appears to have been defined' ) 1214 #endif 1215 ENDIF 1216 1217 CASE('slchltot') 1218 #if defined key_hadocc 1219 ! Surface chlorophyll from HadOCC 1220 zsurfvar(:,:) = HADOCC_CHL(:,:,1) 1221 #elif defined key_medusa 1222 ! Add non-diatom and diatom surface chlorophyll from MEDUSA 1223 zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 1224 #elif defined key_fabm 1225 ! Add all surface chlorophyll groups from ERSEM 1226 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 1227 & trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1228 #else 1229 CALL ctl_stop( ' Trying to run slchltot observation operator', & 1230 & ' but no biogeochemical model appears to have been defined' ) 1231 #endif 1232 llog10 = .TRUE. 1233 1234 CASE('slchldia') 1235 #if defined key_hadocc 1236 CALL ctl_stop( ' Trying to run slchldia observation operator', & 1237 & ' but HadOCC does not explicitly simulate diatoms' ) 1238 #elif defined key_medusa 1239 ! Diatom surface chlorophyll from MEDUSA 1240 zsurfvar(:,:) = trn(:,:,1,jpchd) 1241 #elif defined key_fabm 1242 ! Diatom surface chlorophyll from ERSEM 1243 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) 1244 #else 1245 CALL ctl_stop( ' Trying to run slchldia observation operator', & 1246 & ' but no biogeochemical model appears to have been defined' ) 1247 #endif 1248 llog10 = .TRUE. 1249 1250 CASE('slchlnon') 1251 #if defined key_hadocc 1252 CALL ctl_stop( ' Trying to run slchlnon observation operator', & 1253 & ' but HadOCC does not explicitly simulate non-diatoms' ) 1254 #elif defined key_medusa 1255 ! Non-diatom surface chlorophyll from MEDUSA 1256 zsurfvar(:,:) = trn(:,:,1,jpchn) 1257 #elif defined key_fabm 1258 ! Add all non-diatom surface chlorophyll groups from ERSEM 1259 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 1260 & trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1261 #else 1262 CALL ctl_stop( ' Trying to run slchlnon observation operator', & 1263 & ' but no biogeochemical model appears to have been defined' ) 1264 #endif 1265 llog10 = .TRUE. 1266 1267 CASE('slchldin') 1268 #if defined key_hadocc 1269 CALL ctl_stop( ' Trying to run slchldin observation operator', & 1270 & ' but HadOCC does not explicitly simulate dinoflagellates' ) 1271 #elif defined key_medusa 1272 CALL ctl_stop( ' Trying to run slchldin observation operator', & 1273 & ' but MEDUSA does not explicitly simulate dinoflagellates' ) 1274 #elif defined key_fabm 1275 ! Dinoflagellate surface chlorophyll from ERSEM 1276 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1277 #else 1278 CALL ctl_stop( ' Trying to run slchldin observation operator', & 1279 & ' but no biogeochemical model appears to have been defined' ) 1280 #endif 1281 llog10 = .TRUE. 1282 1283 CASE('slchlmic') 1284 #if defined key_hadocc 1285 CALL ctl_stop( ' Trying to run slchlmic observation operator', & 1286 & ' but HadOCC does not explicitly simulate microphytoplankton' ) 1287 #elif defined key_medusa 1288 CALL ctl_stop( ' Trying to run slchlmic observation operator', & 1289 & ' but MEDUSA does not explicitly simulate microphytoplankton' ) 1290 #elif defined key_fabm 1291 ! Add diatom and dinoflagellate surface chlorophyll from ERSEM 1292 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1293 #else 1294 CALL ctl_stop( ' Trying to run slchlmic observation operator', & 1295 & ' but no biogeochemical model appears to have been defined' ) 1296 #endif 1297 llog10 = .TRUE. 1298 1299 CASE('slchlnan') 1300 #if defined key_hadocc 1301 CALL ctl_stop( ' Trying to run slchlnan observation operator', & 1302 & ' but HadOCC does not explicitly simulate nanophytoplankton' ) 1303 #elif defined key_medusa 1304 CALL ctl_stop( ' Trying to run slchlnan observation operator', & 1305 & ' but MEDUSA does not explicitly simulate nanophytoplankton' ) 1306 #elif defined key_fabm 1307 ! Nanophytoplankton surface chlorophyll from ERSEM 1308 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) 1309 #else 1310 CALL ctl_stop( ' Trying to run slchlnan observation operator', & 1311 & ' but no biogeochemical model appears to have been defined' ) 1312 #endif 1313 llog10 = .TRUE. 1314 1315 CASE('slchlpic') 1316 #if defined key_hadocc 1317 CALL ctl_stop( ' Trying to run slchlpic observation operator', & 1318 & ' but HadOCC does not explicitly simulate picophytoplankton' ) 1319 #elif defined key_medusa 1320 CALL ctl_stop( ' Trying to run slchlpic observation operator', & 1321 & ' but MEDUSA does not explicitly simulate picophytoplankton' ) 1322 #elif defined key_fabm 1323 ! Picophytoplankton surface chlorophyll from ERSEM 1324 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) 1325 #else 1326 CALL ctl_stop( ' Trying to run slchlpic observation operator', & 1327 & ' but no biogeochemical model appears to have been defined' ) 1328 #endif 1329 llog10 = .TRUE. 1330 1331 CASE('schltot') 1332 #if defined key_hadocc 1333 ! Surface chlorophyll from HadOCC 1334 zsurfvar(:,:) = HADOCC_CHL(:,:,1) 1335 #elif defined key_medusa 1336 ! Add non-diatom and diatom surface chlorophyll from MEDUSA 1337 zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 1338 #elif defined key_fabm 1339 ! Add all surface chlorophyll groups from ERSEM 1340 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 1341 & trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1342 #else 1343 CALL ctl_stop( ' Trying to run schltot observation operator', & 1344 & ' but no biogeochemical model appears to have been defined' ) 1345 #endif 1346 1347 CASE('slphytot') 1348 #if defined key_hadocc 1349 ! Surface phytoplankton nitrogen from HadOCC multiplied by C:N ratio 1350 zsurfvar(:,:) = trn(:,:,1,jp_had_phy) * c2n_p 1351 #elif defined key_medusa 1352 ! Add non-diatom and diatom surface phytoplankton nitrogen from MEDUSA 1353 ! multiplied by C:N ratio for each 1354 zsurfvar(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd) 1355 #elif defined key_fabm 1356 ! Add all surface phytoplankton carbon groups from ERSEM 1357 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 1358 & trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c) 1359 #else 1360 CALL ctl_stop( ' Trying to run slphytot observation operator', & 1361 & ' but no biogeochemical model appears to have been defined' ) 1362 #endif 1363 llog10 = .TRUE. 1364 1365 CASE('slphydia') 1366 #if defined key_hadocc 1367 CALL ctl_stop( ' Trying to run slphydia observation operator', & 1368 & ' but HadOCC does not explicitly simulate diatoms' ) 1369 #elif defined key_medusa 1370 ! Diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio 1371 zsurfvar(:,:) = trn(:,:,1,jpphd) * xthetapd 1372 #elif defined key_fabm 1373 ! Diatom surface phytoplankton carbon from ERSEM 1374 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) 1375 #else 1376 CALL ctl_stop( ' Trying to run slphydia observation operator', & 1377 & ' but no biogeochemical model appears to have been defined' ) 1378 #endif 1379 llog10 = .TRUE. 1380 1381 CASE('slphynon') 1382 #if defined key_hadocc 1383 CALL ctl_stop( ' Trying to run slphynon observation operator', & 1384 & ' but HadOCC does not explicitly simulate non-diatoms' ) 1385 #elif defined key_medusa 1386 ! Non-diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio 1387 zsurfvar(:,:) = trn(:,:,1,jpphn) * xthetapn 1388 #elif defined key_fabm 1389 ! Add all non-diatom surface phytoplankton carbon groups from ERSEM 1390 zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 1391 & trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c) 1392 #else 1393 CALL ctl_stop( ' Trying to run slphynon observation operator', & 1394 & ' but no biogeochemical model appears to have been defined' ) 1395 #endif 1396 llog10 = .TRUE. 1397 1398 CASE('sspm') 1399 #if defined key_spm 1400 zsurfvar(:,:) = 0.0 1401 DO jn = 1, jp_spm 1402 zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn) ! sum SPM sizes 1403 END DO 1404 #else 1405 CALL ctl_stop( ' Trying to run sspm observation operator', & 1406 & ' but no spm model appears to have been defined' ) 1407 #endif 1408 1409 CASE('sfco2') 1410 #if defined key_hadocc 1411 zsurfvar(:,:) = HADOCC_FCO2(:,:) ! fCO2 from HadOCC 1412 IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. & 1413 & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN 1414 zsurfvar(:,:) = obfillflt 1415 zsurfmask(:,:) = 0 1416 CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', & 1417 & ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' ) 1418 ENDIF 1419 #elif defined key_medusa && defined key_roam 1420 zsurfvar(:,:) = f2_fco2w(:,:) 1421 #elif defined key_fabm 1422 ! First, get pCO2 from FABM 1423 pco2_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc) 1424 zsurfvar(:,:) = pco2_3d(:,:,1) 1425 ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of: 1426 ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems 1427 ! and data reduction routines, Deep-Sea Research II, 56: 512-522. 1428 ! and 1429 ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas, 1430 ! Marine Chemistry, 2: 203-215. 1431 ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so 1432 ! not explicitly included - atmospheric pressure is not necessarily available so this is 1433 ! the best assumption. 1434 ! Further, the (1-xCO2)^2 term has been neglected. This is common practice 1435 ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes) 1436 ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal 1437 ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway. 1438 zsurfvar(:,:) = zsurfvar(:,:) * EXP((-1636.75 + & 1439 & 12.0408 * (tsn(:,:,1,jp_tem)+rt0) - & 1440 & 0.0327957 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & 1441 & 0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & 1442 & 2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0))) / & 1443 & (82.0578 * (tsn(:,:,1,jp_tem)+rt0))) 1444 #else 1445 CALL ctl_stop( ' Trying to run sfco2 observation operator', & 1446 & ' but no biogeochemical model appears to have been defined' ) 1447 #endif 1448 1449 CASE('spco2') 1450 #if defined key_hadocc 1451 zsurfvar(:,:) = HADOCC_PCO2(:,:) ! pCO2 from HadOCC 1452 IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. & 1453 & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN 1454 zsurfvar(:,:) = obfillflt 1455 zsurfmask(:,:) = 0 1456 CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', & 1457 & ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' ) 1458 ENDIF 1459 #elif defined key_medusa && defined key_roam 1460 zsurfvar(:,:) = f2_pco2w(:,:) 1461 #elif defined key_fabm 1462 pco2_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc) 1463 zsurfvar(:,:) = pco2_3d(:,:,1) 1464 #else 1465 CALL ctl_stop( ' Trying to run spco2 observation operator', & 1466 & ' but no biogeochemical model appears to have been defined' ) 1467 #endif 1468 1469 CASE DEFAULT 1470 1471 CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' ) 1472 1473 END SELECT 1474 1475 IF ( llog10 ) THEN 1476 ! Take the log10 where we can, otherwise exclude 1477 tiny = 1.0e-20 1478 WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt ) 1479 zsurfvar(:,:) = LOG10(zsurfvar(:,:)) 1480 ELSEWHERE 1481 zsurfvar(:,:) = obfillflt 1482 zsurfmask(:,:) = 0 1483 END WHERE 1484 ENDIF 1485 1486 CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj, & 1487 & nit000, idaystp, zsurfvar, zsurfmask, & 1488 & n2dintsurf(jtype), llnightav(jtype), & 1489 & ravglamscl(jtype), ravgphiscl(jtype), & 1490 & lfpindegs(jtype) ) 1491 1492 END DO 1493 1494 CALL wrk_dealloc( jpi, jpj, zsurfvar ) 1495 CALL wrk_dealloc( jpi, jpj, zsurfmask ) 1496 #if defined key_fabm 1497 CALL wrk_dealloc( jpi, jpj, jpk, pco2_3d ) 1498 #endif 1499 1500 ENDIF 1211 IF ( nproftypes > 0 ) THEN 1212 1213 DO jtype = 1, nproftypes 1214 1215 ! Allocate local work arrays 1216 CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar ) 1217 CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask ) 1218 CALL wrk_alloc( jpi, jpj, profdataqc(jtype)%nvar, zglam ) 1219 CALL wrk_alloc( jpi, jpj, profdataqc(jtype)%nvar, zgphi ) 1220 CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofclim ) 1221 1222 ! Defaults which might change 1223 DO jvar = 1, profdataqc(jtype)%nvar 1224 zprofmask(:,:,:,jvar) = tmask(:,:,:) 1225 zglam(:,:,jvar) = glamt(:,:) 1226 zgphi(:,:,jvar) = gphit(:,:) 1227 zprofclim(:,:,:,jvar) = 0._wp 1228 END DO 1229 1230 SELECT CASE ( TRIM(cobstypesprof(jtype)) ) 1231 1232 CASE('prof') 1233 zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem) 1234 zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal) 1235 IF ( ln_output_clim ) THEN 1236 zprofclim(:,:,:,1) = tclim(:,:,:) 1237 zprofclim(:,:,:,2) = sclim(:,:,:) 1238 ENDIF 1239 1240 CASE('vel') 1241 zprofvar(:,:,:,1) = un(:,:,:) 1242 zprofvar(:,:,:,2) = vn(:,:,:) 1243 zprofmask(:,:,:,1) = umask(:,:,:) 1244 zprofmask(:,:,:,2) = vmask(:,:,:) 1245 zglam(:,:,1) = glamu(:,:) 1246 zglam(:,:,2) = glamv(:,:) 1247 zgphi(:,:,1) = gphiu(:,:) 1248 zgphi(:,:,2) = gphiv(:,:) 1249 1250 CASE('plchltot') 1251 #if defined key_hadocc 1252 ! Chlorophyll from HadOCC 1253 zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) 1254 #elif defined key_medusa 1255 ! Add non-diatom and diatom chlorophyll from MEDUSA 1256 zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 1257 #elif defined key_fabm 1258 ! Add all chlorophyll groups from ERSEM 1259 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + & 1260 & trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4) 1261 #else 1262 CALL ctl_stop( ' Trying to run plchltot observation operator', & 1263 & ' but no biogeochemical model appears to have been defined' ) 1264 #endif 1265 ! Take the log10 where we can, otherwise exclude 1266 tiny = 1.0e-20 1267 WHERE(zprofvar(:,:,:,:) > tiny .AND. zprofvar(:,:,:,:) /= obfillflt ) 1268 zprofvar(:,:,:,:) = LOG10(zprofvar(:,:,:,:)) 1269 ELSEWHERE 1270 zprofvar(:,:,:,:) = obfillflt 1271 zprofmask(:,:,:,:) = 0 1272 END WHERE 1273 ! Mask out model below any excluded values, 1274 ! to avoid interpolation issues 1275 DO jvar = 1, profdataqc(jtype)%nvar 1276 DO jj = 1, jpj 1277 DO ji = 1, jpi 1278 depth_loop: DO jk = 1, jpk 1279 IF ( zprofmask(ji,jj,jk,jvar) == 0 ) THEN 1280 zprofmask(ji,jj,jk:jpk,jvar) = 0 1281 EXIT depth_loop 1282 ENDIF 1283 END DO depth_loop 1284 END DO 1285 END DO 1286 END DO 1287 1288 CASE('pchltot') 1289 #if defined key_hadocc 1290 ! Chlorophyll from HadOCC 1291 zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) 1292 #elif defined key_medusa 1293 ! Add non-diatom and diatom chlorophyll from MEDUSA 1294 zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 1295 #elif defined key_fabm 1296 ! Add all chlorophyll groups from ERSEM 1297 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + & 1298 & trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4) 1299 #else 1300 CALL ctl_stop( ' Trying to run pchltot observation operator', & 1301 & ' but no biogeochemical model appears to have been defined' ) 1302 #endif 1303 1304 CASE('pno3') 1305 #if defined key_hadocc 1306 ! Dissolved inorganic nitrogen from HadOCC 1307 zprofvar(:,:,:,1) = trn(:,:,:,jp_had_nut) 1308 #elif defined key_medusa 1309 ! Dissolved inorganic nitrogen from MEDUSA 1310 zprofvar(:,:,:,1) = trn(:,:,:,jpdin) 1311 #elif defined key_fabm 1312 ! Nitrate from ERSEM 1313 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n) 1314 #else 1315 CALL ctl_stop( ' Trying to run pno3 observation operator', & 1316 & ' but no biogeochemical model appears to have been defined' ) 1317 #endif 1318 1319 CASE('psi4') 1320 #if defined key_hadocc 1321 CALL ctl_stop( ' Trying to run psi4 observation operator', & 1322 & ' but HadOCC does not simulate silicate' ) 1323 #elif defined key_medusa 1324 ! Silicate from MEDUSA 1325 zprofvar(:,:,:,1) = trn(:,:,:,jpsil) 1326 #elif defined key_fabm 1327 ! Silicate from ERSEM 1328 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n5s) 1329 #else 1330 CALL ctl_stop( ' Trying to run psi4 observation operator', & 1331 & ' but no biogeochemical model appears to have been defined' ) 1332 #endif 1333 1334 CASE('ppo4') 1335 #if defined key_hadocc 1336 CALL ctl_stop( ' Trying to run ppo4 observation operator', & 1337 & ' but HadOCC does not simulate phosphate' ) 1338 #elif defined key_medusa 1339 CALL ctl_stop( ' Trying to run ppo4 observation operator', & 1340 & ' but MEDUSA does not simulate phosphate' ) 1341 #elif defined key_fabm 1342 ! Phosphate from ERSEM 1343 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n1p) 1344 #else 1345 CALL ctl_stop( ' Trying to run ppo4 observation operator', & 1346 & ' but no biogeochemical model appears to have been defined' ) 1347 #endif 1348 1349 CASE('pdic') 1350 #if defined key_hadocc 1351 ! Dissolved inorganic carbon from HadOCC 1352 zprofvar(:,:,:,1) = trn(:,:,:,jp_had_dic) 1353 #elif defined key_medusa 1354 ! Dissolved inorganic carbon from MEDUSA 1355 zprofvar(:,:,:,1) = trn(:,:,:,jpdic) 1356 #elif defined key_fabm 1357 ! Dissolved inorganic carbon from ERSEM 1358 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3c) 1359 #else 1360 CALL ctl_stop( ' Trying to run pdic observation operator', & 1361 & ' but no biogeochemical model appears to have been defined' ) 1362 #endif 1363 1364 CASE('palk') 1365 #if defined key_hadocc 1366 ! Alkalinity from HadOCC 1367 zprofvar(:,:,:,1) = trn(:,:,:,jp_had_alk) 1368 #elif defined key_medusa 1369 ! Alkalinity from MEDUSA 1370 zprofvar(:,:,:,1) = trn(:,:,:,jpalk) 1371 #elif defined key_fabm 1372 ! Alkalinity from ERSEM 1373 zprofvar(:,:,:,1) = model%get_interior_diagnostic_data(jp_fabm_o3ta) 1374 #else 1375 CALL ctl_stop( ' Trying to run palk observation operator', & 1376 & ' but no biogeochemical model appears to have been defined' ) 1377 #endif 1378 1379 CASE('pph') 1380 #if defined key_hadocc 1381 CALL ctl_stop( ' Trying to run pph observation operator', & 1382 & ' but HadOCC has no pH diagnostic defined' ) 1383 #elif defined key_medusa && defined key_roam 1384 ! pH from MEDUSA 1385 zprofvar(:,:,:,1) = f3_pH(:,:,:) 1386 #elif defined key_fabm 1387 ! pH from ERSEM 1388 zprofvar(:,:,:,1) = model%get_interior_diagnostic_data(jp_fabm_o3ph) 1389 #else 1390 CALL ctl_stop( ' Trying to run pph observation operator', & 1391 & ' but no biogeochemical model appears to have been defined' ) 1392 #endif 1393 1394 CASE('po2') 1395 #if defined key_hadocc 1396 CALL ctl_stop( ' Trying to run po2 observation operator', & 1397 & ' but HadOCC does not simulate oxygen' ) 1398 #elif defined key_medusa 1399 ! Oxygen from MEDUSA 1400 zprofvar(:,:,:,1) = trn(:,:,:,jpoxy) 1401 #elif defined key_fabm 1402 ! Oxygen from ERSEM 1403 zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o2o) 1404 #else 1405 CALL ctl_stop( ' Trying to run po2 observation operator', & 1406 & ' but no biogeochemical model appears to have been defined' ) 1407 #endif 1408 1409 CASE DEFAULT 1410 CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) 1411 1412 END SELECT 1413 1414 DO jvar = 1, profdataqc(jtype)%nvar 1415 CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk, & 1416 & nit000, idaystp, jvar, & 1417 & zprofvar(:,:,:,jvar), & 1418 & zprofclim(:,:,:,jvar), & 1419 & fsdept(:,:,:), fsdepw(:,:,:), & 1420 & zprofmask(:,:,:,jvar), & 1421 & zglam(:,:,jvar), zgphi(:,:,jvar), & 1422 & nn_1dint, nn_2dint_default, & 1423 & kdailyavtypes = nn_profdavtypes ) 1424 END DO 1425 1426 CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar ) 1427 CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask ) 1428 CALL wrk_dealloc( jpi, jpj, profdataqc(jtype)%nvar, zglam ) 1429 CALL wrk_dealloc( jpi, jpj, profdataqc(jtype)%nvar, zgphi ) 1430 CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofclim ) 1431 1432 END DO 1433 1434 ENDIF 1435 1436 IF ( nsurftypes > 0 ) THEN 1437 1438 #if defined key_fabm 1439 CALL wrk_alloc( jpi, jpj, jpk, fabm_3d ) 1440 #endif 1441 1442 DO jtype = 1, nsurftypes 1443 1444 !Allocate local work arrays 1445 CALL wrk_alloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfvar ) 1446 CALL wrk_alloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfclim ) 1447 CALL wrk_alloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfmask ) 1448 CALL wrk_alloc( jpi, jpj, surfdataqc(jtype)%nvar, zglam ) 1449 CALL wrk_alloc( jpi, jpj, surfdataqc(jtype)%nvar, zgphi ) 1450 1451 !Defaults which might be changed 1452 DO jvar = 1, surfdataqc(jtype)%nvar 1453 zsurfmask(:,:,jvar) = tmask(:,:,1) 1454 zsurfclim(:,:,jvar) = 0._wp 1455 zglam(:,:,jvar) = glamt(:,:) 1456 zgphi(:,:,jvar) = gphit(:,:) 1457 END DO 1458 llog10 = .FALSE. 1459 1460 SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 1461 CASE('sst') 1462 zsurfvar(:,:,1) = tsn(:,:,1,jp_tem) 1463 IF ( ln_output_clim ) zsurfclim(:,:,1) = tclim(:,:,1) 1464 CASE('sla') 1465 zsurfvar(:,:,1) = sshn(:,:) 1466 CASE('sss') 1467 zsurfvar(:,:,1) = tsn(:,:,1,jp_sal) 1468 IF ( ln_output_clim ) zsurfclim(:,:,1) = sclim(:,:,1) 1469 CASE('ssv') 1470 zsurfvar(:,:,1) = un(:,:,1) 1471 zsurfvar(:,:,2) = vn(:,:,1) 1472 zsurfmask(:,:,1) = umask(:,:,1) 1473 zsurfmask(:,:,2) = vmask(:,:,1) 1474 zglam(:,:,1) = glamu(:,:) 1475 zglam(:,:,2) = glamv(:,:) 1476 zgphi(:,:,1) = gphiu(:,:) 1477 zgphi(:,:,2) = gphiv(:,:) 1478 CASE('sic') 1479 IF ( kstp == 0 ) THEN 1480 IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN 1481 CALL ctl_warn( 'Sea-ice not initialised on zeroth '// & 1482 & 'time-step but some obs are valid then.' ) 1483 WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & 1484 & ' sea-ice concentration obs will be missed' 1485 ENDIF 1486 ELSE 1487 #if defined key_cice 1488 zsurfvar(:,:,1) = fr_i(:,:) 1489 #elif defined key_lim2 || defined key_lim3 1490 zsurfvar(:,:,1) = 1._wp - frld(:,:) 1491 #else 1492 CALL ctl_stop( ' Trying to run sea-ice concentration observation operator', & 1493 & ' but no sea-ice model appears to have been defined' ) 1494 #endif 1495 ENDIF 1496 1497 CASE('sit') 1498 IF ( kstp == 0 ) THEN 1499 IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN 1500 CALL ctl_warn( 'Sea-ice not initialised on zeroth '// & 1501 & 'time-step but some obs are valid then.' ) 1502 WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & 1503 & ' sea-ice thickness obs will be missed and QC flag set to 4' 1504 ENDIF 1505 ELSE 1506 #if defined key_cice 1507 zsurfvar(:,:,1) = thick_i(:,:) 1508 #elif defined key_lim2 || defined key_lim3 1509 CALL ctl_stop( ' No sea-ice thickness observation operator defined for LIM model' ) 1510 #else 1511 CALL ctl_stop( ' Trying to run sea-ice thickness observation operator', & 1512 & ' but no sea-ice model appears to have been defined' ) 1513 #endif 1514 ENDIF 1515 1516 CASE('fbd') 1517 IF ( kstp == 0 ) THEN 1518 IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN 1519 CALL ctl_warn( 'Sea-ice not initialised on zeroth '// & 1520 & 'time-step but some obs are valid then.' ) 1521 WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & 1522 & ' sea-ice freeboard obs will be missed and QC flag set to 4' 1523 ENDIF 1524 ELSE 1525 #if defined key_cice 1526 zsurfvar(:,:) = thick_i(:,:) 1527 #elif defined key_lim2 || defined key_lim3 1528 CALL ctl_stop( ' No sea-ice freeboard observation operator defined for LIM model' ) 1529 #else 1530 CALL ctl_stop( ' Trying to run sea-ice freeboard observation operator', & 1531 & ' but no sea-ice model appears to have been defined' ) 1532 #endif 1533 ENDIF 1534 1535 CASE('slchltot') 1536 #if defined key_hadocc 1537 ! Surface chlorophyll from HadOCC 1538 zsurfvar(:,:,1) = HADOCC_CHL(:,:,1) 1539 #elif defined key_medusa 1540 ! Add non-diatom and diatom surface chlorophyll from MEDUSA 1541 zsurfvar(:,:,1) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 1542 #elif defined key_fabm 1543 ! Add all surface chlorophyll groups from ERSEM 1544 zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 1545 & trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1546 #else 1547 CALL ctl_stop( ' Trying to run slchltot observation operator', & 1548 & ' but no biogeochemical model appears to have been defined' ) 1549 #endif 1550 llog10 = .TRUE. 1551 1552 CASE('slchldia') 1553 #if defined key_hadocc 1554 CALL ctl_stop( ' Trying to run slchldia observation operator', & 1555 & ' but HadOCC does not explicitly simulate diatoms' ) 1556 #elif defined key_medusa 1557 ! Diatom surface chlorophyll from MEDUSA 1558 zsurfvar(:,:,1) = trn(:,:,1,jpchd) 1559 #elif defined key_fabm 1560 ! Diatom surface chlorophyll from ERSEM 1561 zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) 1562 #else 1563 CALL ctl_stop( ' Trying to run slchldia observation operator', & 1564 & ' but no biogeochemical model appears to have been defined' ) 1565 #endif 1566 llog10 = .TRUE. 1567 1568 CASE('slchlnon') 1569 #if defined key_hadocc 1570 CALL ctl_stop( ' Trying to run slchlnon observation operator', & 1571 & ' but HadOCC does not explicitly simulate non-diatoms' ) 1572 #elif defined key_medusa 1573 ! Non-diatom surface chlorophyll from MEDUSA 1574 zsurfvar(:,:,1) = trn(:,:,1,jpchn) 1575 #elif defined key_fabm 1576 ! Add all non-diatom surface chlorophyll groups from ERSEM 1577 zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 1578 & trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1579 #else 1580 CALL ctl_stop( ' Trying to run slchlnon observation operator', & 1581 & ' but no biogeochemical model appears to have been defined' ) 1582 #endif 1583 llog10 = .TRUE. 1584 1585 CASE('slchldin') 1586 #if defined key_hadocc 1587 CALL ctl_stop( ' Trying to run slchldin observation operator', & 1588 & ' but HadOCC does not explicitly simulate dinoflagellates' ) 1589 #elif defined key_medusa 1590 CALL ctl_stop( ' Trying to run slchldin observation operator', & 1591 & ' but MEDUSA does not explicitly simulate dinoflagellates' ) 1592 #elif defined key_fabm 1593 ! Dinoflagellate surface chlorophyll from ERSEM 1594 zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1595 #else 1596 CALL ctl_stop( ' Trying to run slchldin observation operator', & 1597 & ' but no biogeochemical model appears to have been defined' ) 1598 #endif 1599 llog10 = .TRUE. 1600 1601 CASE('slchlmic') 1602 #if defined key_hadocc 1603 CALL ctl_stop( ' Trying to run slchlmic observation operator', & 1604 & ' but HadOCC does not explicitly simulate microphytoplankton' ) 1605 #elif defined key_medusa 1606 CALL ctl_stop( ' Trying to run slchlmic observation operator', & 1607 & ' but MEDUSA does not explicitly simulate microphytoplankton' ) 1608 #elif defined key_fabm 1609 ! Add diatom and dinoflagellate surface chlorophyll from ERSEM 1610 zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1611 #else 1612 CALL ctl_stop( ' Trying to run slchlmic observation operator', & 1613 & ' but no biogeochemical model appears to have been defined' ) 1614 #endif 1615 llog10 = .TRUE. 1616 1617 CASE('slchlnan') 1618 #if defined key_hadocc 1619 CALL ctl_stop( ' Trying to run slchlnan observation operator', & 1620 & ' but HadOCC does not explicitly simulate nanophytoplankton' ) 1621 #elif defined key_medusa 1622 CALL ctl_stop( ' Trying to run slchlnan observation operator', & 1623 & ' but MEDUSA does not explicitly simulate nanophytoplankton' ) 1624 #elif defined key_fabm 1625 ! Nanophytoplankton surface chlorophyll from ERSEM 1626 zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) 1627 #else 1628 CALL ctl_stop( ' Trying to run slchlnan observation operator', & 1629 & ' but no biogeochemical model appears to have been defined' ) 1630 #endif 1631 llog10 = .TRUE. 1632 1633 CASE('slchlpic') 1634 #if defined key_hadocc 1635 CALL ctl_stop( ' Trying to run slchlpic observation operator', & 1636 & ' but HadOCC does not explicitly simulate picophytoplankton' ) 1637 #elif defined key_medusa 1638 CALL ctl_stop( ' Trying to run slchlpic observation operator', & 1639 & ' but MEDUSA does not explicitly simulate picophytoplankton' ) 1640 #elif defined key_fabm 1641 ! Picophytoplankton surface chlorophyll from ERSEM 1642 zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) 1643 #else 1644 CALL ctl_stop( ' Trying to run slchlpic observation operator', & 1645 & ' but no biogeochemical model appears to have been defined' ) 1646 #endif 1647 llog10 = .TRUE. 1648 1649 CASE('schltot') 1650 #if defined key_hadocc 1651 ! Surface chlorophyll from HadOCC 1652 zsurfvar(:,:,1) = HADOCC_CHL(:,:,1) 1653 #elif defined key_medusa 1654 ! Add non-diatom and diatom surface chlorophyll from MEDUSA 1655 zsurfvar(:,:,1) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 1656 #elif defined key_fabm 1657 ! Add all surface chlorophyll groups from ERSEM 1658 zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 1659 & trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 1660 #else 1661 CALL ctl_stop( ' Trying to run schltot observation operator', & 1662 & ' but no biogeochemical model appears to have been defined' ) 1663 #endif 1664 1665 CASE('slphytot') 1666 #if defined key_hadocc 1667 ! Surface phytoplankton nitrogen from HadOCC multiplied by C:N ratio 1668 zsurfvar(:,:,1) = trn(:,:,1,jp_had_phy) * c2n_p 1669 #elif defined key_medusa 1670 ! Add non-diatom and diatom surface phytoplankton nitrogen from MEDUSA 1671 ! multiplied by C:N ratio for each 1672 zsurfvar(:,:,1) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd) 1673 #elif defined key_fabm 1674 ! Add all surface phytoplankton carbon groups from ERSEM 1675 zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 1676 & trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c) 1677 #else 1678 CALL ctl_stop( ' Trying to run slphytot observation operator', & 1679 & ' but no biogeochemical model appears to have been defined' ) 1680 #endif 1681 llog10 = .TRUE. 1682 1683 CASE('slphydia') 1684 #if defined key_hadocc 1685 CALL ctl_stop( ' Trying to run slphydia observation operator', & 1686 & ' but HadOCC does not explicitly simulate diatoms' ) 1687 #elif defined key_medusa 1688 ! Diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio 1689 zsurfvar(:,:,1) = trn(:,:,1,jpphd) * xthetapd 1690 #elif defined key_fabm 1691 ! Diatom surface phytoplankton carbon from ERSEM 1692 zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) 1693 #else 1694 CALL ctl_stop( ' Trying to run slphydia observation operator', & 1695 & ' but no biogeochemical model appears to have been defined' ) 1696 #endif 1697 llog10 = .TRUE. 1698 1699 CASE('slphynon') 1700 #if defined key_hadocc 1701 CALL ctl_stop( ' Trying to run slphynon observation operator', & 1702 & ' but HadOCC does not explicitly simulate non-diatoms' ) 1703 #elif defined key_medusa 1704 ! Non-diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio 1705 zsurfvar(:,:,1) = trn(:,:,1,jpphn) * xthetapn 1706 #elif defined key_fabm 1707 ! Add all non-diatom surface phytoplankton carbon groups from ERSEM 1708 zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 1709 & trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c) 1710 #else 1711 CALL ctl_stop( ' Trying to run slphynon observation operator', & 1712 & ' but no biogeochemical model appears to have been defined' ) 1713 #endif 1714 llog10 = .TRUE. 1715 1716 CASE('sspm') 1717 #if defined key_spm 1718 zsurfvar(:,:,1) = 0.0 1719 DO jn = 1, jp_spm 1720 zsurfvar(:,:,1) = zsurfvar(:,:,1) + trn(:,:,1,jn) ! sum SPM sizes 1721 END DO 1722 #else 1723 CALL ctl_stop( ' Trying to run sspm observation operator', & 1724 & ' but no spm model appears to have been defined' ) 1725 #endif 1726 1727 CASE('skd490') 1728 #if defined key_hadocc 1729 CALL ctl_stop( ' Trying to run skd490 observation operator', & 1730 & ' but HadOCC does not explicitly simulate Kd490' ) 1731 #elif defined key_medusa 1732 CALL ctl_stop( ' Trying to run skd490 observation operator', & 1733 & ' but MEDUSA does not explicitly simulate Kd490' ) 1734 #elif defined key_fabm 1735 ! light_Kd_band3 diagnostic variable if using spectral optical model 1736 ! light_xEPS diagnostic variable if using standard ERSEM light model 1737 IF ( jp_fabm_kd490 /= -1 ) THEN 1738 fabm_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_kd490) 1739 ELSEIF ( jp_fabm_xeps /= -1 ) THEN 1740 fabm_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_xeps) 1741 ELSE 1742 CALL ctl_stop( ' Trying to run skd490 observation operator', & 1743 & ' but cannot access Kd490 from ERSEM' ) 1744 ENDIF 1745 zsurfvar(:,:,1) = fabm_3d(:,:,1) 1746 #else 1747 CALL ctl_stop( ' Trying to run skd490 observation operator', & 1748 & ' but no biogeochemical model appears to have been defined' ) 1749 #endif 1750 1751 CASE('sfco2') 1752 #if defined key_hadocc 1753 zsurfvar(:,:,1) = HADOCC_FCO2(:,:) ! fCO2 from HadOCC 1754 IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. & 1755 & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN 1756 zsurfvar(:,:) = obfillflt 1757 zsurfmask(:,:,1) = 0 1758 CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', & 1759 & ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' ) 1760 ENDIF 1761 #elif defined key_medusa && defined key_roam 1762 zsurfvar(:,:,1) = f2_fco2w(:,:) 1763 #elif defined key_fabm 1764 ! First, get pCO2 from FABM 1765 fabm_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_o3pc) 1766 zsurfvar(:,:,1) = fabm_3d(:,:,1) 1767 ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of: 1768 ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems 1769 ! and data reduction routines, Deep-Sea Research II, 56: 512-522. 1770 ! and 1771 ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas, 1772 ! Marine Chemistry, 2: 203-215. 1773 ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so 1774 ! not explicitly included - atmospheric pressure is not necessarily available so this is 1775 ! the best assumption. 1776 ! Further, the (1-xCO2)^2 term has been neglected. This is common practice 1777 ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes) 1778 ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal 1779 ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway. 1780 zsurfvar(:,:,1) = zsurfvar(:,:,1) * EXP((-1636.75 + & 1781 & 12.0408 * (tsn(:,:,1,jp_tem)+rt0) - & 1782 & 0.0327957 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & 1783 & 0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & 1784 & 2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0))) / & 1785 & (82.0578 * (tsn(:,:,1,jp_tem)+rt0))) 1786 #else 1787 CALL ctl_stop( ' Trying to run sfco2 observation operator', & 1788 & ' but no biogeochemical model appears to have been defined' ) 1789 #endif 1790 1791 CASE('spco2') 1792 #if defined key_hadocc 1793 zsurfvar(:,:,1) = HADOCC_PCO2(:,:) ! pCO2 from HadOCC 1794 IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. & 1795 & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN 1796 zsurfvar(:,:,1) = obfillflt 1797 zsurfmask(:,:,1) = 0 1798 CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', & 1799 & ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' ) 1800 ENDIF 1801 #elif defined key_medusa && defined key_roam 1802 zsurfvar(:,:,1) = f2_pco2w(:,:) 1803 #elif defined key_fabm 1804 fabm_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_o3pc) 1805 zsurfvar(:,:,1) = fabm_3d(:,:,1) 1806 #else 1807 CALL ctl_stop( ' Trying to run spco2 observation operator', & 1808 & ' but no biogeochemical model appears to have been defined' ) 1809 #endif 1810 1811 CASE DEFAULT 1812 1813 CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' ) 1814 1815 END SELECT 1816 1817 IF ( llog10 ) THEN 1818 ! Take the log10 where we can, otherwise exclude 1819 tiny = 1.0e-20 1820 WHERE(zsurfvar(:,:,1) > tiny .AND. zsurfvar(:,:,1) /= obfillflt ) 1821 zsurfvar(:,:,1) = LOG10(zsurfvar(:,:,1)) 1822 ELSEWHERE 1823 zsurfvar(:,:,1) = obfillflt 1824 zsurfmask(:,:,1) = 0 1825 END WHERE 1826 ENDIF 1827 1828 DO jvar = 1, surfdataqc(jtype)%nvar 1829 1830 IF ( TRIM(cobstypessurf(jtype)) == 'sla' .AND. & 1831 & ln_time_mean_sla_bkg ) THEN 1832 !Number of time-steps in meaning period 1833 imeanstp = NINT( ( MeanPeriodHours * 60. * 60. ) / rdt ) 1834 CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj, & 1835 & nit000, idaystp, jvar, & 1836 & zsurfvar(:,:,jvar), & 1837 & zsurfclim(:,:,jvar), & 1838 & zsurfmask(:,:,jvar), & 1839 & zglam(:,:,jvar), zgphi(:,:,jvar), & 1840 & n2dintsurf(jtype), llnightav(jtype), & 1841 & ravglamscl(jtype), ravgphiscl(jtype), & 1842 & lfpindegs(jtype), kmeanstp = imeanstp ) 1843 1844 ELSE 1845 CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj, & 1846 & nit000, idaystp, jvar, & 1847 & zsurfvar(:,:,jvar), & 1848 & zsurfclim(:,:,jvar), & 1849 & zsurfmask(:,:,jvar), & 1850 & zglam(:,:,jvar), zgphi(:,:,jvar), & 1851 & n2dintsurf(jtype), llnightav(jtype), & 1852 & ravglamscl(jtype), ravgphiscl(jtype), & 1853 & lfpindegs(jtype) ) 1854 ENDIF 1855 1856 END DO 1857 1858 CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfvar ) 1859 CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfclim ) 1860 CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfmask ) 1861 CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zglam ) 1862 CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zgphi ) 1863 1864 END DO 1865 #if defined key_fabm 1866 CALL wrk_dealloc( jpi, jpj, jpk, fabm_3d ) 1867 #endif 1868 1869 ENDIF 1501 1870 1502 1871 END SUBROUTINE dia_obs … … 1549 1918 & ) 1550 1919 1551 CALL obs_rotvel ( profdataqc(jtype), nn_2dint_default, zu, zv )1920 CALL obs_rotvel_pro( profdataqc(jtype), nn_2dint_default, zu, zv ) 1552 1921 1553 1922 DO jo = 1, profdataqc(jtype)%nprof … … 1582 1951 1583 1952 DO jtype = 1, nsurftypes 1953 1954 IF ( TRIM(cobstypessurf(jtype)) == 'vel' ) THEN 1955 ! For velocity data, rotate the model velocities to N/S, E/W 1956 ! using the compressed data structure. 1957 ALLOCATE( & 1958 & zu(surfdataqc(jtype)%nsurf), & 1959 & zv(surfdataqc(jtype)%nsurf) & 1960 & ) 1961 1962 CALL obs_rotvel_surf( surfdataqc(jtype), nn_2dint_default, zu, zv ) 1963 1964 DO jo = 1, surfdataqc(jtype)%nsurf 1965 surfdataqc(jtype)%rmod(jo,1) = zu(jo) 1966 surfdataqc(jtype)%rmod(jo,2) = zv(jo) 1967 END DO 1968 1969 DEALLOCATE( zu ) 1970 DEALLOCATE( zv ) 1971 END IF 1584 1972 1585 1973 CALL obs_surf_decompress( surfdataqc(jtype), & -
branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_level_search.h90
r11203 r15764 13 13 !! ** Method : Straightforward search 14 14 !! 15 !! ** Action : 15 !! ** Action : Will return level associated with T-point below the obs 16 !! depth, except when observation is in the top box will 17 !! return level 2. Also, if obs depth greater than depth 18 !! of last wet T-point (kpk-1) will return level kpk. 16 19 !! 17 20 !! History : … … 43 46 DO ji = 1, kobs 44 47 kobsk(ji) = 1 45 depk: DO jk = 2, kgrd 46 IF ( pgrddep(jk) > =pobsdep(ji) ) EXIT depk48 depk: DO jk = 2, kgrd-1 49 IF ( pgrddep(jk) > pobsdep(ji) ) EXIT depk 47 50 END DO depk 48 51 kobsk(ji) = jk -
branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
r11203 r15764 37 37 USE obs_grid, ONLY : & 38 38 & obs_level_search 39 #if defined key_cice 40 USE ice_constants, ONLY : & ! For conversion from sea ice freeboard to thickness 41 & rhos, rhoi, rhow 42 #endif 39 43 40 44 IMPLICIT NONE … … 60 64 61 65 62 SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, & 63 & kit000, kdaystp, kvar, & 64 & pvar, p gdept, pgdepw, &65 & p mask, &66 & plam, pphi, & 66 SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, & 67 & kit000, kdaystp, kvar, & 68 & pvar, pclim, & 69 & pgdept, pgdepw, pmask, & 70 & plam, pphi, & 67 71 & k1dint, k2dint, kdailyavtypes ) 68 72 … … 137 141 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 138 142 & pvar, & ! Model field for variable 143 & pclim, & ! Climatology field for variable 139 144 & pmask ! Land-sea mask for variable 140 145 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & … … 172 177 REAL(KIND=wp), DIMENSION(kpk) :: & 173 178 & zobsk, & 174 & zobs2k 179 & zobs2k, & 180 & zclm2k 175 181 REAL(KIND=wp), DIMENSION(2,2,1) :: & 176 182 & zweig1, & … … 178 184 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 179 185 & zmask, & 186 & zclim, & 180 187 & zint, & 181 188 & zinm, & … … 187 194 REAL(KIND=wp), DIMENSION(1) :: zmsk 188 195 REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 196 REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner_clim 189 197 190 198 LOGICAL :: ld_dailyav … … 262 270 & ) 263 271 272 IF ( prodatqc%lclim ) ALLOCATE( zclim(2,2,kpk,ipro) ) 273 264 274 DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 265 275 iobs = jobs - prodatqc%nprofup … … 285 295 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept, zgdept ) 286 296 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw, zgdepw ) 297 298 IF ( prodatqc%lclim ) THEN 299 CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pclim, zclim ) 300 ENDIF 287 301 288 302 ! At the end of the day also get interpolated means … … 349 363 inum_obs = iend - ista + 1 350 364 ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 365 IF ( prodatqc%lclim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) ) 351 366 352 367 DO iin=1,2 … … 358 373 & zobs2k, zgdept(iin,ijn,:,iobs), & 359 374 & zmask(iin,ijn,:,iobs)) 375 376 IF ( prodatqc%lclim ) THEN 377 CALL obs_int_z1d_spl( kpk, & 378 & zclim(iin,ijn,:,iobs), & 379 & zclm2k, zgdept(iin,ijn,:,iobs), & 380 & zmask(iin,ijn,:,iobs)) 381 ENDIF 360 382 ENDIF 361 383 … … 372 394 & zmask(iin,ijn,:,iobs)) 373 395 396 IF ( prodatqc%lclim ) THEN 397 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 398 & prodatqc%var(kvar)%vdep(ista:iend), & 399 & zclim(iin,ijn,:,iobs), & 400 & zclm2k, interp_corner_clim(iin,ijn,:), & 401 & zgdept(iin,ijn,:,iobs), & 402 & zmask(iin,ijn,:,iobs)) 403 ENDIF 374 404 ENDDO 375 405 ENDDO … … 386 416 inum_obs = iend - ista + 1 387 417 ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 418 IF ( prodatqc%lclim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) ) 388 419 DO iin=1,2 389 420 DO ijn=1,2 … … 395 426 & zmask(iin,ijn,:,iobs)) 396 427 428 IF ( prodatqc%lclim ) THEN 429 CALL obs_int_z1d_spl( kpk, & 430 & zclim(iin,ijn,:,iobs),& 431 & zclm2k, zgdept(iin,ijn,:,iobs), & 432 & zmask(iin,ijn,:,iobs)) 433 ENDIF 397 434 ENDIF 398 435 … … 408 445 & zgdept(iin,ijn,:,iobs), & 409 446 & zmask(iin,ijn,:,iobs) ) 410 447 448 IF ( prodatqc%lclim ) THEN 449 CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 450 & prodatqc%var(kvar)%vdep(ista:iend), & 451 & zclim(iin,ijn,:,iobs), & 452 & zclm2k,interp_corner_clim(iin,ijn,:), & 453 & zgdept(iin,ijn,:,iobs), & 454 & zmask(iin,ijn,:,iobs) ) 455 ENDIF 411 456 ENDDO 412 457 ENDDO … … 451 496 & prodatqc%var(kvar)%vmod(iend:iend) ) 452 497 498 IF ( prodatqc%lclim ) THEN 499 CALL obs_int_h2d( 1, 1, zweig, interp_corner_clim(:,:,ikn), & 500 & prodatqc%var(kvar)%vclm(iend:iend) ) 501 ENDIF 502 453 503 ! Set QC flag for any observations found below the bottom 454 504 ! needed as the check here is more strict than that in obs_prep … … 458 508 459 509 DEALLOCATE(interp_corner,iv_indic) 510 IF ( prodatqc%lclim ) DEALLOCATE( interp_corner_clim ) 460 511 461 512 ENDIF … … 475 526 & ) 476 527 528 IF ( prodatqc%lclim ) DEALLOCATE( zclim ) 529 477 530 ! At the end of the day also get interpolated means 478 531 IF ( ld_dailyav .AND. idayend == 0 ) THEN … … 486 539 END SUBROUTINE obs_prof_opt 487 540 488 SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, & 489 & kit000, kdaystp, psurf, psurfmask, & 490 & k2dint, ldnightav, plamscl, pphiscl, & 491 & lindegrees ) 541 SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj, & 542 & kit000, kdaystp, kvar, & 543 & psurf, pclim, psurfmask, & 544 & plam, pphi, & 545 & k2dint, ldnightav, plamscl, pphiscl, & 546 & lindegrees, kmeanstp ) 492 547 493 548 !!----------------------------------------------------------------------- … … 522 577 !! ! 15-02 (M. Martin) Combined routine for surface types 523 578 !! ! 17-03 (M. Martin) Added horizontal averaging options 579 !! ! 20-08 (M. Martin) Added surface velocity options 524 580 !!----------------------------------------------------------------------- 525 581 … … 538 594 ! (kit000-1 = restart time) 539 595 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 596 INTEGER, INTENT(IN) :: kvar ! Number of variables in surfdataqc 540 597 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 598 INTEGER, INTENT(IN), OPTIONAL :: & 599 kmeanstp ! Number of time steps for the time meaning 600 ! Averaging is triggered if present and greater than one 541 601 REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 542 602 & psurf, & ! Model surface field 603 & pclim, & ! Climatological surface field 543 604 & psurfmask ! Land-sea mask 605 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 606 & plam, & ! Model longitudes for variable 607 & pphi ! Model latitudes for variable 544 608 LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data 545 609 REAL(KIND=wp), INTENT(IN) :: & … … 559 623 INTEGER :: imodi, imodj 560 624 INTEGER :: idayend 625 INTEGER :: imeanend 561 626 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 562 627 & igrdi, & … … 569 634 REAL(wp) :: zlam 570 635 REAL(wp) :: zphi 571 REAL(wp), DIMENSION(1) :: zext, zobsmask 636 REAL(wp), DIMENSION(1) :: zext, zobsmask, zclm 572 637 REAL(wp) :: zdaystp 638 REAL(wp) :: zmeanstp 573 639 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 574 640 & zweig, & … … 577 643 & zsurfm, & 578 644 & zsurftmp, & 645 & zclim, & 579 646 & zglam, & 580 647 & zgphi, & … … 586 653 & zouttmp, & 587 654 & zmeanday ! to compute model sst in region of 24h daylight (pole) 655 656 LOGICAL :: l_timemean 588 657 589 658 !------------------------------------------------------------------------ … … 594 663 isurf = surfdataqc%nsstp(inrc) 595 664 665 l_timemean = .FALSE. 666 IF ( PRESENT( kmeanstp ) ) THEN 667 IF ( kmeanstp > 1 ) l_timemean = .TRUE. 668 ENDIF 669 596 670 ! Work out the maximum footprint size for the 597 671 ! interpolation/averaging in model grid-points - has to be even. … … 599 673 CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp ) 600 674 675 IF ( l_timemean ) THEN 676 ! Initialize time mean for first timestep 677 imeanend = MOD( kt - kit000 + 1, kmeanstp ) 678 IF (lwp) WRITE(numout,*) 'Obs time mean ', kt, kit000, kmeanstp, imeanend 679 680 ! Added kt == 0 test to catch restart case 681 IF ( ( imeanend == 1 ) .OR. ( kt == 0 ) ) THEN 682 IF (lwp) WRITE(numout,*) 'Reset surfdataqc%vdmean on time-step: ',kt 683 DO jj = 1, jpj 684 DO ji = 1, jpi 685 surfdataqc%vdmean(ji,jj,kvar) = 0.0 686 END DO 687 END DO 688 ENDIF 689 690 ! On each time-step, increment the field for computing time mean 691 IF (lwp) WRITE(numout,*)'Accumulating surfdataqc%vdmean on time-step: ',kt 692 DO jj = 1, jpj 693 DO ji = 1, jpi 694 surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) & 695 & + psurf(ji,jj) 696 END DO 697 END DO 698 699 ! Compute the time mean at the end of time period 700 IF ( imeanend == 0 ) THEN 701 zmeanstp = 1.0 / REAL( kmeanstp ) 702 IF (lwp) WRITE(numout,*)'Calculating surfdataqc%vdmean time mean on time-step: ',kt,' with weight: ',zmeanstp 703 DO jj = 1, jpj 704 DO ji = 1, jpi 705 surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) & 706 & * zmeanstp 707 END DO 708 END DO 709 ENDIF 710 ENDIF !l_timemean 601 711 602 712 IF ( ldnightav ) THEN … … 665 775 666 776 ALLOCATE( & 667 & zweig(imaxifp,imaxjfp,1), & 668 & igrdi(imaxifp,imaxjfp,isurf), & 669 & igrdj(imaxifp,imaxjfp,isurf), & 670 & zglam(imaxifp,imaxjfp,isurf), & 671 & zgphi(imaxifp,imaxjfp,isurf), & 672 & zmask(imaxifp,imaxjfp,isurf), & 673 & zsurf(imaxifp,imaxjfp,isurf), & 674 & zsurftmp(imaxifp,imaxjfp,isurf), & 675 & zglamf(imaxifp+1,imaxjfp+1,isurf), & 676 & zgphif(imaxifp+1,imaxjfp+1,isurf), & 677 & igrdip1(imaxifp+1,imaxjfp+1,isurf), & 678 & igrdjp1(imaxifp+1,imaxjfp+1,isurf) & 679 & ) 777 & zweig(imaxifp,imaxjfp,1), & 778 & igrdi(imaxifp,imaxjfp,isurf), & 779 & igrdj(imaxifp,imaxjfp,isurf), & 780 & zglam(imaxifp,imaxjfp,isurf), & 781 & zgphi(imaxifp,imaxjfp,isurf), & 782 & zmask(imaxifp,imaxjfp,isurf), & 783 & zsurf(imaxifp,imaxjfp,isurf), & 784 & zsurftmp(imaxifp,imaxjfp,isurf) & 785 & ) 786 787 IF ( k2dint > 4 ) THEN 788 ALLOCATE( & 789 & zglamf(imaxifp+1,imaxjfp+1,isurf), & 790 & zgphif(imaxifp+1,imaxjfp+1,isurf), & 791 & igrdip1(imaxifp+1,imaxjfp+1,isurf), & 792 & igrdjp1(imaxifp+1,imaxjfp+1,isurf) & 793 & ) 794 ENDIF 795 796 IF ( surfdataqc%lclim ) ALLOCATE( zclim(imaxifp,imaxjfp,isurf) ) 680 797 681 798 DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf … … 695 812 IF ( imodj > jpjglo ) imodj = jpjglo 696 813 697 igrdip1(ji+1,jj+1,iobs) = imodi 698 igrdjp1(ji+1,jj+1,iobs) = imodj 814 IF ( k2dint > 4 ) THEN 815 igrdip1(ji+1,jj+1,iobs) = imodi 816 igrdjp1(ji+1,jj+1,iobs) = imodj 817 ENDIF 699 818 700 819 IF ( ji >= 1 .AND. jj >= 1 ) THEN … … 713 832 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 714 833 & igrdi, igrdj, psurfmask, zmask ) 715 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 716 & igrdi, igrdj, psurf, zsurf ) 717 CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 718 & igrdip1, igrdjp1, glamf, zglamf ) 719 CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 720 & igrdip1, igrdjp1, gphif, zgphif ) 834 835 ! At the end of the averaging period get interpolated means 836 IF ( l_timemean ) THEN 837 IF ( imeanend == 0 ) THEN 838 ALLOCATE( zsurfm(imaxifp,imaxjfp,isurf) ) 839 IF (lwp) WRITE(numout,*)' Interpolating the time mean values on time step: ',kt 840 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 841 & igrdi, igrdj, surfdataqc%vdmean(:,:,kvar), zsurfm ) 842 ENDIF 843 ELSE 844 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 845 & igrdi, igrdj, psurf, zsurf ) 846 ENDIF 847 848 IF ( k2dint > 4 ) THEN 849 CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 850 & igrdip1, igrdjp1, glamf, zglamf ) 851 CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 852 & igrdip1, igrdjp1, gphif, zgphif ) 853 ENDIF 854 855 IF ( surfdataqc%lclim ) THEN 856 CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 857 & igrdi, igrdj, pclim, zclim ) 858 ENDIF 721 859 722 860 ! At the end of the day get interpolated means … … 758 896 zphi = surfdataqc%rphi(jobs) 759 897 760 IF ( ldnightav .AND. idayend == 0 ) THEN 761 ! Night-time averaged data 762 zsurftmp(:,:,iobs) = zsurfm(:,:,iobs) 763 ELSE 764 zsurftmp(:,:,iobs) = zsurf(:,:,iobs) 898 IF (( ldnightav .AND. idayend == 0 ) .OR. (l_timemean .AND. imeanend == 0)) THEN 899 ! Night-time or N=kmeanstp timestep averaged data 900 zsurftmp(:,:,iobs) = zsurfm(:,:,iobs) 901 ELSE 902 zsurftmp(:,:,iobs) = zsurf(:,:,iobs) 903 ENDIF 904 905 IF ( ( .NOT. l_timemean ) .OR. & 906 & ( l_timemean .AND. imeanend == 0) ) THEN 907 IF ( k2dint <= 4 ) THEN 908 909 ! Get weights to interpolate the model value to the observation point 910 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 911 & zglam(:,:,iobs), zgphi(:,:,iobs), & 912 & zmask(:,:,iobs), zweig, zobsmask ) 913 914 ! Interpolate the model value to the observation point 915 CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 916 917 IF ( surfdataqc%lclim ) THEN 918 CALL obs_int_h2d( 1, 1, zweig, zclim(:,:,iobs), zclm ) 919 ENDIF 920 921 922 ELSE 923 924 ! Get weights to average the model SLA to the observation footprint 925 CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam, zphi, & 926 & zglam(:,:,iobs), zgphi(:,:,iobs), & 927 & zglamf(:,:,iobs), zgphif(:,:,iobs), & 928 & zmask(:,:,iobs), plamscl, pphiscl, & 929 & lindegrees, zweig, zobsmask ) 930 931 ! Average the model SST to the observation footprint 932 CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 933 & zweig, zsurftmp(:,:,iobs), zext ) 934 935 IF ( surfdataqc%lclim ) THEN 936 CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 937 & zweig, zclim(:,:,iobs), zclm ) 938 ENDIF 939 940 ENDIF 941 942 IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN 943 ! ... Remove the MDT from the SSH at the observation point to get the SLA 944 surfdataqc%rext(jobs,1) = zext(1) 945 surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2) 946 ELSE 947 surfdataqc%rmod(jobs,kvar) = zext(1) 948 ENDIF 949 950 IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1) 951 952 IF ( zext(1) == obfillflt ) THEN 953 ! If the observation value is a fill value, set QC flag to bad 954 surfdataqc%nqc(jobs) = 4 955 ENDIF 956 ENDIF 957 958 IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1) 959 960 #if defined key_cice 961 IF ( TRIM(surfdataqc%cvars(1)) == 'FBD' ) THEN 962 ! Convert radar freeboard to true freeboard (add 1/4 snow depth; 1/4 based on ratio of speed of light in vacuum compared to snow (3.0e8 vs 2.4e8 m/s)) 963 surfdataqc%rext(jobs,1) = surfdataqc%robs(jobs,1) 964 surfdataqc%robs(jobs,1) = surfdataqc%rext(jobs,1) + 0.25*surfdataqc%rext(jobs,2) 965 ! If the corrected freeboard observation is outside -0.3 to 3.0 m (CPOM) then set the QC flag to bad 966 IF ((surfdataqc%robs(jobs,1) < -0.3) .OR. (surfdataqc%robs(jobs,1) > 3.0)) THEN 967 surfdataqc%nqc(jobs) = 4 968 ENDIF 969 ! Convert corrected freeboard to ice thickness following Tilling et al. (2016) 970 surfdataqc%robs(jobs,1) = (surfdataqc%robs(jobs,1)*rhow + surfdataqc%rext(jobs,2)*rhos)/ & 971 (rhow - rhoi) 765 972 ENDIF 766 973 767 IF ( k2dint <= 4 ) THEN 768 769 ! Get weights to interpolate the model value to the observation point 770 CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi, & 771 & zglam(:,:,iobs), zgphi(:,:,iobs), & 772 & zmask(:,:,iobs), zweig, zobsmask ) 773 774 ! Interpolate the model value to the observation point 775 CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 776 777 ELSE 778 779 ! Get weights to average the model SLA to the observation footprint 780 CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam, zphi, & 781 & zglam(:,:,iobs), zgphi(:,:,iobs), & 782 & zglamf(:,:,iobs), zgphif(:,:,iobs), & 783 & zmask(:,:,iobs), plamscl, pphiscl, & 784 & lindegrees, zweig, zobsmask ) 785 786 ! Average the model SST to the observation footprint 787 CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 788 & zweig, zsurftmp(:,:,iobs), zext ) 789 790 ENDIF 791 792 IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN 793 ! ... Remove the MDT from the SSH at the observation point to get the SLA 794 surfdataqc%rext(jobs,1) = zext(1) 795 surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2) 796 ELSE 797 surfdataqc%rmod(jobs,1) = zext(1) 798 ENDIF 799 800 IF ( zext(1) == obfillflt ) THEN 801 ! If the observation value is a fill value, set QC flag to bad 802 surfdataqc%nqc(jobs) = 4 803 ENDIF 974 IF ( zext(1) == obfillflt ) THEN 975 ! If the observation value is a fill value, set QC flag to bad 976 surfdataqc%nqc(jobs) = 4 977 ENDIF 978 #endif defined key_cice 804 979 805 980 END DO … … 814 989 & zmask, & 815 990 & zsurf, & 816 & zsurftmp, & 817 & zglamf, & 818 & zgphif, & 819 & igrdip1,& 820 & igrdjp1 & 991 & zsurftmp & 821 992 & ) 822 993 823 ! At the end of the day also deallocate night-time mean array 824 IF ( idayend == 0 .AND. ldnightav ) THEN 825 DEALLOCATE( & 826 & zsurfm & 994 IF ( k2dint > 4 ) THEN 995 DEALLOCATE( & 996 & zglamf, & 997 & zgphif, & 998 & igrdip1,& 999 & igrdjp1 & 827 1000 & ) 828 1001 ENDIF 829 1002 830 surfdataqc%nsurfup = surfdataqc%nsurfup + isurf 1003 IF ( surfdataqc%lclim ) DEALLOCATE( zclim ) 1004 1005 ! At the end of the day also deallocate night-time mean array 1006 IF (( idayend == 0 .AND. ldnightav ) .OR. ( imeanend == 0 .AND. l_timemean )) THEN 1007 DEALLOCATE( & 1008 & zsurfm & 1009 & ) 1010 ENDIF 1011 1012 IF ( kvar == surfdataqc%nvar ) THEN 1013 surfdataqc%nsurfup = surfdataqc%nsurfup + isurf 1014 ENDIF 831 1015 832 1016 END SUBROUTINE obs_surf_opt -
branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90
r11203 r15764 16 16 USE par_kind, ONLY : & ! Precision variables 17 17 & wp 18 USE dom_oce ! ocean space and time domain 18 19 USE in_out_manager ! I/O manager 19 20 USE obs_profiles_def ! Definitions for storage arrays for profiles … … 52 53 CONTAINS 53 54 54 SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject, & 55 kqc_cutoff ) 55 SUBROUTINE obs_pre_surf( surfdata, surfdataqc, & 56 & kpi, kpj, & 57 & zmask, pglam, pgphi, & 58 & ld_nea, ld_bound_reject, & 59 & ld_seaicetypes, kqc_cutoff ) 56 60 !!---------------------------------------------------------------------- 57 61 !! *** ROUTINE obs_pre_sla *** … … 81 85 TYPE(obs_surf), INTENT(INOUT) :: surfdata ! Full set of surface data 82 86 TYPE(obs_surf), INTENT(INOUT) :: surfdataqc ! Subset of surface data not failing screening 87 INTEGER, INTENT(IN) :: kpi, kpj ! Local domain sizes 88 REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,surfdata%nvar) :: & 89 & zmask 90 REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,surfdata%nvar) :: & 91 & pglam, & 92 & pgphi 83 93 LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land 84 94 LOGICAL, INTENT(IN) :: ld_bound_reject ! Switch for rejecting obs near the boundary 95 LOGICAL, INTENT(IN) :: ld_seaicetypes ! Switch to indicate sea ice data 85 96 INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff ! cut off for QC value 86 97 !! * Local declarations … … 93 104 INTEGER :: icycle ! Current assimilation cycle 94 105 ! Counters for observations that 95 INTEGER :: iotdobs ! - outside time domain 96 INTEGER :: iosdsobs ! - outside space domain 97 INTEGER :: ilansobs ! - within a model land cell 98 INTEGER :: inlasobs ! - close to land 99 INTEGER :: igrdobs ! - fail the grid search 100 INTEGER :: ibdysobs ! - close to open boundary 101 ! Global counters for observations that 102 INTEGER :: iotdobsmpp ! - outside time domain 103 INTEGER :: iosdsobsmpp ! - outside space domain 104 INTEGER :: ilansobsmpp ! - within a model land cell 105 INTEGER :: inlasobsmpp ! - close to land 106 INTEGER :: igrdobsmpp ! - fail the grid search 107 INTEGER :: ibdysobsmpp ! - close to open boundary 108 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 109 & llvalid ! SLA data selection 110 INTEGER :: jobs ! Obs. loop variable 111 INTEGER :: jstp ! Time loop variable 112 INTEGER :: inrc ! Time index variable 106 INTEGER :: iotdobs ! - outside time domain 107 INTEGER, DIMENSION(surfdata%nvar) :: iosdsobs ! - outside space domain 108 INTEGER, DIMENSION(surfdata%nvar) :: ilansobs ! - within a model land cell 109 INTEGER, DIMENSION(surfdata%nvar) :: inlasobs ! - close to land 110 INTEGER, DIMENSION(surfdata%nvar) :: ibdysobs ! - close to open boundary 111 INTEGER :: igrdobs ! - fail the grid search 112 ! Global counters for observations that 113 INTEGER :: iotdobsmpp ! - outside time domain 114 INTEGER, DIMENSION(surfdata%nvar) :: iosdsobsmpp ! - outside space domain 115 INTEGER, DIMENSION(surfdata%nvar) :: ilansobsmpp ! - within a model land cell 116 INTEGER, DIMENSION(surfdata%nvar) :: inlasobsmpp ! - close to land 117 INTEGER, DIMENSION(surfdata%nvar) :: ibdysobsmpp ! - close to open boundary 118 INTEGER :: igrdobsmpp ! - fail the grid search 119 120 LOGICAL, DIMENSION(:), ALLOCATABLE :: & 121 & llvalid ! SLA data selection 122 INTEGER :: jobs ! Obs. loop counter 123 INTEGER :: jvar ! Variable loop counter 124 INTEGER :: jstp ! Time loop variable 125 INTEGER :: inrc ! Time index variable 126 CHARACTER(LEN=256) :: cout1 ! Diagnostic output line 113 127 114 128 IF(lwp) WRITE(numout,*)'obs_pre_surf : Preparing the surface observations...' … … 140 154 ! ----------------------------------------------------------------------- 141 155 142 CALL obs_coo_tim( icycle, & 143 & iyea0, imon0, iday0, ihou0, imin0, & 144 & surfdata%nsurf, surfdata%nyea, surfdata%nmon, & 145 & surfdata%nday, surfdata%nhou, surfdata%nmin, & 146 & surfdata%nqc, surfdata%mstp, iotdobs ) 147 148 CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 149 150 ! ----------------------------------------------------------------------- 151 ! Check for surface data failing the grid search 152 ! ----------------------------------------------------------------------- 153 154 CALL obs_coo_grd( surfdata%nsurf, surfdata%mi, surfdata%mj, & 155 & surfdata%nqc, igrdobs ) 156 IF ( ld_seaicetypes ) THEN 157 CALL obs_coo_tim( icycle, & 158 & iyea0, imon0, iday0, ihou0, imin0, & 159 & surfdata%nsurf, surfdata%nyea, surfdata%nmon, & 160 & surfdata%nday, surfdata%nhou, surfdata%nmin, & 161 & surfdata%nqc, surfdata%mstp, iotdobs, & 162 & ld_seaicetypes = ld_seaicetypes ) 163 ELSE 164 CALL obs_coo_tim( icycle, & 165 & iyea0, imon0, iday0, ihou0, imin0, & 166 & surfdata%nsurf, surfdata%nyea, surfdata%nmon, & 167 & surfdata%nday, surfdata%nhou, surfdata%nmin, & 168 & surfdata%nqc, surfdata%mstp, iotdobs ) 169 ENDIF 170 171 CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 172 173 ! ----------------------------------------------------------------------- 174 ! Check for surface data failing the grid search 175 ! ----------------------------------------------------------------------- 176 DO jvar = 1, surfdata%nvar 177 CALL obs_coo_grd( surfdata%nsurf, surfdata%mi(:,jvar), surfdata%mj(:,jvar), & 178 & surfdata%nqc, igrdobs ) 179 END DO 156 180 157 181 CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) … … 161 185 ! ----------------------------------------------------------------------- 162 186 163 CALL obs_coo_spc_2d( surfdata%nsurf, & 164 & jpi, jpj, & 165 & surfdata%mi, surfdata%mj, & 166 & surfdata%rlam, surfdata%rphi, & 167 & glamt, gphit, & 168 & tmask(:,:,1), surfdata%nqc, & 169 & iosdsobs, ilansobs, & 170 & inlasobs, ld_nea, & 171 & ibdysobs, ld_bound_reject, & 172 & iqc_cutoff ) 173 174 CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 175 CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 176 CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 177 CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp ) 187 DO jvar = 1, surfdata%nvar 188 CALL obs_coo_spc_2d( surfdata%nsurf, & 189 & jpi, jpj, & 190 & surfdata%mi(:,jvar), surfdata%mj(:,jvar), & 191 & surfdata%rlam, surfdata%rphi, & 192 & pglam(:,:,jvar), pgphi(:,:,jvar), & 193 & zmask(:,:,jvar), surfdata%nqc(:), & 194 & iosdsobs(jvar), ilansobs(jvar), & 195 & inlasobs(jvar), ld_nea, & 196 & ibdysobs(jvar), ld_bound_reject, & 197 & iqc_cutoff ) 198 CALL obs_mpp_sum_integer( iosdsobs(jvar), iosdsobsmpp(jvar) ) 199 CALL obs_mpp_sum_integer( ilansobs(jvar), ilansobsmpp(jvar) ) 200 CALL obs_mpp_sum_integer( inlasobs(jvar), inlasobsmpp(jvar) ) 201 CALL obs_mpp_sum_integer( ibdysobs(jvar), ibdysobsmpp(jvar) ) 202 END DO 178 203 179 204 ! ----------------------------------------------------------------------- … … 204 229 ! Update the total observation counter array 205 230 206 IF(lwp) THEN 207 WRITE(numout,*) 208 WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data outside time domain = ', & 209 & iotdobsmpp 210 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data that failed grid search = ', & 211 & igrdobsmpp 212 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data outside space domain = ', & 213 & iosdsobsmpp 214 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data at land points = ', & 215 & ilansobsmpp 216 IF (ld_nea) THEN 217 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near land points (removed) = ', & 218 & inlasobsmpp 219 ELSE 220 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near land points (kept) = ', & 221 & inlasobsmpp 222 ENDIF 223 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near open boundary (removed) = ', & 224 & ibdysobsmpp 225 WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data accepted = ', & 226 & surfdataqc%nsurfmpp 227 228 WRITE(numout,*) 229 WRITE(numout,*) ' Number of observations per time step :' 230 WRITE(numout,*) 231 WRITE(numout,'(10X,A,10X,A)')'Time step',surfdataqc%cvars(1) 232 WRITE(numout,'(10X,A,5X,A)')'---------','-----------------' 233 CALL FLUSH(numout) 231 IF(lwp) THEN 232 233 DO jvar = 1, surfdataqc%nvar 234 IF ( jvar == 1 ) THEN 235 cout1=TRIM(surfdataqc%cvars(1)) 236 ELSE 237 WRITE(cout1,'(A,A1,A)') TRIM(cout1), '/', TRIM(surfdataqc%cvars(jvar)) 238 ENDIF 239 END DO 240 241 WRITE(numout,*) 242 WRITE(numout,*) ' '//TRIM(cout1)//' data outside time domain = ', & 243 & iotdobsmpp 244 WRITE(numout,*) ' Remaining '//TRIM(cout1)//' data that failed grid search = ', & 245 & igrdobsmpp 246 247 DO jvar = 1, surfdataqc%nvar 248 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data outside space domain = ', & 249 & iosdsobsmpp(jvar) 250 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data at land points = ', & 251 & ilansobsmpp(jvar) 252 IF (ld_nea) THEN 253 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data near land points (removed) = ', & 254 & inlasobsmpp(jvar) 255 ELSE 256 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data near land points (kept) = ', & 257 & inlasobsmpp(jvar) 258 ENDIF 259 WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data near open boundary (removed) = ', & 260 & ibdysobsmpp(jvar) 261 END DO 262 WRITE(numout,*) ' '//TRIM(cout1)//' data accepted = ', & 263 & surfdataqc%nsurfmpp 264 265 WRITE(numout,*) 266 WRITE(numout,*) ' Number of observations per time step:' 267 WRITE(numout,*) 268 WRITE(numout,'(10X,A,10X,A)')'Time step',TRIM(cout1) 269 WRITE(numout,'(10X,A,5X,A)')'---------','-----------------' 270 CALL FLUSH(numout) 234 271 ENDIF 235 272 … … 434 471 435 472 IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 436 CALL obs_uv_rej ( profdata, iuvchku, iuvchkv, iqc_cutoff )473 CALL obs_uv_rej_pro( profdata, iuvchku, iuvchkv, iqc_cutoff ) 437 474 CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp ) 438 475 CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp ) … … 558 595 & kobsno, & 559 596 & kobsyea, kobsmon, kobsday, kobshou, kobsmin, & 560 & kobsqc, kobsstp, kotdobs 597 & kobsqc, kobsstp, kotdobs, ld_seaicetypes ) 561 598 !!---------------------------------------------------------------------- 562 599 !! *** ROUTINE obs_coo_tim *** … … 606 643 & kobsstp ! Number of time steps up to the 607 644 ! observation time 645 LOGICAL, OPTIONAL, INTENT(IN) :: ld_seaicetypes 608 646 609 647 !! * Local declarations … … 620 658 INTEGER :: iskip 621 659 INTEGER :: idaystp 660 INTEGER :: icecount 622 661 REAL(KIND=wp) :: zminstp 623 662 REAL(KIND=wp) :: zhoustp … … 713 752 kotdobs = kotdobs + 1 714 753 CYCLE 754 ENDIF 755 756 ! Flag sea ice observations falling on initial timestep 757 IF ( PRESENT(ld_seaicetypes) ) THEN 758 759 IF ( ( kobsstp(jobs) == (nit000 - 1) ) ) THEN 760 IF (lwp) WRITE(numout,*)( 'Sea-ice not initialised on zeroth '// & 761 & 'time-step but SIT observation valid then, flagging '// & 762 'in time check subroutine obs_coo_tim.' ) 763 kobsqc(jobs) = IBSET(kobsqc(jobs),13) 764 kotdobs = kotdobs + 1 765 CYCLE 766 ENDIF 715 767 ENDIF 716 768 … … 1159 1211 #endif 1160 1212 REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 1213 & zgdept, & 1161 1214 & zgdepw 1162 1215 REAL(KIND=wp), DIMENSION(2,2,kprofno) :: & 1163 1216 & zglam, & ! Model longitude at grid points 1164 & zgphi ! Model latitude at grid points 1217 & zgphi, & ! Model latitude at grid points 1218 & zbathy ! Index of deepest wet level at grid points 1165 1219 INTEGER, DIMENSION(2,2,kprofno) :: & 1166 1220 & igrdi, & ! Grid i,j … … 1170 1224 INTEGER :: iig, ijg ! i,j of observation on model grid point. 1171 1225 INTEGER :: jobs, jobsp, jk, ji, jj 1226 REAL(KIND=wp) :: maxdepw 1172 1227 1173 1228 ! Get grid point indices … … 1220 1275 CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam ) 1221 1276 CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi ) 1277 CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, REAL(mbathy), zbathy ) 1222 1278 ! Need to know the bathy depth for each observation for sco 1223 1279 CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, fsdepw(:,:,:), & 1224 1280 & zgdepw ) 1281 CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, fsdept(:,:,:), & 1282 & zgdept ) 1225 1283 1226 1284 DO jobs = 1, kprofno … … 1258 1316 DO jobsp = kpstart(jobs), kpend(jobs) 1259 1317 1318 ! Calculate max T and W depths of 2x2 grid 1319 maxdepw=zgdepw(1,1,NINT(zbathy(1,1,jobs))+1,jobs) 1320 DO jj = 1, 2 1321 DO ji = 1, 2 1322 IF ( zgdepw(ji,jj,NINT(zbathy(ji,jj,jobs))+1,jobs) > maxdepw ) THEN 1323 maxdepw = zgdepw(ji,jj,NINT(zbathy(ji,jj,jobs))+1,jobs) 1324 END IF 1325 END DO 1326 END DO 1327 1260 1328 ! Flag if the observation falls outside the model spatial domain 1261 1329 IF ( ( pobslam(jobs) < -180. ) & … … 1264 1332 & .OR. ( pobsphi(jobs) > 90. ) & 1265 1333 & .OR. ( pobsdep(jobsp) < 0.0 ) & 1266 & .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN1334 & .OR. ( pobsdep(jobsp) >= maxdepw ) ) THEN 1267 1335 kobsqc(jobsp) = IBSET(kobsqc(jobsp),11) 1268 1336 kosdobs = kosdobs + 1 … … 1327 1395 ENDIF 1328 1396 1329 ! Set observation depth equal to that of the first model depth1330 IF ( pobsdep(jobsp) <= pdep(1) ) THEN1331 pobsdep(jobsp) = pdep(1)1332 ENDIF1333 1334 1397 #if defined key_bdy 1335 1398 ! Flag if the observation falls close to the boundary rim … … 1405 1468 END SUBROUTINE obs_pro_rej 1406 1469 1407 SUBROUTINE obs_uv_rej ( profdata, knumu, knumv, kqc_cutoff )1408 !!---------------------------------------------------------------------- 1409 !! *** ROUTINE obs_uv_rej ***1470 SUBROUTINE obs_uv_rej_pro( profdata, knumu, knumv, kqc_cutoff ) 1471 !!---------------------------------------------------------------------- 1472 !! *** ROUTINE obs_uv_rej_pro *** 1410 1473 !! 1411 1474 !! ** Purpose : Reject u if v is rejected and vice versa … … 1439 1502 & ( profdata%npvend(jprof,1) /= profdata%npvend(jprof,2) ) ) THEN 1440 1503 1441 CALL ctl_stop('U,V profiles inconsistent in obs_uv_rej ')1504 CALL ctl_stop('U,V profiles inconsistent in obs_uv_rej_pro') 1442 1505 RETURN 1443 1506 … … 1461 1524 END DO 1462 1525 1463 END SUBROUTINE obs_uv_rej 1526 END SUBROUTINE obs_uv_rej_pro 1464 1527 1465 1528 END MODULE obs_prep -
branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_profiles_def.F90
r11203 r15764 72 72 & vdep, & !: Depth coordinate of profile data 73 73 & vobs, & !: Profile data 74 & vmod !: Model counterpart of the profile data vector 74 & vmod, & !: Model counterpart of the profile data vector 75 & vclm !: Climatological counterpart of the profile data vector 75 76 76 77 REAL(KIND=wp), POINTER, DIMENSION(:,:) :: & … … 101 102 INTEGER :: npk 102 103 INTEGER :: nprofup !: Observation counter used in obs_oper 104 105 LOGICAL :: lclim !: Climatology will be calculated for this structure 103 106 104 107 ! Bookkeeping arrays with sizes equal to number of variables … … 198 201 199 202 SUBROUTINE obs_prof_alloc( prof, kvar, kext, kprof, & 200 & ko3dt, kstp, kpi, kpj, kpk )203 & ko3dt, kstp, kpi, kpj, kpk, ldclim ) 201 204 !!---------------------------------------------------------------------- 202 205 !! *** ROUTINE obs_prof_alloc *** … … 221 224 INTEGER, INTENT(IN) :: kpj 222 225 INTEGER, INTENT(IN) :: kpk 226 LOGICAL, INTENT(IN) :: ldclim 223 227 224 228 !!* Local variables … … 236 240 prof%npj = kpj 237 241 prof%npk = kpk 242 243 prof%lclim = ldclim 238 244 239 245 ! Allocate arrays of size number of variables … … 503 509 & ) 504 510 ENDIF 511 IF (prof%lclim) THEN 512 ALLOCATE( & 513 & prof%var(kvar)%vclm(kobs) & 514 & ) 515 ENDIF 505 516 506 517 END SUBROUTINE obs_prof_alloc_var … … 537 548 DEALLOCATE( & 538 549 & prof%var(kvar)%vext & 550 & ) 551 ENDIF 552 IF (prof%lclim) THEN 553 DEALLOCATE( & 554 & prof%var(kvar)%vclm & 539 555 & ) 540 556 ENDIF … … 630 646 & inprof, invpro, & 631 647 & prof%nstp, prof%npi, & 632 & prof%npj, prof%npk ) 648 & prof%npj, prof%npk, & 649 & prof%lclim ) 633 650 ENDIF 634 651 … … 745 762 & prof%var(jvar)%vext(jj,jext) 746 763 END DO 764 IF (newprof%lclim) THEN 765 newprof%var(jvar)%vclm(invpro(jvar)) = & 766 & prof%var(jvar)%vclm(jj) 767 ENDIF 747 768 748 769 ! nvind is the index of the original variable data … … 870 891 & prof%var(jvar)%vext(jj,jext) 871 892 END DO 893 IF (prof%lclim) THEN 894 oldprof%var(jvar)%vclm(jl) = prof%var(jvar)%vclm(jj) 895 ENDIF 872 896 873 897 END DO -
branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90
r11203 r15764 46 46 & kvars, kextr, kstp, ddobsini, ddobsend, & 47 47 & ldvar, ldignmis, ldsatt, & 48 & ldmod, kdailyavtypes )48 & ldmod, ldclim, cdvars, kdailyavtypes ) 49 49 !!--------------------------------------------------------------------- 50 50 !! … … 78 78 LOGICAL, INTENT(IN) :: ldsatt ! Compute salinity at all temperature points 79 79 LOGICAL, INTENT(IN) :: ldmod ! Initialize model from input data 80 LOGICAL, INTENT(IN) :: ldclim ! Set flag to show climatology will be output 80 81 REAL(dp), INTENT(IN) :: ddobsini ! Obs. ini time in YYYYMMDD.HHMMSS 81 82 REAL(dp), INTENT(IN) :: ddobsend ! Obs. end time in YYYYMMDD.HHMMSS 83 CHARACTER(len=8), DIMENSION(kvars), INTENT(IN) :: cdvars 82 84 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 83 85 & kdailyavtypes ! Types of daily average observations … … 86 88 CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_prof' 87 89 CHARACTER(len=8) :: clrefdate 88 CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars 90 CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvarsin 89 91 INTEGER :: jvar 90 92 INTEGER :: ji … … 220 222 ENDIF 221 223 222 IF ( jj == 1 ) THEN 223 ALLOCATE( clvars( inpfiles(jj)%nvar ) ) 224 DO ji = 1, inpfiles(jj)%nvar 225 clvars(ji) = inpfiles(jj)%cname(ji) 226 END DO 224 IF ( jj == 1 ) THEN 225 ALLOCATE( clvarsin( inpfiles(jj)%nvar ) ) 226 DO ji = 1, inpfiles(jj)%nvar 227 clvarsin(ji) = inpfiles(jj)%cname(ji) 228 IF ( clvarsin(ji) /= cdvars(ji) ) THEN 229 CALL ctl_stop( 'Feedback file variables do not match', & 230 & ' expected variable names for this type' ) 231 ENDIF 232 END DO 227 233 ELSE 228 DO ji = 1, inpfiles(jj)%nvar 229 IF ( inpfiles(jj)%cname(ji) /= clvars (ji) ) THEN230 CALL ctl_stop( 'Feedback file variables not consistent', & 231 & ' with previous files for this type' ) 232 ENDIF 233 END DO 234 DO ji = 1, inpfiles(jj)%nvar 235 IF ( inpfiles(jj)%cname(ji) /= clvarsin(ji) ) THEN 236 CALL ctl_stop( 'Feedback file variables not consistent', & 237 & ' with previous files for this type' ) 238 ENDIF 239 END DO 234 240 ENDIF 235 241 … … 500 506 ENDIF 501 507 CALL obs_prof_alloc( profdata, kvars, kextr, iprof, iv3dt, & 502 & kstp, jpi, jpj, jpk )508 & kstp, jpi, jpj, jpk, ldclim ) 503 509 504 510 ! * Read obs/positions, QC, all variable and assign to profdata … … 506 512 profdata%nprof = 0 507 513 profdata%nvprot(:) = 0 508 profdata%cvars(:) = clvars (:)514 profdata%cvars(:) = clvarsin(:) 509 515 iprof = 0 510 516 … … 696 702 profdata%var(jvar)%vmod(ivart(jvar)) = & 697 703 & inpfiles(jj)%padd(ij,ji,1,jvar) 704 ENDIF 705 IF ( profdata%lclim ) THEN 706 profdata%var(jvar)%vclm(ivart(jvar)) = fbrmdi 698 707 ENDIF 699 708 ! Count number of profile var1 data as function of type … … 805 814 ! Deallocate temporary data 806 815 !----------------------------------------------------------------------- 807 DEALLOCATE( ifileidx, iprofidx, zdat, clvars )816 DEALLOCATE( ifileidx, iprofidx, zdat, clvarsin ) 808 817 809 818 !----------------------------------------------------------------------- -
branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90
r11202 r15764 39 39 40 40 SUBROUTINE obs_rea_surf( surfdata, knumfiles, cdfilenames, & 41 & kvars, kextr, kstp, ddobsini, ddobsend, &42 & ldignmis, ldmod, ldnightav )41 & kvars, kextr, kstp, ddobsini, ddobsend, MeanPeriodHours, & 42 & ldignmis, ldmod, ldnightav, ldclim, ln_time_mean_sla_bkg, cdvars ) 43 43 !!--------------------------------------------------------------------- 44 44 !! … … 71 71 LOGICAL, INTENT(IN) :: ldmod ! Initialize model from input data 72 72 LOGICAL, INTENT(IN) :: ldnightav ! Observations represent a night-time average 73 LOGICAL, INTENT(IN) :: ldclim ! Will include climatology at obs points. 74 LOGICAL, INTENT(IN) :: ln_time_mean_sla_bkg ! Will reset times to end of averaging period. 73 75 REAL(dp), INTENT(IN) :: ddobsini ! Obs. ini time in YYYYMMDD.HHMMSS 74 76 REAL(dp), INTENT(IN) :: ddobsend ! Obs. end time in YYYYMMDD.HHMMSS 77 REAL(wp), INTENT(IN) :: MeanPeriodHours ! Averaging period in hours 78 CHARACTER(len=8), DIMENSION(kvars), INTENT(IN) :: cdvars 75 79 76 80 !! * Local declarations 77 81 CHARACTER(LEN=11), PARAMETER :: cpname='obs_rea_surf' 78 82 CHARACTER(len=8) :: clrefdate 79 CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars 83 CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvarsin 80 84 INTEGER :: ji 81 85 INTEGER :: jj 82 86 INTEGER :: jk 87 INTEGER :: jvar 83 88 INTEGER :: iflag 84 89 INTEGER :: inobf … … 103 108 & ityp, & 104 109 & itypmpp 105 INTEGER, DIMENSION(: ), ALLOCATABLE :: &110 INTEGER, DIMENSION(:,:), ALLOCATABLE :: & 106 111 & iobsi, & 107 112 & iobsj, & 108 & iproc, & 113 & iproc 114 INTEGER, DIMENSION(:), ALLOCATABLE :: & 109 115 & iindx, & 110 116 & ifileidx, & … … 122 128 TYPE(obfbdata), POINTER, DIMENSION(:) :: & 123 129 & inpfiles 124 130 CHARACTER(LEN=256) :: cout1 ! Diagnostic output line 131 125 132 ! Local initialization 126 133 iobs = 0 … … 181 188 & ldgrid = .TRUE. ) 182 189 190 IF ( inpfiles(jj)%nvar /= kvars ) THEN 191 CALL ctl_stop( 'Feedback format error: ', & 192 & ' unexpected number of vars in feedback file' ) 193 ENDIF 194 183 195 IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN 184 196 CALL ctl_stop( 'Model not in input data' ) … … 187 199 188 200 IF ( jj == 1 ) THEN 189 ALLOCATE( clvars ( inpfiles(jj)%nvar ) )201 ALLOCATE( clvarsin( inpfiles(jj)%nvar ) ) 190 202 DO ji = 1, inpfiles(jj)%nvar 191 clvars(ji) = inpfiles(jj)%cname(ji) 203 clvarsin(ji) = inpfiles(jj)%cname(ji) 204 IF ( clvarsin(ji) /= cdvars(ji) ) THEN 205 CALL ctl_stop( 'Feedback file variables do not match', & 206 & ' expected variable names for this type' ) 207 ENDIF 192 208 END DO 193 209 ELSE 194 210 DO ji = 1, inpfiles(jj)%nvar 195 IF ( inpfiles(jj)%cname(ji) /= clvars (ji) ) THEN211 IF ( inpfiles(jj)%cname(ji) /= clvarsin(ji) ) THEN 196 212 CALL ctl_stop( 'Feedback file variables not consistent', & 197 213 & ' with previous files for this type' ) … … 265 281 266 282 IF ( inpfiles(jj)%nobs > 0 ) THEN 267 inpfiles(jj)%iproc = -1 268 inpfiles(jj)%iobsi = -1 269 inpfiles(jj)%iobsj = -1 270 ENDIF 283 inpfiles(jj)%iproc(:,:) = -1 284 inpfiles(jj)%iobsi(:,:) = -1 285 inpfiles(jj)%iobsj(:,:) = -1 286 ENDIF 287 288 !If SLA observations are representing a time mean then set the time 289 !of the obs to the end of that meaning period relative to the start of the run 290 IF ( ln_time_mean_sla_bkg .AND. ( TRIM( clvarsin(1) ) == 'SLA' ) ) THEN 291 DO ji = 1, inpfiles(jj)%nobs 292 ! Only do this for obs within time window 293 IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. & 294 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN 295 inpfiles(jj)%ptim(ji) = & 296 & djulini(jj) + (MeanPeriodHours/24.) 297 ENDIF 298 END DO 299 ENDIF 300 271 301 inowin = 0 272 302 DO ji = 1, inpfiles(jj)%nobs … … 278 308 ALLOCATE( zlam(inowin) ) 279 309 ALLOCATE( zphi(inowin) ) 280 ALLOCATE( iobsi(inowin ) )281 ALLOCATE( iobsj(inowin ) )282 ALLOCATE( iproc(inowin ) )310 ALLOCATE( iobsi(inowin,kvars) ) 311 ALLOCATE( iobsj(inowin,kvars) ) 312 ALLOCATE( iproc(inowin,kvars) ) 283 313 inowin = 0 284 314 DO ji = 1, inpfiles(jj)%nobs … … 291 321 END DO 292 322 293 CALL obs_grid_search( inowin, zlam, zphi, iobsi, iobsj, iproc, 'T' ) 323 ! Assume anything other than velocity is on T grid 324 IF ( TRIM( inpfiles(jj)%cname(1) ) == 'UVEL' ) THEN 325 CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,1), iobsj(:,1), & 326 & iproc(:,1), 'U' ) 327 CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,2), iobsj(:,2), & 328 & iproc(:,2), 'V' ) 329 ELSE 330 CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,1), iobsj(:,1), & 331 & iproc(:,1), 'T' ) 332 IF ( kvars > 1 ) THEN 333 DO jvar = 2, kvars 334 iobsi(:,jvar) = iobsi(:,1) 335 iobsj(:,jvar) = iobsj(:,1) 336 iproc(:,jvar) = iproc(:,1) 337 END DO 338 ENDIF 339 ENDIF 294 340 295 341 inowin = 0 … … 298 344 & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN 299 345 inowin = inowin + 1 300 inpfiles(jj)%iproc(ji,1) = iproc(inowin) 301 inpfiles(jj)%iobsi(ji,1) = iobsi(inowin) 302 inpfiles(jj)%iobsj(ji,1) = iobsj(inowin) 346 DO jvar = 1, kvars 347 inpfiles(jj)%iproc(ji,jvar) = iproc(inowin,jvar) 348 inpfiles(jj)%iobsi(ji,jvar) = iobsi(inowin,jvar) 349 inpfiles(jj)%iobsj(ji,jvar) = iobsj(inowin,jvar) 350 END DO 303 351 ENDIF 304 352 END DO … … 359 407 & iindx ) 360 408 361 CALL obs_surf_alloc( surfdata, iobs, kvars, iextr, kstp, jpi, jpj )409 CALL obs_surf_alloc( surfdata, iobs, kvars, iextr, kstp, jpi, jpj, ldclim ) 362 410 363 411 ! Read obs/positions, QC, all variable and assign to surfdata 364 412 365 413 iobs = 0 366 surfdata%cvars(:) = clvars (:)414 surfdata%cvars(:) = clvarsin(:) 367 415 IF ( ldmod .AND. ( TRIM( surfdata%cvars(1) ) == 'SLA' ) ) THEN 368 416 surfdata%cext(1) = 'SSH' 369 417 surfdata%cext(2) = 'MDT' 418 ENDIF 419 IF ( ldmod .AND. ( TRIM( surfdata%cvars(1) ) == 'FBD' ) ) THEN 420 surfdata%cext(1) = 'freeboard' 421 surfdata%cext(2) = 'thick_s' 370 422 ENDIF 371 423 IF ( iextr > kextr ) surfdata%cext(iextr) = 'STD' … … 417 469 418 470 ! Coordinate search parameters 419 surfdata%mi (iobs) = inpfiles(jj)%iobsi(ji,1) 420 surfdata%mj (iobs) = inpfiles(jj)%iobsj(ji,1) 421 471 DO jvar = 1, kvars 472 surfdata%mi(iobs,jvar) = inpfiles(jj)%iobsi(ji,jvar) 473 surfdata%mj(iobs,jvar) = inpfiles(jj)%iobsj(ji,jvar) 474 END DO 475 422 476 ! WMO number 423 477 surfdata%cwmo(iobs) = inpfiles(jj)%cdwmo(ji) … … 449 503 450 504 ! Observed value 451 surfdata%robs(iobs,1) = inpfiles(jj)%pob(1,ji,1) 505 DO jvar = 1, kvars 506 surfdata%robs(iobs,jvar) = inpfiles(jj)%pob(1,ji,jvar) 507 END DO 508 IF ( TRIM(surfdata%cvars(1)) == 'FBD' ) THEN 509 surfdata%rext(iobs,1) = inpfiles(jj)%pob(1,ji,1) 510 surfdata%rext(iobs,2) = fbrmdi 511 ENDIF 452 512 453 513 ! Model and MDT is set to fbrmdi unless read from file 454 514 IF ( ldmod ) THEN 455 surfdata%rmod(iobs,1) = inpfiles(jj)%padd(1,ji,1,1) 515 DO jvar = 1, kvars 516 surfdata%rmod(iobs,jvar) = inpfiles(jj)%padd(1,ji,1,jvar) 517 END DO 456 518 IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) THEN 457 519 surfdata%rext(iobs,1) = inpfiles(jj)%padd(1,ji,2,1) … … 459 521 ENDIF 460 522 ELSE 461 surfdata%rmod(iobs,1) = fbrmdi 523 DO jvar = 1, kvars 524 surfdata%rmod(iobs,jvar) = fbrmdi 525 END DO 462 526 IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) surfdata%rext(iobs,:) = fbrmdi 463 527 ENDIF 464 528 529 ! Initialise climatology if set 530 IF ( surfdata%lclim ) THEN 531 DO jvar = 1, kvars 532 surfdata%rclm(iobs,jvar) = fbrmdi 533 END DO 534 ENDIF 535 465 536 ! STD (obs error standard deviation) read from file and passed through obs operator 466 537 IF ( iadd_std(jj) /= -1 ) THEN … … 483 554 !----------------------------------------------------------------------- 484 555 IF (lwp) THEN 485 556 DO jvar = 1, surfdata%nvar 557 IF ( jvar == 1 ) THEN 558 cout1=TRIM(surfdata%cvars(1)) 559 ELSE 560 WRITE(cout1,'(A,A1,A)') TRIM(cout1), '/', TRIM(surfdata%cvars(jvar)) 561 ENDIF 562 END DO 563 486 564 WRITE(numout,*) 487 WRITE(numout,'(1X,A)')TRIM( surfdata%cvars(1))//' data'565 WRITE(numout,'(1X,A)')TRIM( cout1 )//' data' 488 566 WRITE(numout,'(1X,A)')'--------------' 489 567 DO jj = 1,8 … … 495 573 & '---------------------------------------------------------------' 496 574 WRITE(numout,'(1X,A,I8)') & 497 & 'Total data for variable '//TRIM( surfdata%cvars(1))// &575 & 'Total data for variable '//TRIM( cout1 )// & 498 576 & ' = ', iobsmpp 499 577 WRITE(numout,'(1X,A)') & … … 506 584 ! Deallocate temporary data 507 585 !----------------------------------------------------------------------- 508 DEALLOCATE( ifileidx, isurfidx, zdat, clvars )586 DEALLOCATE( ifileidx, isurfidx, zdat, clvarsin ) 509 587 510 588 !----------------------------------------------------------------------- -
branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_rot_vel.F90
r11203 r15764 6 6 7 7 !!---------------------------------------------------------------------- 8 !! obs_rotvel : Rotate velocity data into N-S,E-W directorions 8 !! obs_rotvel_pro : Rotate profile velocity data into N-S,E-W directions 9 !! obs_rotvel_surf : Rotate surface velocity data into N-S,E-W directions 9 10 !!---------------------------------------------------------------------- 10 11 !! * Modules used … … 17 18 USE obs_utils ! For error handling 18 19 USE obs_profiles_def ! Profile definitions 20 USE obs_surf_def ! Surface definitions 19 21 USE obs_inter_h2d ! Horizontal interpolation 20 22 USE obs_inter_sup ! MPP support routines for interpolation … … 27 29 PRIVATE 28 30 29 PUBLIC obs_rotvel ! Rotate the observations 31 PUBLIC obs_rotvel_pro, & ! Rotate the profile velocity observations 32 & obs_rotvel_surf ! Rotate the surface velocity observations 30 33 31 34 !!---------------------------------------------------------------------- … … 37 40 CONTAINS 38 41 39 SUBROUTINE obs_rotvel ( profdata, k2dint, pu, pv )42 SUBROUTINE obs_rotvel_pro( profdata, k2dint, pu, pv ) 40 43 !!--------------------------------------------------------------------- 41 44 !! … … 228 231 CALL wrk_dealloc(jpi,jpj,zsingu,zcosgu,zsingv,zcosgv) 229 232 230 END SUBROUTINE obs_rotvel 233 END SUBROUTINE obs_rotvel_pro 234 235 SUBROUTINE obs_rotvel_surf( surfdata, k2dint, pu, pv ) 236 !!--------------------------------------------------------------------- 237 !! 238 !! *** ROUTINE obs_rotvel_surf *** 239 !! 240 !! ** Purpose : Rotate surface velocity data into N-S,E-W directorions 241 !! 242 !! ** Method : Interpolation of geo2ocean coefficients on U,V grid 243 !! to observation point followed by a similar computations 244 !! as in geo2ocean. 245 !! 246 !! ** Action : Review if there is a better way to do this. 247 !! 248 !! References : 249 !! 250 !! History : 251 !! ! : 2009-02 (K. Mogensen) : New routine 252 !!---------------------------------------------------------------------- 253 !! * Modules used 254 !! * Arguments 255 TYPE(obs_surf), INTENT(INOUT) :: surfdata ! Surface data to be read 256 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation methed 257 REAL(wp), DIMENSION(*) :: & 258 & pu, & 259 & pv 260 !! * Local declarations 261 REAL(wp), DIMENSION(2,2,1) :: zweig 262 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 263 & zmasku, & 264 & zmaskv, & 265 & zcoslu, & 266 & zsinlu, & 267 & zcoslv, & 268 & zsinlv, & 269 & zglamu, & 270 & zgphiu, & 271 & zglamv, & 272 & zgphiv 273 REAL(wp), DIMENSION(1) :: & 274 & zsinu, & 275 & zcosu, & 276 & zsinv, & 277 & zcosv 278 REAL(wp) :: zsin 279 REAL(wp) :: zcos 280 REAL(wp), DIMENSION(1) :: zobsmask 281 REAL(wp), POINTER, DIMENSION(:,:) :: zsingu,zcosgu,zsingv,zcosgv 282 INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 283 & igrdiu, & 284 & igrdju, & 285 & igrdiv, & 286 & igrdjv 287 INTEGER :: ji 288 INTEGER :: jk 289 290 CALL wrk_alloc(jpi,jpj,zsingu,zcosgu,zsingv,zcosgv) 291 292 !----------------------------------------------------------------------- 293 ! Allocate data for message parsing and interpolation 294 !----------------------------------------------------------------------- 295 296 ALLOCATE( & 297 & igrdiu(2,2,surfdata%nsurf), & 298 & igrdju(2,2,surfdata%nsurf), & 299 & zglamu(2,2,surfdata%nsurf), & 300 & zgphiu(2,2,surfdata%nsurf), & 301 & zmasku(2,2,surfdata%nsurf), & 302 & zcoslu(2,2,surfdata%nsurf), & 303 & zsinlu(2,2,surfdata%nsurf), & 304 & igrdiv(2,2,surfdata%nsurf), & 305 & igrdjv(2,2,surfdata%nsurf), & 306 & zglamv(2,2,surfdata%nsurf), & 307 & zgphiv(2,2,surfdata%nsurf), & 308 & zmaskv(2,2,surfdata%nsurf), & 309 & zcoslv(2,2,surfdata%nsurf), & 310 & zsinlv(2,2,surfdata%nsurf) & 311 & ) 312 313 !----------------------------------------------------------------------- 314 ! Receive the angles on the U and V grids. 315 !----------------------------------------------------------------------- 316 317 CALL obs_rot( zsingu, zcosgu, zsingv, zcosgv ) 318 319 DO ji = 1, surfdata%nsurf 320 igrdiu(1,1,ji) = surfdata%mi(ji,1)-1 321 igrdju(1,1,ji) = surfdata%mj(ji,1)-1 322 igrdiu(1,2,ji) = surfdata%mi(ji,1)-1 323 igrdju(1,2,ji) = surfdata%mj(ji,1) 324 igrdiu(2,1,ji) = surfdata%mi(ji,1) 325 igrdju(2,1,ji) = surfdata%mj(ji,1)-1 326 igrdiu(2,2,ji) = surfdata%mi(ji,1) 327 igrdju(2,2,ji) = surfdata%mj(ji,1) 328 igrdiv(1,1,ji) = surfdata%mi(ji,2)-1 329 igrdjv(1,1,ji) = surfdata%mj(ji,2)-1 330 igrdiv(1,2,ji) = surfdata%mi(ji,2)-1 331 igrdjv(1,2,ji) = surfdata%mj(ji,2) 332 igrdiv(2,1,ji) = surfdata%mi(ji,2) 333 igrdjv(2,1,ji) = surfdata%mj(ji,2)-1 334 igrdiv(2,2,ji) = surfdata%mi(ji,2) 335 igrdjv(2,2,ji) = surfdata%mj(ji,2) 336 END DO 337 338 CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, & 339 & glamu, zglamu ) 340 CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, & 341 & gphiu, zgphiu ) 342 CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, & 343 & umask(:,:,1), zmasku ) 344 CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, & 345 & zsingu, zsinlu ) 346 CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, & 347 & zcosgu, zcoslu ) 348 CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, & 349 & glamv, zglamv ) 350 CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, & 351 & gphiv, zgphiv ) 352 CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, & 353 & vmask(:,:,1), zmaskv ) 354 CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, & 355 & zsingv, zsinlv ) 356 CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, & 357 & zcosgv, zcoslv ) 358 359 DO ji = 1, surfdata%nsurf 360 361 CALL obs_int_h2d_init( 1, 1, k2dint, & 362 & surfdata%rlam(ji), surfdata%rphi(ji), & 363 & zglamu(:,:,ji), zgphiu(:,:,ji), & 364 & zmasku(:,:,ji), zweig, zobsmask ) 365 366 CALL obs_int_h2d( 1, 1, zweig, zsinlu(:,:,ji), zsinu ) 367 368 CALL obs_int_h2d( 1, 1, zweig, zcoslu(:,:,ji), zcosu ) 369 370 CALL obs_int_h2d_init( 1, 1, k2dint, & 371 & surfdata%rlam(ji), surfdata%rphi(ji), & 372 & zglamv(:,:,ji), zgphiv(:,:,ji), & 373 & zmaskv(:,:,ji), zweig, zobsmask ) 374 375 CALL obs_int_h2d( 1, 1, zweig, zsinlv(:,:,ji), zsinv ) 376 377 CALL obs_int_h2d( 1, 1, zweig, zcoslv(:,:,ji), zcosv ) 378 379 ! Assume that the angle at observation point is the 380 ! mean of u and v cosines/sines 381 382 zcos = 0.5_wp * ( zcosu(1) + zcosv(1) ) 383 zsin = 0.5_wp * ( zsinu(1) + zsinv(1) ) 384 385 IF ( ( surfdata%rmod(ji,1) /= fbrmdi ) .AND. & 386 & ( surfdata%rmod(ji,2) /= fbrmdi ) ) THEN 387 pu(ji) = surfdata%rmod(ji,1) * zcos - & 388 & surfdata%rmod(ji,2) * zsin 389 pv(ji) = surfdata%rmod(ji,2) * zcos + & 390 & surfdata%rmod(ji,1) * zsin 391 ELSE 392 pu(ji) = fbrmdi 393 pv(ji) = fbrmdi 394 ENDIF 395 396 END DO 397 398 DEALLOCATE( & 399 & igrdiu, & 400 & igrdju, & 401 & zglamu, & 402 & zgphiu, & 403 & zmasku, & 404 & zcoslu, & 405 & zsinlu, & 406 & igrdiv, & 407 & igrdjv, & 408 & zglamv, & 409 & zgphiv, & 410 & zmaskv, & 411 & zcoslv, & 412 & zsinlv & 413 & ) 414 415 CALL wrk_dealloc(jpi,jpj,zsingu,zcosgu,zsingv,zcosgv) 416 417 END SUBROUTINE obs_rotvel_surf 231 418 232 419 END MODULE obs_rot_vel -
branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90
r11203 r15764 52 52 INTEGER :: nrec !: Number of surface observation records in window 53 53 54 LOGICAL :: lclim !: Climatology will be calculated for this structure 55 54 56 ! Arrays with size equal to the number of surface observations 55 57 56 58 INTEGER, POINTER, DIMENSION(:) :: & 57 & mi, & !: i-th grid coord. for interpolating to surface observation58 & mj, & !: j-th grid coord. for interpolating to surface observation59 59 & mt, & !: time record number for gridded data 60 60 & nsidx,& !: Surface observation number … … 69 69 & ntyp !: Type of surface observation product 70 70 71 INTEGER, POINTER, DIMENSION(:,:) :: & 72 & mi, & !: i-th grid coord. for interpolating to surface observation 73 & mj !: j-th grid coord. for interpolating to surface observation 74 71 75 CHARACTER(len=8), POINTER, DIMENSION(:) :: & 72 76 & cvars !: Variable names … … 84 88 REAL(KIND=wp), POINTER, DIMENSION(:,:) :: & 85 89 & robs, & !: Surface observation 86 & rmod !: Model counterpart of the surface observation vector 90 & rmod, & !: Model counterpart of the surface observation vector 91 & rclm !: Climatological counterpart of the surface observation vector 87 92 88 93 REAL(KIND=wp), POINTER, DIMENSION(:,:) :: & 89 94 & rext !: Extra fields interpolated to observation points 90 95 91 REAL(KIND=wp), POINTER, DIMENSION(:,: ) :: &96 REAL(KIND=wp), POINTER, DIMENSION(:,:,:) :: & 92 97 & vdmean !: Time averaged of model field 93 98 … … 111 116 112 117 ! Is this a gridded product? 113 114 118 LOGICAL :: lgrid 115 119 … … 124 128 CONTAINS 125 129 126 SUBROUTINE obs_surf_alloc( surf, ksurf, kvar, kextra, kstp, kpi, kpj )130 SUBROUTINE obs_surf_alloc( surf, ksurf, kvar, kextra, kstp, kpi, kpj, ldclim ) 127 131 !!---------------------------------------------------------------------- 128 132 !! *** ROUTINE obs_surf_alloc *** … … 143 147 INTEGER, INTENT(IN) :: kpi ! Number of 3D grid points 144 148 INTEGER, INTENT(IN) :: kpj 149 LOGICAL, INTENT(IN) :: ldclim 145 150 146 151 !!* Local variables … … 157 162 surf%npi = kpi 158 163 surf%npj = kpj 164 surf%lclim = ldclim 159 165 160 166 ! Allocate arrays of size number of variables … … 171 177 172 178 ALLOCATE( & 173 & surf%mi(ksurf), &174 & surf%mj(ksurf), &175 179 & surf%mt(ksurf), & 176 180 & surf%nsidx(ksurf), & … … 196 200 197 201 ALLOCATE( & 202 & surf%mi(ksurf,kvar), & 203 & surf%mj(ksurf,kvar), & 198 204 & surf%robs(ksurf,kvar), & 199 205 & surf%rmod(ksurf,kvar) & 200 206 & ) 201 207 208 IF (surf%lclim) ALLOCATE( surf%rclm(ksurf,kvar) ) 209 202 210 ! Allocate arrays of number of extra fields at observation points 203 211 … … 223 231 224 232 ALLOCATE( & 225 & surf%vdmean(kpi,kpj ) &233 & surf%vdmean(kpi,kpj,kvar) & 226 234 & ) 227 235 … … 293 301 & ) 294 302 303 IF (surf%lclim) DEALLOCATE( surf%rclm ) 295 304 ! Deallocate arrays of number of extra fields at observation points 296 305 … … 371 380 IF ( lallocate ) THEN 372 381 CALL obs_surf_alloc( newsurf, insurf, surf%nvar, & 373 & surf%nextra, surf%nstp, surf%npi, surf%npj )382 & surf%nextra, surf%nstp, surf%npi, surf%npj, surf%lclim ) 374 383 ENDIF 375 384 … … 418 427 newsurf%robs(insurf,jk) = surf%robs(ji,jk) 419 428 newsurf%rmod(insurf,jk) = surf%rmod(ji,jk) 429 IF (newsurf%lclim) newsurf%rclm(insurf,jk) = surf%rclm(ji,jk) 420 430 421 431 END DO … … 514 524 oldsurf%robs(jj,jk) = surf%robs(ji,jk) 515 525 oldsurf%rmod(jj,jk) = surf%rmod(ji,jk) 526 IF (surf%lclim) oldsurf%rclm(jj,jk) = surf%rclm(ji,jk) 516 527 517 528 END DO -
branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90
r11203 r15764 28 28 USE obs_const 29 29 USE obs_mpp ! MPP support routines for observation diagnostics 30 USE lib_mpp 30 USE lib_mpp ! MPP routines 31 31 32 32 IMPLICIT NONE … … 95 95 INTEGER :: je 96 96 INTEGER :: iadd 97 INTEGER :: iadd_clm ! 1 if climatology present 97 98 INTEGER :: iext 98 99 REAL(wp) :: zpres 100 101 iadd_clm = 0 102 IF ( profdata%lclim ) iadd_clm = 1 99 103 100 104 IF ( PRESENT( padd ) ) THEN … … 137 141 fbdata%caddunit(1,1) = 'Degrees centigrade' 138 142 fbdata%caddunit(1,2) = 'PSU' 143 IF ( profdata%lclim ) THEN 144 fbdata%caddlong(2,1) = 'Climatology interpolated potential temperature' 145 fbdata%caddlong(2,2) = 'Climatology interpolated practical salinity' 146 fbdata%caddunit(2,1) = 'Degrees centigrade' 147 fbdata%caddunit(2,2) = 'PSU' 148 ENDIF 139 149 fbdata%cgrid(:) = 'T' 140 150 DO je = 1, iext … … 170 180 fbdata%caddunit(1,1) = 'm/s' 171 181 fbdata%caddunit(1,2) = 'm/s' 182 IF ( profdata%lclim ) THEN 183 fbdata%caddlong(2,1) = 'Climatology interpolated zonal velocity' 184 fbdata%caddlong(2,2) = 'Climatology interpolated meridional velocity' 185 fbdata%caddunit(2,1) = 'm/s' 186 fbdata%caddunit(2,2) = 'm/s' 187 ENDIF 172 188 fbdata%cgrid(1) = 'U' 173 189 fbdata%cgrid(2) = 'V' … … 252 268 fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(cllongname) 253 269 fbdata%caddunit(1,1) = clunits 270 IF ( profdata%lclim ) THEN 271 fbdata%caddlong(2,1) = 'Climatological interpolated ' // TRIM(cllongname) 272 fbdata%caddunit(2,1) = clunits 273 ENDIF 254 274 fbdata%cgrid(:) = clgrid 255 275 DO je = 1, iext … … 266 286 267 287 fbdata%caddname(1) = 'Hx' 288 289 IF ( profdata%lclim ) fbdata%caddname(1+iadd_clm) = 'CLM' 268 290 269 291 WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc … … 319 341 DO jk = profdata%npvsta(jo,jvar), profdata%npvend(jo,jvar) 320 342 ik = profdata%var(jvar)%nvlidx(jk) 321 fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk)322 343 fbdata%pob(ik,jo,jvar) = profdata%var(jvar)%vobs(jk) 323 344 fbdata%pdep(ik,jo) = profdata%var(jvar)%vdep(jk) … … 327 348 fbdata%ivlqc(ik,jo,jvar) = IBSET(profdata%var(jvar)%nvqc(jk),2) 328 349 fbdata%ivlqcf(1,ik,jo,jvar) = profdata%var(jvar)%nvqcf(1,jk) 329 fbdata%ivlqcf(2,ik,jo,jvar) = IAND(profdata%var(jvar)%nvqc(jk),b'0000 0000 11111111')350 fbdata%ivlqcf(2,ik,jo,jvar) = IAND(profdata%var(jvar)%nvqc(jk),b'0000000011111111') 330 351 ELSE 331 352 fbdata%ivlqc(ik,jo,jvar) = profdata%var(jvar)%nvqc(jk) … … 333 354 ENDIF 334 355 fbdata%iobsk(ik,jo,jvar) = profdata%var(jvar)%mvk(jk) 335 DO ja = 1, iadd 336 fbdata%padd(ik,jo,1+ja,jvar) = & 356 fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk) 357 IF ( profdata%lclim ) THEN 358 fbdata%padd(ik,jo,1+iadd_clm,jvar) = profdata%var(jvar)%vclm(jk) 359 ENDIF 360 DO ja = 1, iadd 361 fbdata%padd(ik,jo,1+iadd_clm+ja,jvar) = & 337 362 & profdata%var(jvar)%vext(jk,padd%ipoint(ja)) 338 363 END DO … … 416 441 INTEGER :: ja 417 442 INTEGER :: je 443 INTEGER :: jv 418 444 INTEGER :: iadd 419 445 INTEGER :: iext 420 446 INTEGER :: indx_std 421 447 INTEGER :: iadd_std 422 448 INTEGER :: iadd_clm 449 INTEGER :: iadd_mdt 450 451 IF ( PRESENT( pext ) ) THEN 452 iext = pext%inum 453 ELSE 454 iext = 0 455 ENDIF 456 457 ! Set up number of additional variables to be ouput: 458 ! Hx, CLM, STD, MDT... 423 459 IF ( PRESENT( padd ) ) THEN 424 460 iadd = padd%inum 425 461 ELSE 426 462 iadd = 0 427 ENDIF428 429 IF ( PRESENT( pext ) ) THEN430 iext = pext%inum431 ELSE432 iext = 0433 463 ENDIF 434 464 … … 444 474 ENDIF 445 475 476 iadd_clm = 0 477 IF ( surfdata%lclim ) iadd_clm = 1 478 479 iadd_mdt = 0 480 IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) iadd_mdt = 1 481 446 482 CALL init_obfbdata( fbdata ) 447 483 … … 451 487 ! SLA needs special treatment because of MDT, so is all done here 452 488 ! Other variables are done more generically 453 454 CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 455 & 2 + iadd + iadd_std, 1 + iext, .TRUE. ) 489 ! No climatology for SLA, MDT is our best estimate of that and is already output. 490 491 CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 492 & 1 + iadd_mdt + iadd_std + iadd, & 493 & 1 + iext, .TRUE. ) 456 494 457 495 clfiletype = 'slafb' … … 493 531 clgrid = 'T' 494 532 533 CASE('SIT') 534 535 clfiletype = 'sitfb' 536 cllongname = 'Sea ice thickness' 537 clunits = 'm' 538 clgrid = 'T' 539 540 CASE('FBD') 541 542 clfiletype = 'fbdfb' 543 ! Change label from FBD ("freeboard") to SIT ("Sea Ice Thickness") 544 surfdata%cvars(1) = 'SIT' 545 cllongname = 'Sea ice thickness' 546 clunits = 'm' 547 clgrid = 'T' 548 495 549 CASE('SSS') 496 550 … … 500 554 clgrid = 'T' 501 555 502 CASE('SLCHLTOT','LOGCHL','LogChl','logchl') 503 504 clfiletype = 'slchltotfb' 505 cllongname = 'Surface total log10(chlorophyll)' 506 clunits = 'log10(mg/m3)' 507 clgrid = 'T' 508 509 CASE('SLCHLDIA') 510 511 clfiletype = 'slchldiafb' 512 cllongname = 'Surface diatom log10(chlorophyll)' 513 clunits = 'log10(mg/m3)' 514 clgrid = 'T' 515 516 CASE('SLCHLNON') 517 518 clfiletype = 'slchlnonfb' 519 cllongname = 'Surface non-diatom log10(chlorophyll)' 520 clunits = 'log10(mg/m3)' 521 clgrid = 'T' 522 523 CASE('SLCHLDIN') 524 525 clfiletype = 'slchldinfb' 526 cllongname = 'Surface dinoflagellate log10(chlorophyll)' 527 clunits = 'log10(mg/m3)' 528 clgrid = 'T' 529 530 CASE('SLCHLMIC') 531 532 clfiletype = 'slchlmicfb' 533 cllongname = 'Surface microphytoplankton log10(chlorophyll)' 534 clunits = 'log10(mg/m3)' 535 clgrid = 'T' 536 537 CASE('SLCHLNAN') 538 539 clfiletype = 'slchlnanfb' 540 cllongname = 'Surface nanophytoplankton log10(chlorophyll)' 541 clunits = 'log10(mg/m3)' 542 clgrid = 'T' 543 544 CASE('SLCHLPIC') 545 546 clfiletype = 'slchlpicfb' 547 cllongname = 'Surface picophytoplankton log10(chlorophyll)' 548 clunits = 'log10(mg/m3)' 549 clgrid = 'T' 550 551 CASE('SCHLTOT') 552 553 clfiletype = 'schltotfb' 554 cllongname = 'Surface total chlorophyll' 555 clunits = 'mg/m3' 556 clgrid = 'T' 557 558 CASE('SLPHYTOT') 559 560 clfiletype = 'slphytotfb' 561 cllongname = 'Surface total log10(phytoplankton carbon)' 562 clunits = 'log10(mmolC/m3)' 563 clgrid = 'T' 564 565 CASE('SLPHYDIA') 566 567 clfiletype = 'slphydiafb' 568 cllongname = 'Surface diatom log10(phytoplankton carbon)' 569 clunits = 'log10(mmolC/m3)' 570 clgrid = 'T' 571 572 CASE('SLPHYNON') 573 574 clfiletype = 'slphynonfb' 575 cllongname = 'Surface non-diatom log10(phytoplankton carbon)' 576 clunits = 'log10(mmolC/m3)' 577 clgrid = 'T' 578 579 CASE('SSPM') 580 581 clfiletype = 'sspmfb' 582 cllongname = 'Surface suspended particulate matter' 583 clunits = 'g/m3' 584 clgrid = 'T' 585 586 CASE('SFCO2','FCO2','fCO2','fco2') 587 588 clfiletype = 'sfco2fb' 589 cllongname = 'Surface fugacity of carbon dioxide' 590 clunits = 'uatm' 591 clgrid = 'T' 592 593 CASE('SPCO2') 594 595 clfiletype = 'spco2fb' 596 cllongname = 'Surface partial pressure of carbon dioxide' 597 clunits = 'uatm' 598 clgrid = 'T' 556 CASE('UVEL') 557 558 clfiletype = 'ssvfb' 559 cllongname = 'Eastward sea surface velocity' 560 clunits = 'm s-1' 561 clgrid = 'U' 562 cllongname = 'Northward sea surface velocity' 563 clunits = 'm s-1' 564 clgrid = 'V' 565 566 CASE('SLCHLTOT') 567 568 clfiletype = 'slchltotfb' 569 cllongname = 'Surface total log10(chlorophyll)' 570 clunits = 'log10(mg/m3)' 571 clgrid = 'T' 572 573 CASE('SLCHLDIA') 574 575 clfiletype = 'slchldiafb' 576 cllongname = 'Surface diatom log10(chlorophyll)' 577 clunits = 'log10(mg/m3)' 578 clgrid = 'T' 579 580 CASE('SLCHLNON') 581 582 clfiletype = 'slchlnonfb' 583 cllongname = 'Surface non-diatom log10(chlorophyll)' 584 clunits = 'log10(mg/m3)' 585 clgrid = 'T' 586 587 CASE('SLCHLDIN') 588 589 clfiletype = 'slchldinfb' 590 cllongname = 'Surface dinoflagellate log10(chlorophyll)' 591 clunits = 'log10(mg/m3)' 592 clgrid = 'T' 593 594 CASE('SLCHLMIC') 595 596 clfiletype = 'slchlmicfb' 597 cllongname = 'Surface microphytoplankton log10(chlorophyll)' 598 clunits = 'log10(mg/m3)' 599 clgrid = 'T' 600 601 CASE('SLCHLNAN') 602 603 clfiletype = 'slchlnanfb' 604 cllongname = 'Surface nanophytoplankton log10(chlorophyll)' 605 clunits = 'log10(mg/m3)' 606 clgrid = 'T' 607 608 CASE('SLCHLPIC') 609 610 clfiletype = 'slchlpicfb' 611 cllongname = 'Surface picophytoplankton log10(chlorophyll)' 612 clunits = 'log10(mg/m3)' 613 clgrid = 'T' 614 615 CASE('SCHLTOT') 616 617 clfiletype = 'schltotfb' 618 cllongname = 'Surface total chlorophyll' 619 clunits = 'mg/m3' 620 clgrid = 'T' 621 622 CASE('SLPHYTOT') 623 624 clfiletype = 'slphytotfb' 625 cllongname = 'Surface total log10(phytoplankton carbon)' 626 clunits = 'log10(mmolC/m3)' 627 clgrid = 'T' 628 629 CASE('SLPHYDIA') 630 631 clfiletype = 'slphydiafb' 632 cllongname = 'Surface diatom log10(phytoplankton carbon)' 633 clunits = 'log10(mmolC/m3)' 634 clgrid = 'T' 635 636 CASE('SLPHYNON') 637 638 clfiletype = 'slphynonfb' 639 cllongname = 'Surface non-diatom log10(phytoplankton carbon)' 640 clunits = 'log10(mmolC/m3)' 641 clgrid = 'T' 642 643 CASE('SSPM') 644 645 clfiletype = 'sspmfb' 646 cllongname = 'Surface suspended particulate matter' 647 clunits = 'g/m3' 648 clgrid = 'T' 649 650 CASE('SKD490') 651 652 clfiletype = 'skd490fb' 653 cllongname = 'Surface attenuation coefficient of downwelling radiation at 490 nm' 654 clunits = 'm-1' 655 clgrid = 'T' 656 657 CASE('SFCO2') 658 659 clfiletype = 'sfco2fb' 660 cllongname = 'Surface fugacity of carbon dioxide' 661 clunits = 'uatm' 662 clgrid = 'T' 663 664 CASE('SPCO2') 665 666 clfiletype = 'spco2fb' 667 cllongname = 'Surface partial pressure of carbon dioxide' 668 clunits = 'uatm' 669 clgrid = 'T' 599 670 600 671 CASE DEFAULT … … 607 678 ! Remaining variables treated more generically 608 679 609 IF ( TRIM(surfdata%cvars(1)) /= 'SLA' ) THEN 610 611 CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 612 & 1 + iadd + iadd_std, iext, .TRUE. ) 613 614 fbdata%cname(1) = surfdata%cvars(1) 615 fbdata%coblong(1) = cllongname 616 fbdata%cobunit(1) = clunits 617 DO je = 1, iext 618 fbdata%cextname(je) = pext%cdname(je) 619 fbdata%cextlong(je) = pext%cdlong(je,1) 620 fbdata%cextunit(je) = pext%cdunit(je,1) 621 END DO 622 IF ( TRIM(surfdata%cvars(1)) == 'ICECONC' ) THEN 623 fbdata%caddlong(1,1) = 'Model interpolated ICE' 624 ELSE 625 fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(surfdata%cvars(1)) 626 ENDIF 627 fbdata%caddunit(1,1) = clunits 628 fbdata%cgrid(1) = clgrid 629 DO ja = 1, iadd 630 fbdata%caddname(1+iadd_std+ja) = padd%cdname(ja) 631 fbdata%caddlong(1+iadd_std+ja,1) = padd%cdlong(ja,1) 632 fbdata%caddunit(1+iadd_std+ja,1) = padd%cdunit(ja,1) 633 END DO 634 680 IF ( TRIM(surfdata%cvars(1)) /= 'SLA' ) THEN 681 682 CALL alloc_obfbdata( fbdata, surfdata%nvar, surfdata%nsurf, 1, & 683 & 1 + iadd_std + iadd_clm + iadd, iext, .TRUE. ) 684 685 DO jv = 1, surfdata%nvar 686 fbdata%cname(jv) = surfdata%cvars(jv) 687 fbdata%coblong(jv) = cllongname 688 fbdata%cobunit(jv) = clunits 689 END DO 690 DO je = 1, iext 691 fbdata%cextname(je) = pext%cdname(je) 692 fbdata%cextlong(je) = pext%cdlong(je,1) 693 fbdata%cextunit(je) = pext%cdunit(je,1) 694 END DO 695 DO jv = 1, surfdata%nvar 696 IF ( TRIM(surfdata%cvars(1)) == 'ICECONC' ) THEN 697 fbdata%caddlong(1,jv) = 'Model interpolated ICE' 698 ELSE 699 fbdata%caddlong(1,jv) = 'Model interpolated ' // TRIM(surfdata%cvars(jv)) 700 ENDIF 701 fbdata%caddunit(1,jv) = clunits 702 fbdata%cgrid(jv) = clgrid 703 END DO 704 DO ja = 1, iadd 705 fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm+ja) = padd%cdname(ja) 706 DO jv = 1, surfdata%nvar 707 fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = padd%cdlong(ja,jv) 708 fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = padd%cdunit(ja,jv) 709 END DO 710 END DO 635 711 ENDIF 636 712 637 713 fbdata%caddname(1) = 'Hx' 638 IF ( indx_std /= -1 ) THEN 639 IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) iadd_std = iadd_std + 1 640 fbdata%caddname(1+iadd_std) = surfdata%cext(indx_std) 641 fbdata%caddlong(1+iadd_std,1) = 'Obs error standard deviation' 642 fbdata%caddunit(1+iadd_std,1) = fbdata%cobunit(1) 714 IF ( indx_std /= -1 ) THEN 715 fbdata%caddname(1+iadd_mdt+iadd_std) = surfdata%cext(indx_std) 716 DO jv = 1, surfdata%nvar 717 fbdata%caddlong(1+iadd_mdt+iadd_std,1) = 'Obs error standard deviation' 718 fbdata%caddunit(1+iadd_mdt+iadd_std,1) = fbdata%cobunit(1) 719 END DO 720 ENDIF 721 722 IF ( surfdata%lclim ) THEN 723 fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm) = 'CLM' 724 DO jv = 1, surfdata%nvar 725 fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm,jv) = 'Climatology' 726 fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm,jv) = fbdata%cobunit(1) 727 END DO 643 728 ENDIF 644 729 … … 649 734 WRITE(numout,*)'obs_wri_surf :' 650 735 WRITE(numout,*)'~~~~~~~~~~~~~' 651 WRITE(numout,*)'Writing '//TRIM( surfdata%cvars(1))//' feedback file : ',TRIM(clfname)736 WRITE(numout,*)'Writing '//TRIM(clfiletype)//' feedback file : ',TRIM(clfname) 652 737 ENDIF 653 738 … … 663 748 fbdata%ioqc(jo) = 4 664 749 fbdata%ioqcf(1,jo) = 0 665 fbdata%ioqcf(2,jo) = IAND(surfdata%nqc(jo),b'0000 0000 11111111')750 fbdata%ioqcf(2,jo) = IAND(surfdata%nqc(jo),b'0000000011111111') 666 751 ELSE 667 752 fbdata%ioqc(jo) = surfdata%nqc(jo) … … 675 760 fbdata%kindex(jo) = surfdata%nsfil(jo) 676 761 IF (ln_grid_global) THEN 677 fbdata%iobsi(jo,1) = surfdata%mi(jo) 678 fbdata%iobsj(jo,1) = surfdata%mj(jo) 762 DO jv = 1, surfdata%nvar 763 fbdata%iobsi(jo,jv) = surfdata%mi(jo,jv) 764 fbdata%iobsj(jo,jv) = surfdata%mj(jo,jv) 765 END DO 679 766 ELSE 680 fbdata%iobsi(jo,1) = mig(surfdata%mi(jo)) 681 fbdata%iobsj(jo,1) = mjg(surfdata%mj(jo)) 767 DO jv = 1, surfdata%nvar 768 fbdata%iobsi(jo,jv) = mig(surfdata%mi(jo,jv)) 769 fbdata%iobsj(jo,jv) = mjg(surfdata%mj(jo,jv)) 770 END DO 682 771 ENDIF 683 772 CALL greg2jul( 0, & … … 689 778 & fbdata%ptim(jo), & 690 779 & krefdate = 19500101 ) 691 fbdata%padd(1,jo,1,1) = surfdata%rmod(jo,1)692 IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%padd(1,jo,2,1) = surfdata%rext(jo,1)693 fbdata%pob(1,jo,1) = surfdata%robs(jo,1)694 780 fbdata%pdep(1,jo) = 0.0 695 781 fbdata%idqc(1,jo) = 0 696 782 fbdata%idqcf(:,1,jo) = 0 697 IF ( surfdata%nqc(jo) > 255 ) THEN 698 fbdata%ivqc(jo,1) = 4 699 fbdata%ivlqc(1,jo,1) = 4 700 fbdata%ivlqcf(1,1,jo,1) = 0 701 fbdata%ivlqcf(2,1,jo,1) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111') 702 ELSE 703 fbdata%ivqc(jo,1) = surfdata%nqc(jo) 704 fbdata%ivlqc(1,jo,1) = surfdata%nqc(jo) 705 fbdata%ivlqcf(:,1,jo,1) = 0 706 ENDIF 707 fbdata%iobsk(1,jo,1) = 0 708 IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2) 709 IF ( indx_std /= -1 ) THEN 710 fbdata%padd(1,jo,1+iadd_std,1) = surfdata%rext(jo,indx_std) 711 ENDIF 712 713 DO ja = 1, iadd 714 fbdata%padd(1,jo,2+iadd_std+ja,1) = & 715 & surfdata%rext(jo,padd%ipoint(ja)) 716 END DO 717 DO je = 1, iext 718 fbdata%pext(1,jo,1+je) = & 719 & surfdata%rext(jo,pext%ipoint(je)) 783 DO jv = 1, surfdata%nvar 784 fbdata%pob(1,jo,jv) = surfdata%robs(jo,jv) 785 786 IF ( surfdata%nqc(jo) > 255 ) THEN 787 fbdata%ivqc(jo,jv) = 4 788 fbdata%ivlqc(1,jo,jv) = 4 789 fbdata%ivlqcf(1,1,jo,jv) = 0 790 fbdata%ivlqcf(2,1,jo,jv) = IAND(surfdata%nqc(jo),b'0000000011111111') 791 ELSE 792 fbdata%ivqc(jo,jv) = surfdata%nqc(jo) 793 fbdata%ivlqc(1,jo,jv) = surfdata%nqc(jo) 794 fbdata%ivlqcf(:,1,jo,jv) = 0 795 ENDIF 796 fbdata%iobsk(1,jo,jv) = 0 797 798 ! Additional variables. 799 ! Hx is always the first additional variable 800 fbdata%padd(1,jo,1,jv) = surfdata%rmod(jo,jv) 801 ! MDT is output as an additional variable if SLA obs type 802 IF ( TRIM(surfdata%cvars(jv)) == 'SLA' ) THEN 803 fbdata%padd(1,jo,1+iadd_mdt,jv) = surfdata%rext(jo,1) 804 ENDIF 805 ! STD is output as an additional variable if available 806 IF ( indx_std /= -1 ) THEN 807 fbdata%padd(1,jo,1+iadd_mdt+iadd_std,jv) = surfdata%rext(jo,indx_std) 808 ENDIF 809 ! CLM is output as an additional variable if available 810 IF ( surfdata%lclim ) THEN 811 fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm,jv) = surfdata%rclm(jo,jv) 812 ENDIF 813 ! Then other additional variables are output 814 DO ja = 1, iadd 815 fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = & 816 & surfdata%rext(jo,padd%ipoint(ja)) 817 END DO 818 END DO 819 ! Extra variables 820 IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2) 821 DO je = 1, iext 822 fbdata%pext(1,jo,1+je) = & 823 & surfdata%rext(jo,pext%ipoint(je)) 720 824 END DO 721 825 END DO -
branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_z1d.h90
r11203 r15764 62 62 z1dm = ( pdep(kkco(jdep)) - pobsdep(jdep) ) 63 63 z1dp = ( pobsdep(jdep) - pdep(kkco(jdep)-1) ) 64 IF ( pobsmask(kkco(jdep)) == 0.0_wp ) z1dp = 0.0_wp65 66 zsum = z1dm + z1dp67 64 68 IF ( k1dint == 0 ) THEN 69 70 !----------------------------------------------------------------- 71 ! Linear interpolation 72 !----------------------------------------------------------------- 73 pobs(jdep) = ( z1dm * pobsk(kkco(jdep)-1) & 74 & + z1dp * pobsk(kkco(jdep) ) ) / zsum 75 76 ELSEIF ( k1dint == 1 ) THEN 77 78 !----------------------------------------------------------------- 79 ! Cubic spline interpolation 80 !----------------------------------------------------------------- 81 zsum2 = zsum * zsum 82 pobs(jdep) = ( z1dm * pobsk (kkco(jdep)-1) & 83 & + z1dp * pobsk (kkco(jdep) ) & 84 & + ( z1dm * ( z1dm * z1dm - zsum2 ) * pobs2k(kkco(jdep)-1) & 85 & + z1dp * ( z1dp * z1dp - zsum2 ) * pobs2k(kkco(jdep) ) & 86 & ) / 6.0_wp & 87 & ) / zsum 88 65 ! Where both levels are masked, return a fill value 66 IF ( ( pobsmask(kkco(jdep)-1) == 0.0_wp ) .AND. (pobsmask(kkco(jdep)) == 0.0_wp) ) THEN 67 pobs(jdep) = 99999. 68 ELSE 69 70 ! Where upper level is masked (e.g., under ice cavity), only use deeper level 71 ! otherwise where ob is at or above upper level model T-point, 72 ! use upper model level rather than extrapolate 73 IF ( pobsmask(kkco(jdep)-1) == 0.0_wp ) THEN 74 z1dm = 0.0_wp 75 ELSE IF ( pobsdep(jdep) <= pdep(kkco(jdep)-1) ) THEN 76 z1dp = 0.0_wp 77 END IF 78 79 ! Where deeper level is masked (e.g., near sea bed), only use upper level 80 ! otherwise where ob is at or below deeper level model T-point, 81 ! use deeper model level rather than extrapolate 82 IF ( pobsmask(kkco(jdep)) == 0.0_wp ) THEN 83 z1dp = 0.0_wp 84 ELSE IF ( pobsdep(jdep) >= pdep(kkco(jdep)) ) THEN 85 z1dm = 0.0_wp 86 END IF 87 88 zsum = z1dm + z1dp 89 90 IF ( k1dint == 0 ) THEN 91 92 !----------------------------------------------------------------- 93 ! Linear interpolation 94 !----------------------------------------------------------------- 95 pobs(jdep) = ( z1dm * pobsk(kkco(jdep)-1) & 96 & + z1dp * pobsk(kkco(jdep) ) ) / zsum 97 98 ELSEIF ( k1dint == 1 ) THEN 99 100 !----------------------------------------------------------------- 101 ! Cubic spline interpolation 102 !----------------------------------------------------------------- 103 zsum2 = zsum * zsum 104 pobs(jdep) = ( z1dm * pobsk (kkco(jdep)-1) & 105 & + z1dp * pobsk (kkco(jdep) ) & 106 & + ( z1dm * ( z1dm * z1dm - zsum2 ) * pobs2k(kkco(jdep)-1) & 107 & + z1dp * ( z1dp * z1dp - zsum2 ) * pobs2k(kkco(jdep) ) & 108 & ) / 6.0_wp & 109 & ) / zsum 110 111 ENDIF 89 112 ENDIF 90 113 END DO -
branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
r11203 r15764 56 56 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: resto !: restoring coeff. on T and S (s-1) 57 57 58 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tclim !: temperature climatology on each time step(Celcius) 59 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sclim !: salinity climatology on each time step (psu) 60 58 61 !! * Substitutions 59 62 # include "domzgr_substitute.h90" … … 70 73 !! *** FUNCTION tra_dmp_alloc *** 71 74 !!---------------------------------------------------------------------- 72 ALLOCATE( strdmp(jpi,jpj,jpk) , ttrdmp(jpi,jpj,jpk), resto(jpi,jpj,jpk), STAT= tra_dmp_alloc ) 75 ALLOCATE( strdmp(jpi,jpj,jpk) , ttrdmp(jpi,jpj,jpk), resto(jpi,jpj,jpk), & 76 & tclim(jpi,jpj,jpk) , sclim(jpi,jpj,jpk), STAT=tra_dmp_alloc ) 73 77 ! 74 78 IF( lk_mpp ) CALL mpp_sum ( tra_dmp_alloc ) … … 110 114 ! !== input T-S data at kt ==! 111 115 CALL dta_tsd( kt, zts_dta ) ! read and interpolates T-S data at kt 116 117 tclim(:,:,:) = zts_dta(:,:,:,jp_tem) 118 sclim(:,:,:) = zts_dta(:,:,:,jp_sal) 112 119 ! 113 120 SELECT CASE ( nn_zdmp ) !== type of damping ==!
Note: See TracChangeset
for help on using the changeset viewer.