New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 7713 for branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 – NEMO

Ignore:
Timestamp:
2017-02-22T12:40:19+01:00 (7 years ago)
Author:
dford
Message:

Add observation operator code for surface log10(chlorophyll), SPM, pCO2 and fCO2, for use with FABM-ERSEM, HadOCC and MEDUSA. See internal Met Office NEMO ticket 660.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

    r6301 r7713  
    2323   !!   obs_vel_opt :    Compute the model counterpart of zonal and meridional 
    2424   !!                    components of velocity from observations. 
     25   !!   obs_logchl_opt : Compute the model counterpart of log10(chlorophyll) 
     26   !!                    observations 
     27   !!   obs_spm_opt :    Compute the model counterpart of spm 
     28   !!                    observations 
     29   !!   obs_fco2_opt :   Compute the model counterpart of fco2 
     30   !!                    observations 
     31   !!   obs_pco2_opt :   Compute the model counterpart of pco2 
     32   !!                    observations 
    2533   !!---------------------------------------------------------------------- 
    2634 
     
    6371      &   obs_sss_opt, &  ! Compute the model counterpart of SSS observations 
    6472      &   obs_seaice_opt, & 
    65       &   obs_vel_opt     ! Compute the model counterpart of velocity profile data 
     73      &   obs_vel_opt, &  ! Compute the model counterpart of velocity profile data 
     74      &   obs_logchl_opt, & ! Compute the model counterpart of logchl data 
     75      &   obs_spm_opt, &  ! Compute the model counterpart of spm data 
     76      &   obs_fco2_opt, & ! Compute the model counterpart of fco2 data 
     77      &   obs_pco2_opt    ! Compute the model counterpart of pco2 data 
    6678 
    6779   INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types 
     
    20522064   END SUBROUTINE obs_vel_opt 
    20532065 
     2066   SUBROUTINE obs_logchl_opt( logchldatqc, kt, kpi, kpj, kit000, & 
     2067      &                    plogchln, plogchlmask, k2dint ) 
     2068 
     2069      !!----------------------------------------------------------------------- 
     2070      !! 
     2071      !!                     ***  ROUTINE obs_logchl_opt  *** 
     2072      !! 
     2073      !! ** Purpose : Compute the model counterpart of logchl 
     2074      !!              data by interpolating from the model grid to the  
     2075      !!              observation point. 
     2076      !! 
     2077      !! ** Method  : Linearly interpolate to each observation point using  
     2078      !!              the model values at the corners of the surrounding grid box. 
     2079      !! 
     2080      !!    The now model logchl is first computed at the obs (lon, lat) point. 
     2081      !! 
     2082      !!    Several horizontal interpolation schemes are available: 
     2083      !!        - distance-weighted (great circle) (k2dint = 0) 
     2084      !!        - distance-weighted (small angle)  (k2dint = 1) 
     2085      !!        - bilinear (geographical grid)     (k2dint = 2) 
     2086      !!        - bilinear (quadrilateral grid)    (k2dint = 3) 
     2087      !!        - polynomial (quadrilateral grid)  (k2dint = 4) 
     2088      !! 
     2089      !! 
     2090      !! ** Action  : 
     2091      !! 
     2092      !! History : 
     2093      !!       
     2094      !!----------------------------------------------------------------------- 
     2095 
     2096      !! * Modules used 
     2097      USE obs_surf_def  ! Definition of storage space for surface observations 
     2098 
     2099      IMPLICIT NONE 
     2100 
     2101      !! * Arguments 
     2102      TYPE(obs_surf), INTENT(INOUT) :: logchldatqc     ! Subset of surface data not failing screening 
     2103      INTEGER, INTENT(IN) :: kt       ! Time step 
     2104      INTEGER, INTENT(IN) :: kpi      ! Model grid parameters 
     2105      INTEGER, INTENT(IN) :: kpj 
     2106      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step  
     2107                                      !   (kit000-1 = restart time) 
     2108      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header) 
     2109      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
     2110         & plogchln,  &    ! Model logchl field 
     2111         & plogchlmask     ! Land-sea mask 
     2112          
     2113      !! * Local declarations 
     2114      INTEGER :: ji 
     2115      INTEGER :: jj 
     2116      INTEGER :: jobs 
     2117      INTEGER :: inrc 
     2118      INTEGER :: ilogchl 
     2119      INTEGER :: iobs 
     2120        
     2121      REAL(KIND=wp) :: zlam 
     2122      REAL(KIND=wp) :: zphi 
     2123      REAL(KIND=wp) :: zext(1), zobsmask(1) 
     2124      REAL(kind=wp), DIMENSION(2,2,1) :: & 
     2125         & zweig 
     2126      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
     2127         & zmask, & 
     2128         & zlogchll, & 
     2129         & zglam, & 
     2130         & zgphi 
     2131      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
     2132         & igrdi, & 
     2133         & igrdj 
     2134 
     2135      !------------------------------------------------------------------------ 
     2136      ! Local initialization  
     2137      !------------------------------------------------------------------------ 
     2138      ! ... Record and data counters 
     2139      inrc = kt - kit000 + 2 
     2140      ilogchl = logchldatqc%nsstp(inrc) 
     2141 
     2142      ! Get the data for interpolation 
     2143       
     2144      ALLOCATE( & 
     2145         & igrdi(2,2,ilogchl), & 
     2146         & igrdj(2,2,ilogchl), & 
     2147         & zglam(2,2,ilogchl), & 
     2148         & zgphi(2,2,ilogchl), & 
     2149         & zmask(2,2,ilogchl), & 
     2150         & zlogchll(2,2,ilogchl)  & 
     2151         & ) 
     2152       
     2153      DO jobs = logchldatqc%nsurfup + 1, logchldatqc%nsurfup + ilogchl 
     2154         iobs = jobs - logchldatqc%nsurfup 
     2155         igrdi(1,1,iobs) = logchldatqc%mi(jobs)-1 
     2156         igrdj(1,1,iobs) = logchldatqc%mj(jobs)-1 
     2157         igrdi(1,2,iobs) = logchldatqc%mi(jobs)-1 
     2158         igrdj(1,2,iobs) = logchldatqc%mj(jobs) 
     2159         igrdi(2,1,iobs) = logchldatqc%mi(jobs) 
     2160         igrdj(2,1,iobs) = logchldatqc%mj(jobs)-1 
     2161         igrdi(2,2,iobs) = logchldatqc%mi(jobs) 
     2162         igrdj(2,2,iobs) = logchldatqc%mj(jobs) 
     2163      END DO 
     2164       
     2165      CALL obs_int_comm_2d( 2, 2, ilogchl, & 
     2166         &                  igrdi, igrdj, glamt, zglam ) 
     2167      CALL obs_int_comm_2d( 2, 2, ilogchl, & 
     2168         &                  igrdi, igrdj, gphit, zgphi ) 
     2169      CALL obs_int_comm_2d( 2, 2, ilogchl, & 
     2170         &                  igrdi, igrdj, plogchlmask, zmask ) 
     2171      CALL obs_int_comm_2d( 2, 2, ilogchl, & 
     2172         &                  igrdi, igrdj, plogchln, zlogchll ) 
     2173       
     2174      DO jobs = logchldatqc%nsurfup + 1, logchldatqc%nsurfup + ilogchl 
     2175          
     2176         iobs = jobs - logchldatqc%nsurfup 
     2177          
     2178         IF ( kt /= logchldatqc%mstp(jobs) ) THEN 
     2179             
     2180            IF(lwp) THEN 
     2181               WRITE(numout,*) 
     2182               WRITE(numout,*) ' E R R O R : Observation',              & 
     2183                  &            ' time step is not consistent with the', & 
     2184                  &            ' model time step' 
     2185               WRITE(numout,*) ' =========' 
     2186               WRITE(numout,*) 
     2187               WRITE(numout,*) ' Record  = ', jobs,                & 
     2188                  &            ' kt      = ', kt,                  & 
     2189                  &            ' mstp    = ', logchldatqc%mstp(jobs), & 
     2190                  &            ' ntyp    = ', logchldatqc%ntyp(jobs) 
     2191            ENDIF 
     2192            CALL ctl_stop( 'obs_logchl_opt', 'Inconsistent time' ) 
     2193             
     2194         ENDIF 
     2195          
     2196         zlam = logchldatqc%rlam(jobs) 
     2197         zphi = logchldatqc%rphi(jobs) 
     2198          
     2199         ! Get weights to interpolate the model logchl to the observation point 
     2200         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
     2201            &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     2202            &                   zmask(:,:,iobs), zweig, zobsmask ) 
     2203          
     2204         ! ... Interpolate the model logchl to the observation point 
     2205         CALL obs_int_h2d( 1, 1,      & 
     2206            &              zweig, zlogchll(:,:,iobs),  zext ) 
     2207          
     2208         logchldatqc%rmod(jobs,1) = zext(1) 
     2209          
     2210      END DO 
     2211       
     2212      ! Deallocate the data for interpolation 
     2213      DEALLOCATE( & 
     2214         & igrdi,    & 
     2215         & igrdj,    & 
     2216         & zglam,    & 
     2217         & zgphi,    & 
     2218         & zmask,    & 
     2219         & zlogchll  & 
     2220         & ) 
     2221       
     2222      logchldatqc%nsurfup = logchldatqc%nsurfup + ilogchl 
     2223 
     2224   END SUBROUTINE obs_logchl_opt 
     2225 
     2226   SUBROUTINE obs_spm_opt( spmdatqc, kt, kpi, kpj, kit000, & 
     2227      &                    pspmn, pspmmask, k2dint ) 
     2228 
     2229      !!----------------------------------------------------------------------- 
     2230      !! 
     2231      !!                     ***  ROUTINE obs_spm_opt  *** 
     2232      !! 
     2233      !! ** Purpose : Compute the model counterpart of spm 
     2234      !!              data by interpolating from the model grid to the  
     2235      !!              observation point. 
     2236      !! 
     2237      !! ** Method  : Linearly interpolate to each observation point using  
     2238      !!              the model values at the corners of the surrounding grid box. 
     2239      !! 
     2240      !!    The now model spm is first computed at the obs (lon, lat) point. 
     2241      !! 
     2242      !!    Several horizontal interpolation schemes are available: 
     2243      !!        - distance-weighted (great circle) (k2dint = 0) 
     2244      !!        - distance-weighted (small angle)  (k2dint = 1) 
     2245      !!        - bilinear (geographical grid)     (k2dint = 2) 
     2246      !!        - bilinear (quadrilateral grid)    (k2dint = 3) 
     2247      !!        - polynomial (quadrilateral grid)  (k2dint = 4) 
     2248      !! 
     2249      !! 
     2250      !! ** Action  : 
     2251      !! 
     2252      !! History : 
     2253      !!       
     2254      !!----------------------------------------------------------------------- 
     2255 
     2256      !! * Modules used 
     2257      USE obs_surf_def  ! Definition of storage space for surface observations 
     2258 
     2259      IMPLICIT NONE 
     2260 
     2261      !! * Arguments 
     2262      TYPE(obs_surf), INTENT(INOUT) :: spmdatqc     ! Subset of surface data not failing screening 
     2263      INTEGER, INTENT(IN) :: kt       ! Time step 
     2264      INTEGER, INTENT(IN) :: kpi      ! Model grid parameters 
     2265      INTEGER, INTENT(IN) :: kpj 
     2266      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step  
     2267                                      !   (kit000-1 = restart time) 
     2268      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header) 
     2269      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
     2270         & pspmn,  &    ! Model spm field 
     2271         & pspmmask     ! Land-sea mask 
     2272          
     2273      !! * Local declarations 
     2274      INTEGER :: ji 
     2275      INTEGER :: jj 
     2276      INTEGER :: jobs 
     2277      INTEGER :: inrc 
     2278      INTEGER :: ispm 
     2279      INTEGER :: iobs 
     2280        
     2281      REAL(KIND=wp) :: zlam 
     2282      REAL(KIND=wp) :: zphi 
     2283      REAL(KIND=wp) :: zext(1), zobsmask(1) 
     2284      REAL(kind=wp), DIMENSION(2,2,1) :: & 
     2285         & zweig 
     2286      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
     2287         & zmask, & 
     2288         & zspml, & 
     2289         & zglam, & 
     2290         & zgphi 
     2291      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
     2292         & igrdi, & 
     2293         & igrdj 
     2294 
     2295      !------------------------------------------------------------------------ 
     2296      ! Local initialization  
     2297      !------------------------------------------------------------------------ 
     2298      ! ... Record and data counters 
     2299      inrc = kt - kit000 + 2 
     2300      ispm = spmdatqc%nsstp(inrc) 
     2301 
     2302      ! Get the data for interpolation 
     2303       
     2304      ALLOCATE( & 
     2305         & igrdi(2,2,ispm), & 
     2306         & igrdj(2,2,ispm), & 
     2307         & zglam(2,2,ispm), & 
     2308         & zgphi(2,2,ispm), & 
     2309         & zmask(2,2,ispm), & 
     2310         & zspml(2,2,ispm)  & 
     2311         & ) 
     2312       
     2313      DO jobs = spmdatqc%nsurfup + 1, spmdatqc%nsurfup + ispm 
     2314         iobs = jobs - spmdatqc%nsurfup 
     2315         igrdi(1,1,iobs) = spmdatqc%mi(jobs)-1 
     2316         igrdj(1,1,iobs) = spmdatqc%mj(jobs)-1 
     2317         igrdi(1,2,iobs) = spmdatqc%mi(jobs)-1 
     2318         igrdj(1,2,iobs) = spmdatqc%mj(jobs) 
     2319         igrdi(2,1,iobs) = spmdatqc%mi(jobs) 
     2320         igrdj(2,1,iobs) = spmdatqc%mj(jobs)-1 
     2321         igrdi(2,2,iobs) = spmdatqc%mi(jobs) 
     2322         igrdj(2,2,iobs) = spmdatqc%mj(jobs) 
     2323      END DO 
     2324       
     2325      CALL obs_int_comm_2d( 2, 2, ispm, & 
     2326         &                  igrdi, igrdj, glamt, zglam ) 
     2327      CALL obs_int_comm_2d( 2, 2, ispm, & 
     2328         &                  igrdi, igrdj, gphit, zgphi ) 
     2329      CALL obs_int_comm_2d( 2, 2, ispm, & 
     2330         &                  igrdi, igrdj, pspmmask, zmask ) 
     2331      CALL obs_int_comm_2d( 2, 2, ispm, & 
     2332         &                  igrdi, igrdj, pspmn, zspml ) 
     2333       
     2334      DO jobs = spmdatqc%nsurfup + 1, spmdatqc%nsurfup + ispm 
     2335          
     2336         iobs = jobs - spmdatqc%nsurfup 
     2337          
     2338         IF ( kt /= spmdatqc%mstp(jobs) ) THEN 
     2339             
     2340            IF(lwp) THEN 
     2341               WRITE(numout,*) 
     2342               WRITE(numout,*) ' E R R O R : Observation',              & 
     2343                  &            ' time step is not consistent with the', & 
     2344                  &            ' model time step' 
     2345               WRITE(numout,*) ' =========' 
     2346               WRITE(numout,*) 
     2347               WRITE(numout,*) ' Record  = ', jobs,                & 
     2348                  &            ' kt      = ', kt,                  & 
     2349                  &            ' mstp    = ', spmdatqc%mstp(jobs), & 
     2350                  &            ' ntyp    = ', spmdatqc%ntyp(jobs) 
     2351            ENDIF 
     2352            CALL ctl_stop( 'obs_spm_opt', 'Inconsistent time' ) 
     2353             
     2354         ENDIF 
     2355          
     2356         zlam = spmdatqc%rlam(jobs) 
     2357         zphi = spmdatqc%rphi(jobs) 
     2358          
     2359         ! Get weights to interpolate the model spm to the observation point 
     2360         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
     2361            &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     2362            &                   zmask(:,:,iobs), zweig, zobsmask ) 
     2363          
     2364         ! ... Interpolate the model spm to the observation point 
     2365         CALL obs_int_h2d( 1, 1,      & 
     2366            &              zweig, zspml(:,:,iobs),  zext ) 
     2367          
     2368         spmdatqc%rmod(jobs,1) = zext(1) 
     2369          
     2370      END DO 
     2371       
     2372      ! Deallocate the data for interpolation 
     2373      DEALLOCATE( & 
     2374         & igrdi,    & 
     2375         & igrdj,    & 
     2376         & zglam,    & 
     2377         & zgphi,    & 
     2378         & zmask,    & 
     2379         & zspml  & 
     2380         & ) 
     2381       
     2382      spmdatqc%nsurfup = spmdatqc%nsurfup + ispm 
     2383 
     2384   END SUBROUTINE obs_spm_opt 
     2385 
     2386   SUBROUTINE obs_fco2_opt( fco2datqc, kt, kpi, kpj, kit000, & 
     2387      &                    pfco2n, pfco2mask, k2dint ) 
     2388 
     2389      !!----------------------------------------------------------------------- 
     2390      !! 
     2391      !!                     ***  ROUTINE obs_fco2_opt  *** 
     2392      !! 
     2393      !! ** Purpose : Compute the model counterpart of fco2 
     2394      !!              data by interpolating from the model grid to the  
     2395      !!              observation point. 
     2396      !! 
     2397      !! ** Method  : Linearly interpolate to each observation point using  
     2398      !!              the model values at the corners of the surrounding grid box. 
     2399      !! 
     2400      !!    The now model fco2 is first computed at the obs (lon, lat) point. 
     2401      !! 
     2402      !!    Several horizontal interpolation schemes are available: 
     2403      !!        - distance-weighted (great circle) (k2dint = 0) 
     2404      !!        - distance-weighted (small angle)  (k2dint = 1) 
     2405      !!        - bilinear (geographical grid)     (k2dint = 2) 
     2406      !!        - bilinear (quadrilateral grid)    (k2dint = 3) 
     2407      !!        - polynomial (quadrilateral grid)  (k2dint = 4) 
     2408      !! 
     2409      !! 
     2410      !! ** Action  : 
     2411      !! 
     2412      !! History : 
     2413      !!       
     2414      !!----------------------------------------------------------------------- 
     2415 
     2416      !! * Modules used 
     2417      USE obs_surf_def  ! Definition of storage space for surface observations 
     2418 
     2419      IMPLICIT NONE 
     2420 
     2421      !! * Arguments 
     2422      TYPE(obs_surf), INTENT(INOUT) :: fco2datqc     ! Subset of surface data not failing screening 
     2423      INTEGER, INTENT(IN) :: kt       ! Time step 
     2424      INTEGER, INTENT(IN) :: kpi      ! Model grid parameters 
     2425      INTEGER, INTENT(IN) :: kpj 
     2426      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step  
     2427                                      !   (kit000-1 = restart time) 
     2428      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header) 
     2429      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
     2430         & pfco2n,  &    ! Model fco2 field 
     2431         & pfco2mask     ! Land-sea mask 
     2432          
     2433      !! * Local declarations 
     2434      INTEGER :: ji 
     2435      INTEGER :: jj 
     2436      INTEGER :: jobs 
     2437      INTEGER :: inrc 
     2438      INTEGER :: ifco2 
     2439      INTEGER :: iobs 
     2440        
     2441      REAL(KIND=wp) :: zlam 
     2442      REAL(KIND=wp) :: zphi 
     2443      REAL(KIND=wp) :: zext(1), zobsmask(1) 
     2444      REAL(kind=wp), DIMENSION(2,2,1) :: & 
     2445         & zweig 
     2446      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
     2447         & zmask, & 
     2448         & zfco2l, & 
     2449         & zglam, & 
     2450         & zgphi 
     2451      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
     2452         & igrdi, & 
     2453         & igrdj 
     2454 
     2455      !------------------------------------------------------------------------ 
     2456      ! Local initialization  
     2457      !------------------------------------------------------------------------ 
     2458      ! ... Record and data counters 
     2459      inrc = kt - kit000 + 2 
     2460      ifco2 = fco2datqc%nsstp(inrc) 
     2461 
     2462      ! Get the data for interpolation 
     2463       
     2464      ALLOCATE( & 
     2465         & igrdi(2,2,ifco2), & 
     2466         & igrdj(2,2,ifco2), & 
     2467         & zglam(2,2,ifco2), & 
     2468         & zgphi(2,2,ifco2), & 
     2469         & zmask(2,2,ifco2), & 
     2470         & zfco2l(2,2,ifco2)  & 
     2471         & ) 
     2472       
     2473      DO jobs = fco2datqc%nsurfup + 1, fco2datqc%nsurfup + ifco2 
     2474         iobs = jobs - fco2datqc%nsurfup 
     2475         igrdi(1,1,iobs) = fco2datqc%mi(jobs)-1 
     2476         igrdj(1,1,iobs) = fco2datqc%mj(jobs)-1 
     2477         igrdi(1,2,iobs) = fco2datqc%mi(jobs)-1 
     2478         igrdj(1,2,iobs) = fco2datqc%mj(jobs) 
     2479         igrdi(2,1,iobs) = fco2datqc%mi(jobs) 
     2480         igrdj(2,1,iobs) = fco2datqc%mj(jobs)-1 
     2481         igrdi(2,2,iobs) = fco2datqc%mi(jobs) 
     2482         igrdj(2,2,iobs) = fco2datqc%mj(jobs) 
     2483      END DO 
     2484       
     2485      CALL obs_int_comm_2d( 2, 2, ifco2, & 
     2486         &                  igrdi, igrdj, glamt, zglam ) 
     2487      CALL obs_int_comm_2d( 2, 2, ifco2, & 
     2488         &                  igrdi, igrdj, gphit, zgphi ) 
     2489      CALL obs_int_comm_2d( 2, 2, ifco2, & 
     2490         &                  igrdi, igrdj, pfco2mask, zmask ) 
     2491      CALL obs_int_comm_2d( 2, 2, ifco2, & 
     2492         &                  igrdi, igrdj, pfco2n, zfco2l ) 
     2493       
     2494      DO jobs = fco2datqc%nsurfup + 1, fco2datqc%nsurfup + ifco2 
     2495          
     2496         iobs = jobs - fco2datqc%nsurfup 
     2497          
     2498         IF ( kt /= fco2datqc%mstp(jobs) ) THEN 
     2499             
     2500            IF(lwp) THEN 
     2501               WRITE(numout,*) 
     2502               WRITE(numout,*) ' E R R O R : Observation',              & 
     2503                  &            ' time step is not consistent with the', & 
     2504                  &            ' model time step' 
     2505               WRITE(numout,*) ' =========' 
     2506               WRITE(numout,*) 
     2507               WRITE(numout,*) ' Record  = ', jobs,                & 
     2508                  &            ' kt      = ', kt,                  & 
     2509                  &            ' mstp    = ', fco2datqc%mstp(jobs), & 
     2510                  &            ' ntyp    = ', fco2datqc%ntyp(jobs) 
     2511            ENDIF 
     2512            CALL ctl_stop( 'obs_fco2_opt', 'Inconsistent time' ) 
     2513             
     2514         ENDIF 
     2515          
     2516         zlam = fco2datqc%rlam(jobs) 
     2517         zphi = fco2datqc%rphi(jobs) 
     2518          
     2519         ! Get weights to interpolate the model fco2 to the observation point 
     2520         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
     2521            &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     2522            &                   zmask(:,:,iobs), zweig, zobsmask ) 
     2523          
     2524         ! ... Interpolate the model fco2 to the observation point 
     2525         CALL obs_int_h2d( 1, 1,      & 
     2526            &              zweig, zfco2l(:,:,iobs),  zext ) 
     2527          
     2528         fco2datqc%rmod(jobs,1) = zext(1) 
     2529          
     2530      END DO 
     2531       
     2532      ! Deallocate the data for interpolation 
     2533      DEALLOCATE( & 
     2534         & igrdi,    & 
     2535         & igrdj,    & 
     2536         & zglam,    & 
     2537         & zgphi,    & 
     2538         & zmask,    & 
     2539         & zfco2l  & 
     2540         & ) 
     2541       
     2542      fco2datqc%nsurfup = fco2datqc%nsurfup + ifco2 
     2543 
     2544   END SUBROUTINE obs_fco2_opt 
     2545 
     2546   SUBROUTINE obs_pco2_opt( pco2datqc, kt, kpi, kpj, kit000, & 
     2547      &                    ppco2n, ppco2mask, k2dint ) 
     2548 
     2549      !!----------------------------------------------------------------------- 
     2550      !! 
     2551      !!                     ***  ROUTINE obs_pco2_opt  *** 
     2552      !! 
     2553      !! ** Purpose : Compute the model counterpart of pco2 
     2554      !!              data by interpolating from the model grid to the  
     2555      !!              observation point. 
     2556      !! 
     2557      !! ** Method  : Linearly interpolate to each observation point using  
     2558      !!              the model values at the corners of the surrounding grid box. 
     2559      !! 
     2560      !!    The now model pco2 is first computed at the obs (lon, lat) point. 
     2561      !! 
     2562      !!    Several horizontal interpolation schemes are available: 
     2563      !!        - distance-weighted (great circle) (k2dint = 0) 
     2564      !!        - distance-weighted (small angle)  (k2dint = 1) 
     2565      !!        - bilinear (geographical grid)     (k2dint = 2) 
     2566      !!        - bilinear (quadrilateral grid)    (k2dint = 3) 
     2567      !!        - polynomial (quadrilateral grid)  (k2dint = 4) 
     2568      !! 
     2569      !! 
     2570      !! ** Action  : 
     2571      !! 
     2572      !! History : 
     2573      !!       
     2574      !!----------------------------------------------------------------------- 
     2575 
     2576      !! * Modules used 
     2577      USE obs_surf_def  ! Definition of storage space for surface observations 
     2578 
     2579      IMPLICIT NONE 
     2580 
     2581      !! * Arguments 
     2582      TYPE(obs_surf), INTENT(INOUT) :: pco2datqc     ! Subset of surface data not failing screening 
     2583      INTEGER, INTENT(IN) :: kt       ! Time step 
     2584      INTEGER, INTENT(IN) :: kpi      ! Model grid parameters 
     2585      INTEGER, INTENT(IN) :: kpj 
     2586      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step  
     2587                                      !   (kit000-1 = restart time) 
     2588      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header) 
     2589      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
     2590         & ppco2n,  &    ! Model pco2 field 
     2591         & ppco2mask     ! Land-sea mask 
     2592          
     2593      !! * Local declarations 
     2594      INTEGER :: ji 
     2595      INTEGER :: jj 
     2596      INTEGER :: jobs 
     2597      INTEGER :: inrc 
     2598      INTEGER :: ipco2 
     2599      INTEGER :: iobs 
     2600        
     2601      REAL(KIND=wp) :: zlam 
     2602      REAL(KIND=wp) :: zphi 
     2603      REAL(KIND=wp) :: zext(1), zobsmask(1) 
     2604      REAL(kind=wp), DIMENSION(2,2,1) :: & 
     2605         & zweig 
     2606      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
     2607         & zmask, & 
     2608         & zpco2l, & 
     2609         & zglam, & 
     2610         & zgphi 
     2611      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
     2612         & igrdi, & 
     2613         & igrdj 
     2614 
     2615      !------------------------------------------------------------------------ 
     2616      ! Local initialization  
     2617      !------------------------------------------------------------------------ 
     2618      ! ... Record and data counters 
     2619      inrc = kt - kit000 + 2 
     2620      ipco2 = pco2datqc%nsstp(inrc) 
     2621 
     2622      ! Get the data for interpolation 
     2623       
     2624      ALLOCATE( & 
     2625         & igrdi(2,2,ipco2), & 
     2626         & igrdj(2,2,ipco2), & 
     2627         & zglam(2,2,ipco2), & 
     2628         & zgphi(2,2,ipco2), & 
     2629         & zmask(2,2,ipco2), & 
     2630         & zpco2l(2,2,ipco2)  & 
     2631         & ) 
     2632       
     2633      DO jobs = pco2datqc%nsurfup + 1, pco2datqc%nsurfup + ipco2 
     2634         iobs = jobs - pco2datqc%nsurfup 
     2635         igrdi(1,1,iobs) = pco2datqc%mi(jobs)-1 
     2636         igrdj(1,1,iobs) = pco2datqc%mj(jobs)-1 
     2637         igrdi(1,2,iobs) = pco2datqc%mi(jobs)-1 
     2638         igrdj(1,2,iobs) = pco2datqc%mj(jobs) 
     2639         igrdi(2,1,iobs) = pco2datqc%mi(jobs) 
     2640         igrdj(2,1,iobs) = pco2datqc%mj(jobs)-1 
     2641         igrdi(2,2,iobs) = pco2datqc%mi(jobs) 
     2642         igrdj(2,2,iobs) = pco2datqc%mj(jobs) 
     2643      END DO 
     2644       
     2645      CALL obs_int_comm_2d( 2, 2, ipco2, & 
     2646         &                  igrdi, igrdj, glamt, zglam ) 
     2647      CALL obs_int_comm_2d( 2, 2, ipco2, & 
     2648         &                  igrdi, igrdj, gphit, zgphi ) 
     2649      CALL obs_int_comm_2d( 2, 2, ipco2, & 
     2650         &                  igrdi, igrdj, ppco2mask, zmask ) 
     2651      CALL obs_int_comm_2d( 2, 2, ipco2, & 
     2652         &                  igrdi, igrdj, ppco2n, zpco2l ) 
     2653       
     2654      DO jobs = pco2datqc%nsurfup + 1, pco2datqc%nsurfup + ipco2 
     2655          
     2656         iobs = jobs - pco2datqc%nsurfup 
     2657          
     2658         IF ( kt /= pco2datqc%mstp(jobs) ) THEN 
     2659             
     2660            IF(lwp) THEN 
     2661               WRITE(numout,*) 
     2662               WRITE(numout,*) ' E R R O R : Observation',              & 
     2663                  &            ' time step is not consistent with the', & 
     2664                  &            ' model time step' 
     2665               WRITE(numout,*) ' =========' 
     2666               WRITE(numout,*) 
     2667               WRITE(numout,*) ' Record  = ', jobs,                & 
     2668                  &            ' kt      = ', kt,                  & 
     2669                  &            ' mstp    = ', pco2datqc%mstp(jobs), & 
     2670                  &            ' ntyp    = ', pco2datqc%ntyp(jobs) 
     2671            ENDIF 
     2672            CALL ctl_stop( 'obs_pco2_opt', 'Inconsistent time' ) 
     2673             
     2674         ENDIF 
     2675          
     2676         zlam = pco2datqc%rlam(jobs) 
     2677         zphi = pco2datqc%rphi(jobs) 
     2678          
     2679         ! Get weights to interpolate the model pco2 to the observation point 
     2680         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
     2681            &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     2682            &                   zmask(:,:,iobs), zweig, zobsmask ) 
     2683          
     2684         ! ... Interpolate the model pco2 to the observation point 
     2685         CALL obs_int_h2d( 1, 1,      & 
     2686            &              zweig, zpco2l(:,:,iobs),  zext ) 
     2687          
     2688         pco2datqc%rmod(jobs,1) = zext(1) 
     2689          
     2690      END DO 
     2691       
     2692      ! Deallocate the data for interpolation 
     2693      DEALLOCATE( & 
     2694         & igrdi,    & 
     2695         & igrdj,    & 
     2696         & zglam,    & 
     2697         & zgphi,    & 
     2698         & zmask,    & 
     2699         & zpco2l  & 
     2700         & ) 
     2701       
     2702      pco2datqc%nsurfup = pco2datqc%nsurfup + ipco2 
     2703 
     2704   END SUBROUTINE obs_pco2_opt 
     2705 
    20542706END MODULE obs_oper 
    20552707 
Note: See TracChangeset for help on using the changeset viewer.