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 – NEMO

Changeset 7713


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.

Location:
branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS
Files:
4 edited
16 copied

Legend:

Unmodified
Added
Removed

  • branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90

    r6990 r7713  
    2121   USE par_oce 
    2222   USE dom_oce                  ! Ocean space and time domain variables 
     23   USE obs_const, ONLY: obfillflt ! Fill value 
    2324   USE obs_fbm, ONLY: ln_cl4    ! Class 4 diagnostic switch 
    2425   USE obs_read_prof            ! Reading and allocation of observations (Coriolis) 
    … …  
    2930   USE obs_read_seaice          ! Reading and allocation of Sea Ice observations   
    3031   USE obs_read_vel             ! Reading and allocation of velocity component observations 
     32   USE obs_read_logchl          ! Reading and allocation of logchl observations 
     33   USE obs_read_spm             ! Reading and allocation of spm observations 
     34   USE obs_read_fco2            ! Reading and allocation of fco2 observations 
     35   USE obs_read_pco2            ! Reading and allocation of pco2 observations 
    3136   USE obs_prep                 ! Preparation of obs. (grid search etc). 
    3237   USE obs_oper                 ! Observation operators 
    … …  
    4045   USE obs_sst                  ! SST data storage 
    4146   USE obs_seaice               ! Sea Ice data storage 
     47   USE obs_logchl               ! logchl data storage 
     48   USE obs_spm                  ! spm data storage 
     49   USE obs_fco2                 ! fco2 data storage 
     50   USE obs_pco2                 ! pco2 data storage 
    4251   USE obs_types                ! Definitions for observation types 
    4352   USE mpp_map                  ! MPP mapping 
    … …  
    8190   LOGICAL, PUBLIC :: ln_velhradcp   !: Logical switch for raw high freq netCDF ADCP vel. data  
    8291   LOGICAL, PUBLIC :: ln_velfb       !: Logical switch for velocities from feedback files 
     92   LOGICAL, PUBLIC :: ln_logchl      !: Logical switch for log10(chlorophyll) 
     93   LOGICAL, PUBLIC :: ln_logchlfb    !: Logical switch for logchl from feedback files 
     94   LOGICAL, PUBLIC :: ln_spm         !: Logical switch for spm 
     95   LOGICAL, PUBLIC :: ln_spmfb       !: Logical switch for spm from feedback files 
     96   LOGICAL, PUBLIC :: ln_fco2        !: Logical switch for fco2 
     97   LOGICAL, PUBLIC :: ln_fco2fb      !: Logical switch for fco2 from feedback files 
     98   LOGICAL, PUBLIC :: ln_pco2        !: Logical switch for pco2 
     99   LOGICAL, PUBLIC :: ln_pco2fb      !: Logical switch for pco2 from feedback files 
    83100   LOGICAL, PUBLIC :: ln_ssh         !: Logical switch for sea surface height 
    84101   LOGICAL, PUBLIC :: ln_sss         !: Logical switch for sea surface salinity 
    … …  
    164181      CHARACTER(len=128) :: velhradcpfiles(MaxNumFiles) 
    165182      CHARACTER(len=128) :: velfbfiles(MaxNumFiles) 
     183      CHARACTER(len=128) :: logchlfiles(MaxNumFiles) 
     184      CHARACTER(len=128) :: logchlfbfiles(MaxNumFiles) 
     185      CHARACTER(len=128) :: spmfiles(MaxNumFiles) 
     186      CHARACTER(len=128) :: spmfbfiles(MaxNumFiles) 
     187      CHARACTER(len=128) :: fco2files(MaxNumFiles) 
     188      CHARACTER(len=128) :: fco2fbfiles(MaxNumFiles) 
     189      CHARACTER(len=128) :: pco2files(MaxNumFiles) 
     190      CHARACTER(len=128) :: pco2fbfiles(MaxNumFiles) 
    166191      CHARACTER(LEN=128) :: reysstname 
    167192      CHARACTER(LEN=12)  :: reysstfmt 
    … …  
    189214         &            ln_velhradcp, velhradcpfiles,                   & 
    190215         &            ln_velfb, velfbfiles, ln_velfb_av,              & 
     216         &            ln_logchl, ln_logchlfb,                         & 
     217         &            logchlfiles, logchlfbfiles,                     & 
     218         &            ln_spm, ln_spmfb,                               & 
     219         &            spmfiles, spmfbfiles,                           & 
     220         &            ln_fco2, ln_fco2fb,                             & 
     221         &            fco2files, fco2fbfiles,                         & 
     222         &            ln_pco2, ln_pco2fb,                             & 
     223         &            pco2files, pco2fbfiles,                         & 
    191224         &            ln_profb_enatim, ln_ignmis, ln_cl4,             & 
    192225         &            ln_sstbias, sstbias_files 
    … …  
    210243      INTEGER :: jnumvelhradcp    
    211244      INTEGER :: jnumvelfb 
     245      INTEGER :: jnumlogchl 
     246      INTEGER :: jnumlogchlfb 
     247      INTEGER :: jnumspm 
     248      INTEGER :: jnumspmfb 
     249      INTEGER :: jnumfco2 
     250      INTEGER :: jnumfco2fb 
     251      INTEGER :: jnumpco2 
     252      INTEGER :: jnumpco2fb 
    212253      INTEGER :: ji 
    213254      INTEGER :: jset 
    … …  
    218259      ! Read namelist parameters 
    219260      !----------------------------------------------------------------------- 
     261 
     262      ln_logchl   = .FALSE. 
     263      ln_logchlfb = .FALSE. 
     264      ln_spm      = .FALSE. 
     265      ln_spmfb    = .FALSE. 
     266      ln_fco2     = .FALSE. 
     267      ln_fco2fb   = .FALSE. 
     268      ln_pco2     = .FALSE. 
     269      ln_pco2fb   = .FALSE. 
    220270       
    221271      !Initalise all values in namelist arrays 
    … …  
    238288      velcurfiles(:) = '' 
    239289      veladcpfiles(:) = '' 
     290      logchlfiles(:) = '' 
     291      logchlfbfiles(:) = '' 
     292      spmfiles(:) = '' 
     293      spmfbfiles(:) = '' 
     294      fco2files(:) = '' 
     295      fco2fbfiles(:) = '' 
     296      pco2files(:) = '' 
     297      pco2fbfiles(:) = '' 
    240298      sstbias_files(:) = '' 
    241299      endailyavtypes(:) = -1 
    … …  
    337395         jnumvelfb = COUNT(lmask) 
    338396         lmask(:) = .FALSE. 
     397      ENDIF 
     398      IF (ln_logchl) THEN 
     399         lmask(:) = .FALSE. 
     400         WHERE (logchlfiles(:) /= '') lmask(:) = .TRUE. 
     401         jnumlogchl = COUNT(lmask) 
     402      ENDIF 
     403      IF (ln_logchlfb) THEN 
     404         lmask(:) = .FALSE. 
     405         WHERE (logchlfbfiles(:) /= '') lmask(:) = .TRUE. 
     406         jnumlogchlfb = COUNT(lmask) 
     407      ENDIF 
     408      IF (ln_spm) THEN 
     409         lmask(:) = .FALSE. 
     410         WHERE (spmfiles(:) /= '') lmask(:) = .TRUE. 
     411         jnumspm = COUNT(lmask) 
     412      ENDIF 
     413      IF (ln_spmfb) THEN 
     414         lmask(:) = .FALSE. 
     415         WHERE (spmfbfiles(:) /= '') lmask(:) = .TRUE. 
     416         jnumspmfb = COUNT(lmask) 
     417      ENDIF 
     418      IF (ln_fco2) THEN 
     419         lmask(:) = .FALSE. 
     420         WHERE (fco2files(:) /= '') lmask(:) = .TRUE. 
     421         jnumfco2 = COUNT(lmask) 
     422      ENDIF 
     423      IF (ln_fco2fb) THEN 
     424         lmask(:) = .FALSE. 
     425         WHERE (fco2fbfiles(:) /= '') lmask(:) = .TRUE. 
     426         jnumfco2fb = COUNT(lmask) 
     427      ENDIF 
     428      IF (ln_pco2) THEN 
     429         lmask(:) = .FALSE. 
     430         WHERE (pco2files(:) /= '') lmask(:) = .TRUE. 
     431         jnumpco2 = COUNT(lmask) 
     432      ENDIF 
     433      IF (ln_pco2fb) THEN 
     434         lmask(:) = .FALSE. 
     435         WHERE (pco2fbfiles(:) /= '') lmask(:) = .TRUE. 
     436         jnumpco2fb = COUNT(lmask) 
    339437      ENDIF 
    340438       
    … …  
    368466         WRITE(numout,*) '             Logical switch for velocity high freq. ADCP  ln_velhradcp = ', ln_velhradcp 
    369467         WRITE(numout,*) '             Logical switch for feedback velocity data        ln_velfb = ', ln_velfb 
     468         WRITE(numout,*) '             Logical switch for logchl observations          ln_logchl = ', ln_logchl 
     469         WRITE(numout,*) '             Logical switch for feedback logchl data       ln_logchlfb = ', ln_logchlfb 
     470         WRITE(numout,*) '             Logical switch for spm observations                ln_spm = ', ln_spm 
     471         WRITE(numout,*) '             Logical switch for feedback spm data             ln_spmfb = ', ln_spmfb 
     472         WRITE(numout,*) '             Logical switch for fco2 observations              ln_fco2 = ', ln_fco2 
     473         WRITE(numout,*) '             Logical switch for pco2 observations              ln_pco2 = ', ln_pco2 
     474         WRITE(numout,*) '             Logical switch for feedback pco2 data           ln_pco2fb = ', ln_pco2fb 
     475         WRITE(numout,*) '             Logical switch for feedback fco2 data           ln_fco2fb = ', ln_fco2fb 
    370476         WRITE(numout,*) '             Global distribtion of observations         ln_grid_global = ',ln_grid_global 
    371477         WRITE(numout,*) & 
    … …  
    464570                     TRIM(velfbfiles(ji)) 
    465571               ENDIF 
     572            END DO 
     573         ENDIF 
     574         IF (ln_logchl) THEN 
     575            DO ji = 1, jnumlogchl 
     576               WRITE(numout,'(1X,2A)') '             logchl input observation file name        logchlfiles = ', & 
     577                  TRIM(logchlfiles(ji)) 
     578            END DO 
     579         ENDIF 
     580         IF (ln_logchlfb) THEN 
     581            DO ji = 1, jnumlogchlfb 
     582               WRITE(numout,'(1X,2A)') '        Feedback logchl input observation file name  logchlfbfiles = ', & 
     583                  TRIM(logchlfbfiles(ji)) 
     584            END DO 
     585         ENDIF 
     586         IF (ln_spm) THEN 
     587            DO ji = 1, jnumspm 
     588               WRITE(numout,'(1X,2A)') '             spm input observation file name              spmfiles = ', & 
     589                  TRIM(spmfiles(ji)) 
     590            END DO 
     591         ENDIF 
     592         IF (ln_spmfb) THEN 
     593            DO ji = 1, jnumspmfb 
     594               WRITE(numout,'(1X,2A)') '             Feedback spm input observation file name   spmfbfiles = ', & 
     595                  TRIM(spmfbfiles(ji)) 
     596            END DO 
     597         ENDIF 
     598         IF (ln_fco2) THEN 
     599            DO ji = 1, jnumfco2 
     600               WRITE(numout,'(1X,2A)') '             fco2 input observation file name            fco2files = ', & 
     601                  TRIM(fco2files(ji)) 
     602            END DO 
     603         ENDIF 
     604         IF (ln_fco2fb) THEN 
     605            DO ji = 1, jnumfco2fb 
     606               WRITE(numout,'(1X,2A)') '            Feedback fco2 input observation file name  fco2fbfiles = ', & 
     607                  TRIM(fco2fbfiles(ji)) 
     608            END DO 
     609         ENDIF 
     610         IF (ln_pco2) THEN 
     611            DO ji = 1, jnumpco2 
     612               WRITE(numout,'(1X,2A)') '             pco2 input observation file name            pco2files = ', & 
     613                  TRIM(pco2files(ji)) 
     614            END DO 
     615         ENDIF 
     616         IF (ln_pco2fb) THEN 
     617            DO ji = 1, jnumpco2fb 
     618               WRITE(numout,'(1X,2A)') '            Feedback pco2 input observation file name  pco2fbfiles = ', & 
     619                  TRIM(pco2fbfiles(ji)) 
    466620            END DO 
    467621         ENDIF 
    … …  
    501655         & ( .NOT. ln_vel3d ).AND.                                         & 
    502656         & ( .NOT. ln_ssh ).AND.( .NOT. ln_sst ).AND.( .NOT. ln_sss ).AND. & 
    503          & ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ) ) THEN 
     657         & ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ).AND.( .NOT. ln_logchl ).AND. & 
     658         & ( .NOT. ln_spm ).AND.( .NOT. ln_fco2 ).AND.( .NOT. ln_pco2 ) ) THEN 
    504659         IF(lwp) WRITE(numout,cform_war) 
    505660         IF(lwp) WRITE(numout,*) ' key_diaobs is activated but logical flags', & 
    506             &                    ' ln_t3d, ln_s3d, ln_sla, ln_ssh, ln_sst, ln_sss, ln_seaice, ln_vel3d are all set to .FALSE.' 
     661            &                    ' ln_t3d, ln_s3d, ln_sla, ln_ssh, ln_sst, ln_sss, ln_seaice, ln_vel3d,', & 
     662            &                    ' ln_logchl, ln_spm, ln_fco2, ln_pco2 are all set to .FALSE.' 
    507663         nwarn = nwarn + 1 
    508664      ENDIF 
    … …  
    10021158 
    10031159      ENDIF 
     1160 
     1161      !  - log10(chlorophyll) 
     1162       
     1163      IF ( ln_logchl ) THEN 
     1164 
     1165         ! Set the number of variables for logchl to 1 
     1166         nlogchlvars = 1 
     1167 
     1168         ! Set the number of extra variables for logchl to 0 
     1169         nlogchlextr = 0 
     1170          
     1171         IF ( ln_logchlfb ) THEN 
     1172            nlogchlsets = jnumlogchlfb 
     1173         ELSE 
     1174            nlogchlsets = 1 
     1175         ENDIF 
     1176 
     1177         ALLOCATE(logchldata(nlogchlsets)) 
     1178         ALLOCATE(logchldatqc(nlogchlsets)) 
     1179         logchldata(:)%nsurf=0 
     1180         logchldatqc(:)%nsurf=0 
     1181 
     1182         nlogchlsets = 0 
     1183 
     1184         IF ( ln_logchlfb ) THEN             ! Feedback file format 
     1185 
     1186            DO jset = 1, jnumlogchlfb 
     1187             
     1188               nlogchlsets = nlogchlsets + 1 
     1189 
     1190               CALL obs_rea_logchl( 0, logchldata(nlogchlsets), 1, & 
     1191                  &                 logchlfbfiles(jset:jset), & 
     1192                  &                 nlogchlvars, nlogchlextr, nitend-nit000+2, & 
     1193                  &                 dobsini, dobsend, ln_ignmis, .FALSE. ) 
     1194 
     1195               CALL obs_pre_logchl( logchldata(nlogchlsets), logchldatqc(nlogchlsets), & 
     1196                  &                 ln_logchl, ln_nea ) 
     1197             
     1198            ENDDO 
     1199 
     1200         ELSE                              ! Original file format 
     1201 
     1202            nlogchlsets = nlogchlsets + 1 
     1203 
     1204            CALL obs_rea_logchl( 1, logchldata(nlogchlsets), jnumlogchl, & 
     1205               &                 logchlfiles(1:jnumlogchl), & 
     1206               &                 nlogchlvars, nlogchlextr, nitend-nit000+2, & 
     1207               &                 dobsini, dobsend, ln_ignmis, .FALSE. ) 
     1208 
     1209            CALL obs_pre_logchl( logchldata(nlogchlsets), logchldatqc(nlogchlsets), & 
     1210               &                 ln_logchl, ln_nea ) 
     1211 
     1212         ENDIF 
     1213  
     1214      ENDIF 
     1215 
     1216      !  - spm 
     1217       
     1218      IF ( ln_spm ) THEN 
     1219 
     1220         ! Set the number of variables for spm to 1 
     1221         nspmvars = 1 
     1222 
     1223         ! Set the number of extra variables for spm to 0 
     1224         nspmextr = 0 
     1225          
     1226         IF ( ln_spmfb ) THEN 
     1227            nspmsets = jnumspmfb 
     1228         ELSE 
     1229            nspmsets = 1 
     1230         ENDIF 
     1231 
     1232         ALLOCATE(spmdata(nspmsets)) 
     1233         ALLOCATE(spmdatqc(nspmsets)) 
     1234         spmdata(:)%nsurf=0 
     1235         spmdatqc(:)%nsurf=0 
     1236 
     1237         nspmsets = 0 
     1238 
     1239         IF ( ln_spmfb ) THEN             ! Feedback file format 
     1240 
     1241            DO jset = 1, jnumspmfb 
     1242             
     1243               nspmsets = nspmsets + 1 
     1244 
     1245               CALL obs_rea_spm( 0, spmdata(nspmsets), 1, & 
     1246                  &                 spmfbfiles(jset:jset), & 
     1247                  &                 nspmvars, nspmextr, nitend-nit000+2, & 
     1248                  &                 dobsini, dobsend, ln_ignmis, .FALSE. ) 
     1249 
     1250               CALL obs_pre_spm( spmdata(nspmsets), spmdatqc(nspmsets), & 
     1251                  &                 ln_spm, ln_nea ) 
     1252             
     1253            ENDDO 
     1254 
     1255         ELSE                              ! Original file format 
     1256 
     1257            nspmsets = nspmsets + 1 
     1258 
     1259            CALL obs_rea_spm( 1, spmdata(nspmsets), jnumspm, & 
     1260               &                 spmfiles(1:jnumspm), & 
     1261               &                 nspmvars, nspmextr, nitend-nit000+2, & 
     1262               &                 dobsini, dobsend, ln_ignmis, .FALSE. ) 
     1263 
     1264            CALL obs_pre_spm( spmdata(nspmsets), spmdatqc(nspmsets), & 
     1265               &                 ln_spm, ln_nea ) 
     1266 
     1267         ENDIF 
     1268  
     1269      ENDIF 
     1270 
     1271      !  - fco2 
     1272       
     1273      IF ( ln_fco2 ) THEN 
     1274 
     1275         ! Set the number of variables for fco2 to 1 
     1276         nfco2vars = 1 
     1277 
     1278         ! Set the number of extra variables for fco2 to 0 
     1279         nfco2extr = 0 
     1280          
     1281         IF ( ln_fco2fb ) THEN 
     1282            nfco2sets = jnumfco2fb 
     1283         ELSE 
     1284            nfco2sets = 1 
     1285         ENDIF 
     1286 
     1287         ALLOCATE(fco2data(nfco2sets)) 
     1288         ALLOCATE(fco2datqc(nfco2sets)) 
     1289         fco2data(:)%nsurf=0 
     1290         fco2datqc(:)%nsurf=0 
     1291 
     1292         nfco2sets = 0 
     1293 
     1294         IF ( ln_fco2fb ) THEN             ! Feedback file format 
     1295 
     1296            DO jset = 1, jnumfco2fb 
     1297             
     1298               nfco2sets = nfco2sets + 1 
     1299 
     1300               CALL obs_rea_fco2( 0, fco2data(nfco2sets), 1, & 
     1301                  &                 fco2fbfiles(jset:jset), & 
     1302                  &                 nfco2vars, nfco2extr, nitend-nit000+2, & 
     1303                  &                 dobsini, dobsend, ln_ignmis, .FALSE. ) 
     1304 
     1305               CALL obs_pre_fco2( fco2data(nfco2sets), fco2datqc(nfco2sets), & 
     1306                  &                 ln_fco2, ln_nea ) 
     1307             
     1308            ENDDO 
     1309 
     1310         ELSE                              ! Original file format 
     1311 
     1312            nfco2sets = nfco2sets + 1 
     1313 
     1314            CALL obs_rea_fco2( 1, fco2data(nfco2sets), jnumfco2, & 
     1315               &                 fco2files(1:jnumfco2), & 
     1316               &                 nfco2vars, nfco2extr, nitend-nit000+2, & 
     1317               &                 dobsini, dobsend, ln_ignmis, .FALSE. ) 
     1318 
     1319            CALL obs_pre_fco2( fco2data(nfco2sets), fco2datqc(nfco2sets), & 
     1320               &                 ln_fco2, ln_nea ) 
     1321 
     1322         ENDIF 
     1323  
     1324      ENDIF 
     1325 
     1326      !  - pco2 
     1327       
     1328      IF ( ln_pco2 ) THEN 
     1329 
     1330         ! Set the number of variables for pco2 to 1 
     1331         npco2vars = 1 
     1332 
     1333         ! Set the number of extra variables for pco2 to 0 
     1334         npco2extr = 0 
     1335          
     1336         IF ( ln_pco2fb ) THEN 
     1337            npco2sets = jnumpco2fb 
     1338         ELSE 
     1339            npco2sets = 1 
     1340         ENDIF 
     1341 
     1342         ALLOCATE(pco2data(npco2sets)) 
     1343         ALLOCATE(pco2datqc(npco2sets)) 
     1344         pco2data(:)%nsurf=0 
     1345         pco2datqc(:)%nsurf=0 
     1346 
     1347         npco2sets = 0 
     1348 
     1349         IF ( ln_pco2fb ) THEN             ! Feedback file format 
     1350 
     1351            DO jset = 1, jnumpco2fb 
     1352             
     1353               npco2sets = npco2sets + 1 
     1354 
     1355               CALL obs_rea_pco2( 0, pco2data(npco2sets), 1, & 
     1356                  &                 pco2fbfiles(jset:jset), & 
     1357                  &                 npco2vars, npco2extr, nitend-nit000+2, & 
     1358                  &                 dobsini, dobsend, ln_ignmis, .FALSE. ) 
     1359 
     1360               CALL obs_pre_pco2( pco2data(npco2sets), pco2datqc(npco2sets), & 
     1361                  &                 ln_pco2, ln_nea ) 
     1362             
     1363            ENDDO 
     1364 
     1365         ELSE                              ! Original file format 
     1366 
     1367            npco2sets = npco2sets + 1 
     1368 
     1369            CALL obs_rea_pco2( 1, pco2data(npco2sets), jnumpco2, & 
     1370               &                 pco2files(1:jnumpco2), & 
     1371               &                 npco2vars, npco2extr, nitend-nit000+2, & 
     1372               &                 dobsini, dobsend, ln_ignmis, .FALSE. ) 
     1373 
     1374            CALL obs_pre_pco2( pco2data(npco2sets), pco2datqc(npco2sets), & 
     1375               &                 ln_pco2, ln_nea ) 
     1376 
     1377         ENDIF 
     1378  
     1379      ENDIF 
    10041380      
    10051381   END SUBROUTINE dia_obs_init 
    … …  
    10191395      !!               - Sea surface salinity 
    10201396      !!               - Velocity component (U,V) profiles 
     1397      !!               - Sea surface log10(chlorophyll) 
     1398      !!               - Sea surface spm 
     1399      !!               - Sea surface fco2 
     1400      !!               - Sea surface pco2 
    10211401      !! 
    10221402      !! ** Action  :  
    … …  
    10431423         & tmask, umask, vmask                             
    10441424      USE phycst, ONLY : &              ! Physical constants 
    1045          & rday                          
     1425         & rday, & 
     1426         & rt0 
    10461427      USE oce, ONLY : &                 ! Ocean dynamics and tracers variables 
    10471428         & tsn,  &              
    … …  
    10561437         & frld 
    10571438#endif 
     1439#if defined key_hadocc 
     1440      USE trc, ONLY :  &                ! HadOCC chlorophyll, fCO2 and pCO2 
     1441         & HADOCC_CHL, & 
     1442         & HADOCC_FCO2, & 
     1443         & HADOCC_PCO2, & 
     1444         & HADOCC_FILL_FLT 
     1445#elif defined key_medusa && defined key_foam_medusa 
     1446      USE trc, ONLY :  &                ! MEDUSA chlorophyll, fCO2 and pCO2 
     1447         & MEDUSA_CHL, & 
     1448         & MEDUSA_FCO2, & 
     1449         & MEDUSA_PCO2, & 
     1450         & MEDUSA_FILL_FLT 
     1451#elif defined key_fabm 
     1452      USE fabm 
     1453      USE par_fabm 
     1454#endif 
     1455#if defined key_spm 
     1456      USE par_spm, ONLY: &              ! ERSEM/SPM sediments 
     1457         & jp_spm 
     1458      USE trc, ONLY :  & 
     1459         & trn 
     1460#endif 
    10581461      IMPLICIT NONE 
    10591462 
    … …  
    10671470      INTEGER :: jseaiceset             ! sea ice data set loop variable 
    10681471      INTEGER :: jveloset               ! velocity profile data loop variable 
     1472      INTEGER :: jlogchlset             ! logchl data set loop variable 
     1473      INTEGER :: jspmset                ! spm data set loop variable 
     1474      INTEGER :: jfco2set               ! fco2 data set loop variable 
     1475      INTEGER :: jpco2set               ! pco2 data set loop variable 
    10691476      INTEGER :: jvar                   ! Variable number     
    10701477#if ! defined key_lim2 && ! defined key_lim3 
    10711478      REAL(wp), POINTER, DIMENSION(:,:) :: frld    
     1479#endif 
     1480      REAL(wp) :: tiny                  ! small number 
     1481      REAL(wp), DIMENSION(jpi,jpj) :: & 
     1482         logchl                         ! array for log chlorophyll 
     1483      REAL(wp), DIMENSION(jpi,jpj) :: & 
     1484         maskchl                        ! array for special chlorophyll mask 
     1485      REAL(wp), DIMENSION(jpi,jpj) :: & 
     1486         spm                            ! array for spm 
     1487      REAL(wp), DIMENSION(jpi,jpj) :: & 
     1488         fco2                           ! array for fco2 
     1489      REAL(wp), DIMENSION(jpi,jpj) :: & 
     1490         maskfco2                       ! array for special fco2 mask 
     1491      REAL(wp), DIMENSION(jpi,jpj) :: & 
     1492         pco2                           ! array for pco2 
     1493      REAL(wp), DIMENSION(jpi,jpj) :: & 
     1494         maskpco2                       ! array for special pco2 mask 
     1495      INTEGER :: jn                     ! loop index 
     1496#if defined key_fabm 
     1497      REAL(wp), DIMENSION(jpi,jpj,jpk) :: logchl_3d 
     1498      REAL(wp), DIMENSION(jpi,jpj,jpk) :: pco2_3d 
    10721499#endif 
    10731500      CHARACTER(LEN=20) :: datestr=" ",timestr=" " 
    … …  
    11801607      ENDIF 
    11811608 
     1609      IF ( ln_logchl ) THEN 
     1610 
     1611#if defined key_hadocc 
     1612         logchl(:,:) = HADOCC_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC 
     1613#elif defined key_medusa && defined key_foam_medusa 
     1614         logchl(:,:) = MEDUSA_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC 
     1615#elif defined key_fabm 
     1616         logchl_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabmdia_chltot) 
     1617         logchl(:,:) = logchl_3d(:,:,1) 
     1618#else 
     1619         CALL ctl_stop( ' Trying to run logchl observation operator', & 
     1620            &           ' but no biogeochemical model appears to have been defined' ) 
     1621#endif 
     1622 
     1623         maskchl(:,:) = tmask(:,:,1)         ! create a special mask to exclude certain things 
     1624 
     1625         ! Take the log10 where we can, otherwise exclude 
     1626         tiny = 1.0e-20 
     1627         WHERE(logchl(:,:) > tiny .AND. logchl(:,:) /= obfillflt ) 
     1628            logchl(:,:)  = LOG10(logchl(:,:)) 
     1629         ELSEWHERE 
     1630            logchl(:,:)  = obfillflt 
     1631            maskchl(:,:) = 0 
     1632         END WHERE 
     1633 
     1634         DO jlogchlset = 1, nlogchlsets 
     1635             CALL obs_logchl_opt( logchldatqc(jlogchlset),             & 
     1636               &                  kstp, jpi, jpj, nit000, logchl(:,:), & 
     1637               &                  maskchl(:,:), n2dint ) 
     1638         END DO          
     1639      ENDIF  
     1640 
     1641      IF ( ln_spm ) THEN 
     1642#if defined key_spm 
     1643         spm(:,:) = 0.0 
     1644         DO jn = 1, jp_spm 
     1645            spm(:,:) = spm(:,:) + trn(:,:,1,jn)   ! sum SPM sizes 
     1646         END DO 
     1647#else 
     1648         CALL ctl_stop( ' Trying to run spm observation operator', & 
     1649            &           ' but no spm model appears to have been defined' ) 
     1650#endif 
     1651 
     1652         DO jspmset = 1, nspmsets 
     1653             CALL obs_spm_opt( spmdatqc(jspmset),                & 
     1654               &               kstp, jpi, jpj, nit000, spm(:,:), & 
     1655               &               tmask(:,:,1), n2dint ) 
     1656         END DO          
     1657      ENDIF 
     1658 
     1659      IF ( ln_fco2 ) THEN 
     1660         maskfco2(:,:) = tmask(:,:,1)         ! create a special mask to exclude certain things 
     1661#if defined key_hadocc 
     1662         fco2(:,:) = HADOCC_FCO2(:,:)    ! fCO2 from HadOCC 
     1663         IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ).AND.( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN 
     1664            fco2(:,:) = obfillflt 
     1665            maskfco2(:,:) = 0 
     1666            CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', & 
     1667               &           ' on timestep ' // TRIM(STR(kstp)),                              & 
     1668               &           ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' ) 
     1669         ENDIF 
     1670#elif defined key_medusa && defined key_foam_medusa 
     1671         fco2(:,:) = MEDUSA_FCO2(:,:)    ! fCO2 from MEDUSA 
     1672         IF ( ( MINVAL( MEDUSA_FCO2 ) == MEDUSA_FILL_FLT ).AND.( MAXVAL( MEDUSA_FCO2 ) == MEDUSA_FILL_FLT ) ) THEN 
     1673            fco2(:,:) = obfillflt 
     1674            maskfco2(:,:) = 0 
     1675            CALL ctl_warn( ' MEDUSA fCO2 values masked out for observation operator', & 
     1676               &           ' on timestep ' // TRIM(STR(kstp)),                              & 
     1677               &           ' as MEDUSA_FCO2(:,:) == MEDUSA_FILL_FLT' ) 
     1678         ENDIF 
     1679#elif defined key_fabm 
     1680         ! First, get pCO2 from FABM 
     1681         pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc) 
     1682         pco2(:,:) = pco2_3d(:,:,1) 
     1683         ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of: 
     1684         ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems 
     1685         ! and data reduction routines, Deep-Sea Research II, 56: 512-522. 
     1686         ! and 
     1687         ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas, 
     1688         ! Marine Chemistry, 2: 203-215. 
     1689         ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so 
     1690         ! not explicitly included - atmospheric pressure is not necessarily available so this is 
     1691         ! the best assumption. 
     1692         ! Further, the (1-xCO2)^2 term has been neglected. This is common practice 
     1693         ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes) 
     1694         ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal 
     1695         ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway. 
     1696         fco2(:,:) = pco2(:,:) * EXP((-1636.75                                                                               + & 
     1697            &                         12.0408      * (tsn(:,:,1,jp_tem)+rt0)                                                 - & 
     1698            &                         0.0327957    * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)                         + & 
     1699            &                         0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & 
     1700            &                         2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0)))                                        / & 
     1701            &                        (82.0578 * (tsn(:,:,1,jp_tem)+rt0))) 
     1702#else 
     1703         CALL ctl_stop( ' Trying to run fco2 observation operator', & 
     1704            &           ' but no biogeochemical model appears to have been defined' ) 
     1705#endif 
     1706 
     1707         DO jfco2set = 1, nfco2sets 
     1708             CALL obs_fco2_opt( fco2datqc(jfco2set),                      & 
     1709               &                kstp, jpi, jpj, nit000, fco2(:,:), & 
     1710               &                maskfco2(:,:), n2dint ) 
     1711         END DO 
     1712      ENDIF 
     1713 
     1714      IF ( ln_pco2 ) THEN 
     1715         maskpco2(:,:) = tmask(:,:,1)         ! create a special mask to exclude certain things 
     1716#if defined key_hadocc 
     1717         pco2(:,:) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC 
     1718         IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ).AND.( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN 
     1719            pco2(:,:) = obfillflt 
     1720            maskpco2(:,:) = 0 
     1721            CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', & 
     1722               &           ' on timestep ' // TRIM(STR(kstp)),                              & 
     1723               &           ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' ) 
     1724         ENDIF 
     1725#elif defined key_medusa && defined key_foam_medusa 
     1726         pco2(:,:) = MEDUSA_PCO2(:,:)    ! pCO2 from MEDUSA 
     1727         IF ( ( MINVAL( MEDUSA_PCO2 ) == MEDUSA_FILL_FLT ).AND.( MAXVAL( MEDUSA_PCO2 ) == MEDUSA_FILL_FLT ) ) THEN 
     1728            pco2(:,:) = obfillflt 
     1729            maskpco2(:,:) = 0 
     1730            CALL ctl_warn( ' MEDUSA pCO2 values masked out for observation operator', & 
     1731               &           ' on timestep ' // TRIM(STR(kstp)),                              & 
     1732               &           ' as MEDUSA_PCO2(:,:) == MEDUSA_FILL_FLT' ) 
     1733         ENDIF 
     1734#elif defined key_fabm 
     1735         pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc) 
     1736         pco2(:,:) = pco2_3d(:,:,1) 
     1737#else 
     1738         CALL ctl_stop( ' Trying to run pCO2 observation operator', & 
     1739            &           ' but no biogeochemical model appears to have been defined' ) 
     1740#endif 
     1741 
     1742         DO jpco2set = 1, npco2sets 
     1743             CALL obs_pco2_opt( pco2datqc(jpco2set),                      & 
     1744               &                kstp, jpi, jpj, nit000, pco2(:,:), & 
     1745               &                maskpco2(:,:), n2dint ) 
     1746         END DO 
     1747      ENDIF 
     1748 
    11821749#if ! defined key_lim2 && ! defined key_lim3 
    11831750      CALL wrk_dealloc(jpi,jpj,frld)  
    … …  
    12121779      INTEGER :: jsstset                  ! SST data set loop variable 
    12131780      INTEGER :: jseaiceset               ! Sea Ice data set loop variable 
     1781      INTEGER :: jlogchlset               ! logchl data set loop variable 
     1782      INTEGER :: jspmset                  ! spm data set loop variable 
     1783      INTEGER :: jfco2set                 ! fco2 data set loop variable 
     1784      INTEGER :: jpco2set                 ! pco2 data set loop variable 
    12141785      INTEGER :: jset 
    12151786      INTEGER :: jfbini 
    12161787      CHARACTER(LEN=20) :: datestr=" ",timestr=" " 
    1217       CHARACTER(LEN=10) :: cdtmp 
     1788      CHARACTER(LEN=20) :: cdtmp 
    12181789      !----------------------------------------------------------------------- 
    12191790      ! Depending on switches call various observation output routines 
    … …  
    14562027         ENDIF 
    14572028          
     2029      ENDIF 
     2030 
     2031      !  - log10(chlorophyll) 
     2032      IF ( ln_logchl ) THEN 
     2033 
     2034         ! Copy data from logchldatqc to logchldata structures 
     2035         DO jlogchlset = 1, nlogchlsets 
     2036 
     2037            CALL obs_surf_decompress( logchldatqc(jlogchlset), & 
     2038                 &                    logchldata(jlogchlset), .TRUE., numout ) 
     2039 
     2040         END DO 
     2041          
     2042         ! Mark as bad observations with no valid model counterpart due to activities in dia_obs 
     2043         ! Seem to need to set to fill value rather than marking as bad to be effective, so do both 
     2044         DO jlogchlset = 1, nlogchlsets 
     2045            WHERE ( logchldata(jlogchlset)%rmod(:,1) == obfillflt ) 
     2046               logchldata(jlogchlset)%nqc(:)    = 1 
     2047               logchldata(jlogchlset)%robs(:,1) = obfillflt 
     2048            END WHERE 
     2049         END DO 
     2050 
     2051         ! Write the logchl data 
     2052         DO jlogchlset = 1, nlogchlsets 
     2053       
     2054            WRITE(cdtmp,'(A,I2.2)')'logchlfb_',jlogchlset 
     2055            CALL obs_wri_logchl( cdtmp, logchldata(jlogchlset) ) 
     2056 
     2057         END DO 
     2058 
     2059      ENDIF 
     2060 
     2061      !  - spm 
     2062      IF ( ln_spm ) THEN 
     2063 
     2064         ! Copy data from spmdatqc to spmdata structures 
     2065         DO jspmset = 1, nspmsets 
     2066 
     2067            CALL obs_surf_decompress( spmdatqc(jspmset), & 
     2068                 &                    spmdata(jspmset), .TRUE., numout ) 
     2069 
     2070         END DO 
     2071 
     2072         ! Write the spm data 
     2073         DO jspmset = 1, nspmsets 
     2074       
     2075            WRITE(cdtmp,'(A,I2.2)')'spmfb_',jspmset 
     2076            CALL obs_wri_spm( cdtmp, spmdata(jspmset) ) 
     2077 
     2078         END DO 
     2079 
     2080      ENDIF 
     2081 
     2082      !  - fco2 
     2083      IF ( ln_fco2 ) THEN 
     2084 
     2085         ! Copy data from fco2datqc to fco2data structures 
     2086         DO jfco2set = 1, nfco2sets 
     2087 
     2088            CALL obs_surf_decompress( fco2datqc(jfco2set), & 
     2089                 &                    fco2data(jfco2set), .TRUE., numout ) 
     2090 
     2091         END DO 
     2092          
     2093         ! Mark as bad observations with no valid model counterpart due to fCO2 not being in the restart 
     2094         ! Seem to need to set to fill value rather than marking as bad to be effective, so do both 
     2095         DO jfco2set = 1, nfco2sets 
     2096            WHERE ( fco2data(jfco2set)%rmod(:,1) == obfillflt ) 
     2097               fco2data(jfco2set)%nqc(:)    = 1 
     2098               fco2data(jfco2set)%robs(:,1) = obfillflt 
     2099            END WHERE 
     2100         END DO 
     2101 
     2102         ! Write the fco2 data 
     2103         DO jfco2set = 1, nfco2sets 
     2104       
     2105            WRITE(cdtmp,'(A,I2.2)')'fco2fb_',jfco2set 
     2106            CALL obs_wri_fco2( cdtmp, fco2data(jfco2set) ) 
     2107 
     2108         END DO 
     2109 
     2110      ENDIF 
     2111 
     2112      !  - pco2 
     2113      IF ( ln_pco2 ) THEN 
     2114 
     2115         ! Copy data from pco2datqc to pco2data structures 
     2116         DO jpco2set = 1, npco2sets 
     2117 
     2118            CALL obs_surf_decompress( pco2datqc(jpco2set), & 
     2119                 &                    pco2data(jpco2set), .TRUE., numout ) 
     2120 
     2121         END DO 
     2122          
     2123         ! Mark as bad observations with no valid model counterpart due to pco2 not being in the restart 
     2124         ! Seem to need to set to fill value rather than marking as bad to be effective, so do both 
     2125         DO jpco2set = 1, npco2sets 
     2126            WHERE ( pco2data(jpco2set)%rmod(:,1) == obfillflt ) 
     2127               pco2data(jpco2set)%nqc(:)    = 1 
     2128               pco2data(jpco2set)%robs(:,1) = obfillflt 
     2129            END WHERE 
     2130         END DO 
     2131 
     2132         ! Write the pco2 data 
     2133         DO jpco2set = 1, npco2sets 
     2134       
     2135            WRITE(cdtmp,'(A,I2.2)')'pco2fb_',jpco2set 
     2136            CALL obs_wri_pco2( cdtmp, pco2data(jpco2set) ) 
     2137 
     2138         END DO 
     2139 
    14582140      ENDIF 
    14592141 

  • 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 

  • branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90

    r6990 r7713  
    1212   !!   obs_pre_seaice : First level check and screening of sea ice observations 
    1313   !!   obs_pre_vel  : First level check and screening of velocity obs. 
     14   !!   obs_pre_logchl : First level check and screening of logchl obs. 
     15   !!   obs_pre_spm  : First level check and screening of spm obs. 
     16   !!   obs_pre_fco2 : First level check and screening of fco2 obs. 
     17   !!   obs_pre_pco2 : First level check and screening of pco2 obs. 
    1418   !!   obs_scr      : Basic screening of the observations 
    1519   !!   obs_coo_tim  : Compute number of time steps to the observation time 
    … …  
    4549      & obs_pre_seaice, & ! First level check and screening of sea ice data 
    4650      & obs_pre_vel, &     ! First level check and screening of velocity profiles 
     51      & obs_pre_logchl, & ! First level check and screening of logchl data 
     52      & obs_pre_spm, &    ! First level check and screening of spm data 
     53      & obs_pre_fco2, &   ! First level check and screening of fco2 data 
     54      & obs_pre_pco2, &   ! First level check and screening of pco2 data 
    4755      & calc_month_len     ! Calculate the number of days in the months of a year   
    4856 
    … …  
    13761384   END SUBROUTINE obs_pre_vel 
    13771385 
     1386   SUBROUTINE obs_pre_logchl( logchldata, logchldatqc, ld_logchl, ld_nea ) 
     1387      !!---------------------------------------------------------------------- 
     1388      !!                    ***  ROUTINE obs_pre_logchl  *** 
     1389      !! 
     1390      !! ** Purpose : First level check and screening of logchl observations 
     1391      !! 
     1392      !! ** Method  : First level check and screening of logchl observations 
     1393      !! 
     1394      !! ** Action  :  
     1395      !! 
     1396      !! References : 
     1397      !!    
     1398      !! History : 
     1399      !!---------------------------------------------------------------------- 
     1400      !! * Modules used 
     1401      USE domstp              ! Domain: set the time-step 
     1402      USE par_oce             ! Ocean parameters 
     1403      USE dom_oce, ONLY : &   ! Geographical information 
     1404         & glamt,   & 
     1405         & gphit,   & 
     1406         & tmask 
     1407      !! * Arguments 
     1408      TYPE(obs_surf), INTENT(INOUT) :: logchldata     ! Full set of logchl data 
     1409      TYPE(obs_surf), INTENT(INOUT) :: logchldatqc    ! Subset of logchl data not failing screening 
     1410      LOGICAL :: ld_logchl     ! Switch for logchl data 
     1411      LOGICAL :: ld_nea        ! Switch for rejecting observation near land 
     1412      !! * Local declarations 
     1413      INTEGER :: iyea0         ! Initial date 
     1414      INTEGER :: imon0         !  - (year, month, day, hour, minute) 
     1415      INTEGER :: iday0     
     1416      INTEGER :: ihou0     
     1417      INTEGER :: imin0 
     1418      INTEGER :: icycle       ! Current assimilation cycle 
     1419                              ! Counters for observations that 
     1420      INTEGER :: iotdobs      !  - outside time domain 
     1421      INTEGER :: iosdsobs     !  - outside space domain 
     1422      INTEGER :: ilansobs     !  - within a model land cell 
     1423      INTEGER :: inlasobs     !  - close to land 
     1424      INTEGER :: igrdobs      !  - fail the grid search 
     1425      INTEGER :: ibdysobs     !  - close to open boundary 
     1426                              ! Global counters for observations that 
     1427      INTEGER :: iotdobsmpp   !  - outside time domain 
     1428      INTEGER :: iosdsobsmpp  !  - outside space domain 
     1429      INTEGER :: ilansobsmpp  !  - within a model land cell 
     1430      INTEGER :: inlasobsmpp  !  - close to land 
     1431      INTEGER :: igrdobsmpp   !  - fail the grid search 
     1432      INTEGER :: ibdysobsmpp  !  - close to open boundary 
     1433      LOGICAL, DIMENSION(:), ALLOCATABLE :: &  
     1434         & llvalid            ! data selection 
     1435      INTEGER :: jobs         ! Obs. loop variable 
     1436      INTEGER :: jstp         ! Time loop variable 
     1437      INTEGER :: inrc         ! Time index variable 
     1438      INTEGER :: irec         ! Record index 
     1439 
     1440      IF (lwp) WRITE(numout,*)'obs_pre_logchl : Preparing the logchl observations...' 
     1441 
     1442      ! Initial date initialization (year, month, day, hour, minute) 
     1443      iyea0 =   ndate0 / 10000 
     1444      imon0 = ( ndate0 - iyea0 * 10000 ) / 100 
     1445      iday0 =   ndate0 - iyea0 * 10000 - imon0 * 100 
     1446      ihou0 = 0 
     1447      imin0 = 0 
     1448 
     1449      icycle = no     ! Assimilation cycle 
     1450 
     1451      ! Diagnostics counters for various failures. 
     1452 
     1453      iotdobs  = 0 
     1454      igrdobs  = 0 
     1455      iosdsobs = 0 
     1456      ilansobs = 0 
     1457      inlasobs = 0 
     1458      ibdysobs = 0 
     1459 
     1460      ! ----------------------------------------------------------------------- 
     1461      ! Find time coordinate for logchl data 
     1462      ! ----------------------------------------------------------------------- 
     1463 
     1464      CALL obs_coo_tim( icycle, & 
     1465         &              iyea0,   imon0,   iday0,   ihou0,   imin0,      & 
     1466         &              logchldata%nsurf,   logchldata%nyea, logchldata%nmon, & 
     1467         &              logchldata%nday,    logchldata%nhou, logchldata%nmin, & 
     1468         &              logchldata%nqc,     logchldata%mstp, iotdobs        ) 
     1469      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 
     1470      ! ----------------------------------------------------------------------- 
     1471      ! Check for logchl data failing the grid search 
     1472      ! ----------------------------------------------------------------------- 
     1473 
     1474      CALL obs_coo_grd( logchldata%nsurf,   logchldata%mi, logchldata%mj, & 
     1475         &              logchldata%nqc,     igrdobs                         ) 
     1476      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) 
     1477 
     1478      ! ----------------------------------------------------------------------- 
     1479      ! Check for land points.  
     1480      ! ----------------------------------------------------------------------- 
     1481 
     1482      CALL obs_coo_spc_2d( logchldata%nsurf,                 & 
     1483         &                 jpi,             jpj,             & 
     1484         &                 logchldata%mi,   logchldata%mj,   &  
     1485         &                 logchldata%rlam, logchldata%rphi, & 
     1486         &                 glamt,           gphit,           & 
     1487         &                 tmask(:,:,1),    logchldata%nqc,  & 
     1488         &                 iosdsobs,        ilansobs,        & 
     1489         &                 inlasobs,        ld_nea,          & 
     1490         &                 ibdysobs,        ln_bound_reject  )  
     1491          
     1492      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 
     1493      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 
     1494      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 
     1495      CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp ) 
     1496 
     1497      ! ----------------------------------------------------------------------- 
     1498      ! Copy useful data from the logchldata data structure to 
     1499      ! the logchldatqc data structure  
     1500      ! ----------------------------------------------------------------------- 
     1501 
     1502      ! Allocate the selection arrays 
     1503 
     1504      ALLOCATE( llvalid(logchldata%nsurf) ) 
     1505       
     1506      ! We want all data which has qc flags <= 0 
     1507 
     1508      llvalid(:)  = ( logchldata%nqc(:)  <= 10 ) 
     1509 
     1510      ! The actual copying 
     1511 
     1512      CALL obs_surf_compress( logchldata,     logchldatqc,       .TRUE.,  numout, & 
     1513         &                    lvalid=llvalid ) 
     1514 
     1515      ! Dellocate the selection arrays 
     1516      DEALLOCATE( llvalid ) 
     1517 
     1518      ! ----------------------------------------------------------------------- 
     1519      ! Print information about what observations are left after qc 
     1520      ! ----------------------------------------------------------------------- 
     1521 
     1522      ! Update the total observation counter array 
     1523       
     1524      IF(lwp) THEN 
     1525         WRITE(numout,*) 
     1526         WRITE(numout,*) 'obs_pre_logchl :' 
     1527         WRITE(numout,*) '~~~~~~~~~~~' 
     1528         WRITE(numout,*) 
     1529         WRITE(numout,*) ' logchl data outside time domain                  = ', & 
     1530            &            iotdobsmpp 
     1531         WRITE(numout,*) ' Remaining logchl data that failed grid search    = ', & 
     1532            &            igrdobsmpp 
     1533         WRITE(numout,*) ' Remaining logchl data outside space domain       = ', & 
     1534            &            iosdsobsmpp 
     1535         WRITE(numout,*) ' Remaining logchl data at land points             = ', & 
     1536            &            ilansobsmpp 
     1537         IF (ld_nea) THEN 
     1538            WRITE(numout,*) ' Remaining logchl data near land points (removed) = ', & 
     1539               &            inlasobsmpp 
     1540         ELSE 
     1541            WRITE(numout,*) ' Remaining logchl data near land points (kept)    = ', & 
     1542               &            inlasobsmpp 
     1543         ENDIF 
     1544         WRITE(numout,*) ' Remaining logchl data near open boundary (removed) = ', & 
     1545           &            ibdysobsmpp 
     1546         WRITE(numout,*) ' logchl data accepted                             = ', & 
     1547            &            logchldatqc%nsurfmpp 
     1548 
     1549         WRITE(numout,*) 
     1550         WRITE(numout,*) ' Number of observations per time step :' 
     1551         WRITE(numout,*) 
     1552         WRITE(numout,1997) 
     1553         WRITE(numout,1998) 
     1554      ENDIF 
     1555       
     1556      DO jobs = 1, logchldatqc%nsurf 
     1557         inrc = logchldatqc%mstp(jobs) + 2 - nit000 
     1558         logchldatqc%nsstp(inrc)  = logchldatqc%nsstp(inrc) + 1 
     1559      END DO 
     1560       
     1561      CALL obs_mpp_sum_integers( logchldatqc%nsstp, logchldatqc%nsstpmpp, & 
     1562         &                       nitend - nit000 + 2 ) 
     1563 
     1564      IF ( lwp ) THEN 
     1565         DO jstp = nit000 - 1, nitend 
     1566            inrc = jstp - nit000 + 2 
     1567            WRITE(numout,1999) jstp, logchldatqc%nsstpmpp(inrc) 
     1568         END DO 
     1569      ENDIF 
     1570 
     15711997  FORMAT(10X,'Time step',5X,'logchl data') 
     15721998  FORMAT(10X,'---------',5X,'------------') 
     15731999  FORMAT(10X,I9,5X,I17) 
     1574       
     1575   END SUBROUTINE obs_pre_logchl 
     1576 
     1577   SUBROUTINE obs_pre_spm( spmdata, spmdatqc, ld_spm, ld_nea ) 
     1578      !!---------------------------------------------------------------------- 
     1579      !!                    ***  ROUTINE obs_pre_spm  *** 
     1580      !! 
     1581      !! ** Purpose : First level check and screening of spm observations 
     1582      !! 
     1583      !! ** Method  : First level check and screening of spm observations 
     1584      !! 
     1585      !! ** Action  :  
     1586      !! 
     1587      !! References : 
     1588      !!    
     1589      !! History : 
     1590      !!---------------------------------------------------------------------- 
     1591      !! * Modules used 
     1592      USE domstp              ! Domain: set the time-step 
     1593      USE par_oce             ! Ocean parameters 
     1594      USE dom_oce, ONLY : &   ! Geographical information 
     1595         & glamt,   & 
     1596         & gphit,   & 
     1597         & tmask 
     1598      !! * Arguments 
     1599      TYPE(obs_surf), INTENT(INOUT) :: spmdata     ! Full set of spm data 
     1600      TYPE(obs_surf), INTENT(INOUT) :: spmdatqc    ! Subset of spm data not failing screening 
     1601      LOGICAL :: ld_spm     ! Switch for spm data 
     1602      LOGICAL :: ld_nea        ! Switch for rejecting observation near land 
     1603      !! * Local declarations 
     1604      INTEGER :: iyea0         ! Initial date 
     1605      INTEGER :: imon0         !  - (year, month, day, hour, minute) 
     1606      INTEGER :: iday0     
     1607      INTEGER :: ihou0     
     1608      INTEGER :: imin0 
     1609      INTEGER :: icycle       ! Current assimilation cycle 
     1610                              ! Counters for observations that 
     1611      INTEGER :: iotdobs      !  - outside time domain 
     1612      INTEGER :: iosdsobs     !  - outside space domain 
     1613      INTEGER :: ilansobs     !  - within a model land cell 
     1614      INTEGER :: inlasobs     !  - close to land 
     1615      INTEGER :: igrdobs      !  - fail the grid search 
     1616      INTEGER :: ibdysobs     !  - close to open boundary 
     1617                              ! Global counters for observations that 
     1618      INTEGER :: iotdobsmpp   !  - outside time domain 
     1619      INTEGER :: iosdsobsmpp  !  - outside space domain 
     1620      INTEGER :: ilansobsmpp  !  - within a model land cell 
     1621      INTEGER :: inlasobsmpp  !  - close to land 
     1622      INTEGER :: igrdobsmpp   !  - fail the grid search 
     1623      INTEGER :: ibdysobsmpp  !  - close to open boundary 
     1624      LOGICAL, DIMENSION(:), ALLOCATABLE :: &  
     1625         & llvalid            ! data selection 
     1626      INTEGER :: jobs         ! Obs. loop variable 
     1627      INTEGER :: jstp         ! Time loop variable 
     1628      INTEGER :: inrc         ! Time index variable 
     1629      INTEGER :: irec         ! Record index 
     1630 
     1631      IF (lwp) WRITE(numout,*)'obs_pre_spm : Preparing the spm observations...' 
     1632 
     1633      ! Initial date initialization (year, month, day, hour, minute) 
     1634      iyea0 =   ndate0 / 10000 
     1635      imon0 = ( ndate0 - iyea0 * 10000 ) / 100 
     1636      iday0 =   ndate0 - iyea0 * 10000 - imon0 * 100 
     1637      ihou0 = 0 
     1638      imin0 = 0 
     1639 
     1640      icycle = no     ! Assimilation cycle 
     1641 
     1642      ! Diagnostics counters for various failures. 
     1643 
     1644      iotdobs  = 0 
     1645      igrdobs  = 0 
     1646      iosdsobs = 0 
     1647      ilansobs = 0 
     1648      inlasobs = 0 
     1649      ibdysobs = 0 
     1650 
     1651      ! ----------------------------------------------------------------------- 
     1652      ! Find time coordinate for spm data 
     1653      ! ----------------------------------------------------------------------- 
     1654 
     1655      CALL obs_coo_tim( icycle, & 
     1656         &              iyea0,   imon0,   iday0,   ihou0,   imin0,      & 
     1657         &              spmdata%nsurf,   spmdata%nyea, spmdata%nmon, & 
     1658         &              spmdata%nday,    spmdata%nhou, spmdata%nmin, & 
     1659         &              spmdata%nqc,     spmdata%mstp, iotdobs        ) 
     1660      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 
     1661      ! ----------------------------------------------------------------------- 
     1662      ! Check for spm data failing the grid search 
     1663      ! ----------------------------------------------------------------------- 
     1664 
     1665      CALL obs_coo_grd( spmdata%nsurf,   spmdata%mi, spmdata%mj, & 
     1666         &              spmdata%nqc,     igrdobs                         ) 
     1667      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) 
     1668 
     1669      ! ----------------------------------------------------------------------- 
     1670      ! Check for land points.  
     1671      ! ----------------------------------------------------------------------- 
     1672 
     1673      CALL obs_coo_spc_2d( spmdata%nsurf,                 & 
     1674         &                 jpi,             jpj,             & 
     1675         &                 spmdata%mi,   spmdata%mj,   &  
     1676         &                 spmdata%rlam, spmdata%rphi, & 
     1677         &                 glamt,           gphit,           & 
     1678         &                 tmask(:,:,1),    spmdata%nqc,  & 
     1679         &                 iosdsobs,        ilansobs,        & 
     1680         &                 inlasobs,        ld_nea,          & 
     1681         &                 ibdysobs,        ln_bound_reject  )  
     1682          
     1683      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 
     1684      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 
     1685      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 
     1686      CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp ) 
     1687 
     1688      ! ----------------------------------------------------------------------- 
     1689      ! Copy useful data from the spmdata data structure to 
     1690      ! the spmdatqc data structure  
     1691      ! ----------------------------------------------------------------------- 
     1692 
     1693      ! Allocate the selection arrays 
     1694 
     1695      ALLOCATE( llvalid(spmdata%nsurf) ) 
     1696       
     1697      ! We want all data which has qc flags <= 0 
     1698 
     1699      llvalid(:)  = ( spmdata%nqc(:)  <= 10 ) 
     1700 
     1701      ! The actual copying 
     1702 
     1703      CALL obs_surf_compress( spmdata,     spmdatqc,       .TRUE.,  numout, & 
     1704         &                    lvalid=llvalid ) 
     1705 
     1706      ! Dellocate the selection arrays 
     1707      DEALLOCATE( llvalid ) 
     1708 
     1709      ! ----------------------------------------------------------------------- 
     1710      ! Print information about what observations are left after qc 
     1711      ! ----------------------------------------------------------------------- 
     1712 
     1713      ! Update the total observation counter array 
     1714       
     1715      IF(lwp) THEN 
     1716         WRITE(numout,*) 
     1717         WRITE(numout,*) 'obs_pre_spm :' 
     1718         WRITE(numout,*) '~~~~~~~~~~~' 
     1719         WRITE(numout,*) 
     1720         WRITE(numout,*) ' spm data outside time domain                  = ', & 
     1721            &            iotdobsmpp 
     1722         WRITE(numout,*) ' Remaining spm data that failed grid search    = ', & 
     1723            &            igrdobsmpp 
     1724         WRITE(numout,*) ' Remaining spm data outside space domain       = ', & 
     1725            &            iosdsobsmpp 
     1726         WRITE(numout,*) ' Remaining spm data at land points             = ', & 
     1727            &            ilansobsmpp 
     1728         IF (ld_nea) THEN 
     1729            WRITE(numout,*) ' Remaining spm data near land points (removed) = ', & 
     1730               &            inlasobsmpp 
     1731         ELSE 
     1732            WRITE(numout,*) ' Remaining spm data near land points (kept)    = ', & 
     1733               &            inlasobsmpp 
     1734         ENDIF 
     1735         WRITE(numout,*) ' Remaining spm data near open boundary (removed) = ', & 
     1736            &            ibdysobsmpp 
     1737         WRITE(numout,*) ' spm data accepted                             = ', & 
     1738            &            spmdatqc%nsurfmpp 
     1739 
     1740         WRITE(numout,*) 
     1741         WRITE(numout,*) ' Number of observations per time step :' 
     1742         WRITE(numout,*) 
     1743         WRITE(numout,1997) 
     1744         WRITE(numout,1998) 
     1745      ENDIF 
     1746       
     1747      DO jobs = 1, spmdatqc%nsurf 
     1748         inrc = spmdatqc%mstp(jobs) + 2 - nit000 
     1749         spmdatqc%nsstp(inrc)  = spmdatqc%nsstp(inrc) + 1 
     1750      END DO 
     1751       
     1752      CALL obs_mpp_sum_integers( spmdatqc%nsstp, spmdatqc%nsstpmpp, & 
     1753         &                       nitend - nit000 + 2 ) 
     1754 
     1755      IF ( lwp ) THEN 
     1756         DO jstp = nit000 - 1, nitend 
     1757            inrc = jstp - nit000 + 2 
     1758            WRITE(numout,1999) jstp, spmdatqc%nsstpmpp(inrc) 
     1759         END DO 
     1760      ENDIF 
     1761 
     17621997  FORMAT(10X,'Time step',5X,'spm data') 
     17631998  FORMAT(10X,'---------',5X,'------------') 
     17641999  FORMAT(10X,I9,5X,I17) 
     1765       
     1766   END SUBROUTINE obs_pre_spm 
     1767 
     1768   SUBROUTINE obs_pre_fco2( fco2data, fco2datqc, ld_fco2, ld_nea ) 
     1769      !!---------------------------------------------------------------------- 
     1770      !!                    ***  ROUTINE obs_pre_fco2  *** 
     1771      !! 
     1772      !! ** Purpose : First level check and screening of fco2 observations 
     1773      !! 
     1774      !! ** Method  : First level check and screening of fco2 observations 
     1775      !! 
     1776      !! ** Action  :  
     1777      !! 
     1778      !! References : 
     1779      !!    
     1780      !! History : 
     1781      !!---------------------------------------------------------------------- 
     1782      !! * Modules used 
     1783      USE domstp              ! Domain: set the time-step 
     1784      USE par_oce             ! Ocean parameters 
     1785      USE dom_oce, ONLY : &   ! Geographical information 
     1786         & glamt,   & 
     1787         & gphit,   & 
     1788         & tmask 
     1789      !! * Arguments 
     1790      TYPE(obs_surf), INTENT(INOUT) :: fco2data     ! Full set of fco2 data 
     1791      TYPE(obs_surf), INTENT(INOUT) :: fco2datqc    ! Subset of fco2 data not failing screening 
     1792      LOGICAL :: ld_fco2     ! Switch for fco2 data 
     1793      LOGICAL :: ld_nea        ! Switch for rejecting observation near land 
     1794      !! * Local declarations 
     1795      INTEGER :: iyea0         ! Initial date 
     1796      INTEGER :: imon0         !  - (year, month, day, hour, minute) 
     1797      INTEGER :: iday0     
     1798      INTEGER :: ihou0     
     1799      INTEGER :: imin0 
     1800      INTEGER :: icycle       ! Current assimilation cycle 
     1801                              ! Counters for observations that 
     1802      INTEGER :: iotdobs      !  - outside time domain 
     1803      INTEGER :: iosdsobs     !  - outside space domain 
     1804      INTEGER :: ilansobs     !  - within a model land cell 
     1805      INTEGER :: inlasobs     !  - close to land 
     1806      INTEGER :: igrdobs      !  - fail the grid search 
     1807      INTEGER :: ibdysobs     !  - close to open boundary 
     1808                              ! Global counters for observations that 
     1809      INTEGER :: iotdobsmpp   !  - outside time domain 
     1810      INTEGER :: iosdsobsmpp  !  - outside space domain 
     1811      INTEGER :: ilansobsmpp  !  - within a model land cell 
     1812      INTEGER :: inlasobsmpp  !  - close to land 
     1813      INTEGER :: igrdobsmpp   !  - fail the grid search 
     1814      INTEGER :: ibdysobsmpp  !  - close to open boundary 
     1815      LOGICAL, DIMENSION(:), ALLOCATABLE :: &  
     1816         & llvalid            ! data selection 
     1817      INTEGER :: jobs         ! Obs. loop variable 
     1818      INTEGER :: jstp         ! Time loop variable 
     1819      INTEGER :: inrc         ! Time index variable 
     1820      INTEGER :: irec         ! Record index 
     1821 
     1822      IF (lwp) WRITE(numout,*)'obs_pre_fco2 : Preparing the fco2 observations...' 
     1823 
     1824      ! Initial date initialization (year, month, day, hour, minute) 
     1825      iyea0 =   ndate0 / 10000 
     1826      imon0 = ( ndate0 - iyea0 * 10000 ) / 100 
     1827      iday0 =   ndate0 - iyea0 * 10000 - imon0 * 100 
     1828      ihou0 = 0 
     1829      imin0 = 0 
     1830 
     1831      icycle = no     ! Assimilation cycle 
     1832 
     1833      ! Diagnostics counters for various failures. 
     1834 
     1835      iotdobs  = 0 
     1836      igrdobs  = 0 
     1837      iosdsobs = 0 
     1838      ilansobs = 0 
     1839      inlasobs = 0 
     1840      ibdysobs = 0 
     1841 
     1842      ! ----------------------------------------------------------------------- 
     1843      ! Find time coordinate for fco2 data 
     1844      ! ----------------------------------------------------------------------- 
     1845 
     1846      CALL obs_coo_tim( icycle, & 
     1847         &              iyea0,   imon0,   iday0,   ihou0,   imin0,      & 
     1848         &              fco2data%nsurf,   fco2data%nyea, fco2data%nmon, & 
     1849         &              fco2data%nday,    fco2data%nhou, fco2data%nmin, & 
     1850         &              fco2data%nqc,     fco2data%mstp, iotdobs        ) 
     1851      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 
     1852      ! ----------------------------------------------------------------------- 
     1853      ! Check for fco2 data failing the grid search 
     1854      ! ----------------------------------------------------------------------- 
     1855 
     1856      CALL obs_coo_grd( fco2data%nsurf,   fco2data%mi, fco2data%mj, & 
     1857         &              fco2data%nqc,     igrdobs                         ) 
     1858      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) 
     1859 
     1860      ! ----------------------------------------------------------------------- 
     1861      ! Check for land points.  
     1862      ! ----------------------------------------------------------------------- 
     1863 
     1864      CALL obs_coo_spc_2d( fco2data%nsurf,                 & 
     1865         &                 jpi,             jpj,             & 
     1866         &                 fco2data%mi,   fco2data%mj,   &  
     1867         &                 fco2data%rlam, fco2data%rphi, & 
     1868         &                 glamt,           gphit,           & 
     1869         &                 tmask(:,:,1),    fco2data%nqc,  & 
     1870         &                 iosdsobs,        ilansobs,        & 
     1871         &                 inlasobs,        ld_nea,          & 
     1872         &                 ibdysobs,        ln_bound_reject  )  
     1873          
     1874      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 
     1875      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 
     1876      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 
     1877      CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp ) 
     1878 
     1879      ! ----------------------------------------------------------------------- 
     1880      ! Copy useful data from the fco2data data structure to 
     1881      ! the fco2datqc data structure  
     1882      ! ----------------------------------------------------------------------- 
     1883 
     1884      ! Allocate the selection arrays 
     1885 
     1886      ALLOCATE( llvalid(fco2data%nsurf) ) 
     1887       
     1888      ! We want all data which has qc flags <= 0 
     1889 
     1890      llvalid(:)  = ( fco2data%nqc(:)  <= 10 ) 
     1891 
     1892      ! The actual copying 
     1893 
     1894      CALL obs_surf_compress( fco2data,     fco2datqc,       .TRUE.,  numout, & 
     1895         &                    lvalid=llvalid ) 
     1896 
     1897      ! Dellocate the selection arrays 
     1898      DEALLOCATE( llvalid ) 
     1899 
     1900      ! ----------------------------------------------------------------------- 
     1901      ! Print information about what observations are left after qc 
     1902      ! ----------------------------------------------------------------------- 
     1903 
     1904      ! Update the total observation counter array 
     1905       
     1906      IF(lwp) THEN 
     1907         WRITE(numout,*) 
     1908         WRITE(numout,*) 'obs_pre_fco2 :' 
     1909         WRITE(numout,*) '~~~~~~~~~~~' 
     1910         WRITE(numout,*) 
     1911         WRITE(numout,*) ' fco2 data outside time domain                  = ', & 
     1912            &            iotdobsmpp 
     1913         WRITE(numout,*) ' Remaining fco2 data that failed grid search    = ', & 
     1914            &            igrdobsmpp 
     1915         WRITE(numout,*) ' Remaining fco2 data outside space domain       = ', & 
     1916            &            iosdsobsmpp 
     1917         WRITE(numout,*) ' Remaining fco2 data at land points             = ', & 
     1918            &            ilansobsmpp 
     1919         IF (ld_nea) THEN 
     1920            WRITE(numout,*) ' Remaining fco2 data near land points (removed) = ', & 
     1921               &            inlasobsmpp 
     1922         ELSE 
     1923            WRITE(numout,*) ' Remaining fco2 data near land points (kept)    = ', & 
     1924               &            inlasobsmpp 
     1925         ENDIF 
     1926         WRITE(numout,*) ' Remaining fco2 data near open boundary (removed) = ', & 
     1927           &            ibdysobsmpp 
     1928         WRITE(numout,*) ' fco2 data accepted                             = ', & 
     1929            &            fco2datqc%nsurfmpp 
     1930 
     1931         WRITE(numout,*) 
     1932         WRITE(numout,*) ' Number of observations per time step :' 
     1933         WRITE(numout,*) 
     1934         WRITE(numout,1997) 
     1935         WRITE(numout,1998) 
     1936      ENDIF 
     1937       
     1938      DO jobs = 1, fco2datqc%nsurf 
     1939         inrc = fco2datqc%mstp(jobs) + 2 - nit000 
     1940         fco2datqc%nsstp(inrc)  = fco2datqc%nsstp(inrc) + 1 
     1941      END DO 
     1942       
     1943      CALL obs_mpp_sum_integers( fco2datqc%nsstp, fco2datqc%nsstpmpp, & 
     1944         &                       nitend - nit000 + 2 ) 
     1945 
     1946      IF ( lwp ) THEN 
     1947         DO jstp = nit000 - 1, nitend 
     1948            inrc = jstp - nit000 + 2 
     1949            WRITE(numout,1999) jstp, fco2datqc%nsstpmpp(inrc) 
     1950         END DO 
     1951      ENDIF 
     1952 
     19531997  FORMAT(10X,'Time step',5X,'fco2 data') 
     19541998  FORMAT(10X,'---------',5X,'------------') 
     19551999  FORMAT(10X,I9,5X,I17) 
     1956       
     1957   END SUBROUTINE obs_pre_fco2 
     1958 
     1959   SUBROUTINE obs_pre_pco2( pco2data, pco2datqc, ld_pco2, ld_nea ) 
     1960      !!---------------------------------------------------------------------- 
     1961      !!                    ***  ROUTINE obs_pre_pco2  *** 
     1962      !! 
     1963      !! ** Purpose : First level check and screening of pco2 observations 
     1964      !! 
     1965      !! ** Method  : First level check and screening of pco2 observations 
     1966      !! 
     1967      !! ** Action  :  
     1968      !! 
     1969      !! References : 
     1970      !!    
     1971      !! History : 
     1972      !!---------------------------------------------------------------------- 
     1973      !! * Modules used 
     1974      USE domstp              ! Domain: set the time-step 
     1975      USE par_oce             ! Ocean parameters 
     1976      USE dom_oce, ONLY : &   ! Geographical information 
     1977         & glamt,   & 
     1978         & gphit,   & 
     1979         & tmask 
     1980      !! * Arguments 
     1981      TYPE(obs_surf), INTENT(INOUT) :: pco2data     ! Full set of pco2 data 
     1982      TYPE(obs_surf), INTENT(INOUT) :: pco2datqc    ! Subset of pco2 data not failing screening 
     1983      LOGICAL :: ld_pco2     ! Switch for pco2 data 
     1984      LOGICAL :: ld_nea        ! Switch for rejecting observation near land 
     1985      !! * Local declarations 
     1986      INTEGER :: iyea0         ! Initial date 
     1987      INTEGER :: imon0         !  - (year, month, day, hour, minute) 
     1988      INTEGER :: iday0     
     1989      INTEGER :: ihou0     
     1990      INTEGER :: imin0 
     1991      INTEGER :: icycle       ! Current assimilation cycle 
     1992                              ! Counters for observations that 
     1993      INTEGER :: iotdobs      !  - outside time domain 
     1994      INTEGER :: iosdsobs     !  - outside space domain 
     1995      INTEGER :: ilansobs     !  - within a model land cell 
     1996      INTEGER :: inlasobs     !  - close to land 
     1997      INTEGER :: igrdobs      !  - fail the grid search 
     1998      INTEGER :: ibdysobs     !  - close to open boundary 
     1999                              ! Global counters for observations that 
     2000      INTEGER :: iotdobsmpp   !  - outside time domain 
     2001      INTEGER :: iosdsobsmpp  !  - outside space domain 
     2002      INTEGER :: ilansobsmpp  !  - within a model land cell 
     2003      INTEGER :: inlasobsmpp  !  - close to land 
     2004      INTEGER :: igrdobsmpp   !  - fail the grid search 
     2005      INTEGER :: ibdysobsmpp  !  - close to open boundary 
     2006      LOGICAL, DIMENSION(:), ALLOCATABLE :: &  
     2007         & llvalid            ! data selection 
     2008      INTEGER :: jobs         ! Obs. loop variable 
     2009      INTEGER :: jstp         ! Time loop variable 
     2010      INTEGER :: inrc         ! Time index variable 
     2011      INTEGER :: irec         ! Record index 
     2012 
     2013      IF (lwp) WRITE(numout,*)'obs_pre_pco2 : Preparing the pco2 observations...' 
     2014 
     2015      ! Initial date initialization (year, month, day, hour, minute) 
     2016      iyea0 =   ndate0 / 10000 
     2017      imon0 = ( ndate0 - iyea0 * 10000 ) / 100 
     2018      iday0 =   ndate0 - iyea0 * 10000 - imon0 * 100 
     2019      ihou0 = 0 
     2020      imin0 = 0 
     2021 
     2022      icycle = no     ! Assimilation cycle 
     2023 
     2024      ! Diagnostics counters for various failures. 
     2025 
     2026      iotdobs  = 0 
     2027      igrdobs  = 0 
     2028      iosdsobs = 0 
     2029      ilansobs = 0 
     2030      inlasobs = 0 
     2031      ibdysobs = 0 
     2032 
     2033      ! ----------------------------------------------------------------------- 
     2034      ! Find time coordinate for pco2 data 
     2035      ! ----------------------------------------------------------------------- 
     2036 
     2037      CALL obs_coo_tim( icycle, & 
     2038         &              iyea0,   imon0,   iday0,   ihou0,   imin0,      & 
     2039         &              pco2data%nsurf,   pco2data%nyea, pco2data%nmon, & 
     2040         &              pco2data%nday,    pco2data%nhou, pco2data%nmin, & 
     2041         &              pco2data%nqc,     pco2data%mstp, iotdobs        ) 
     2042      CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp ) 
     2043      ! ----------------------------------------------------------------------- 
     2044      ! Check for pco2 data failing the grid search 
     2045      ! ----------------------------------------------------------------------- 
     2046 
     2047      CALL obs_coo_grd( pco2data%nsurf,   pco2data%mi, pco2data%mj, & 
     2048         &              pco2data%nqc,     igrdobs                         ) 
     2049      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) 
     2050 
     2051      ! ----------------------------------------------------------------------- 
     2052      ! Check for land points.  
     2053      ! ----------------------------------------------------------------------- 
     2054 
     2055      CALL obs_coo_spc_2d( pco2data%nsurf,                 & 
     2056         &                 jpi,             jpj,             & 
     2057         &                 pco2data%mi,   pco2data%mj,   &  
     2058         &                 pco2data%rlam, pco2data%rphi, & 
     2059         &                 glamt,           gphit,           & 
     2060         &                 tmask(:,:,1),    pco2data%nqc,  & 
     2061         &                 iosdsobs,        ilansobs,        & 
     2062         &                 inlasobs,        ld_nea,          & 
     2063         &                 ibdysobs,        ln_bound_reject  )  
     2064          
     2065      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 
     2066      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 
     2067      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 
     2068      CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp ) 
     2069 
     2070      ! ----------------------------------------------------------------------- 
     2071      ! Copy useful data from the pco2data data structure to 
     2072      ! the pco2datqc data structure  
     2073      ! ----------------------------------------------------------------------- 
     2074 
     2075      ! Allocate the selection arrays 
     2076 
     2077      ALLOCATE( llvalid(pco2data%nsurf) ) 
     2078       
     2079      ! We want all data which has qc flags <= 0 
     2080 
     2081      llvalid(:)  = ( pco2data%nqc(:)  <= 10 ) 
     2082 
     2083      ! The actual copying 
     2084 
     2085      CALL obs_surf_compress( pco2data,     pco2datqc,       .TRUE.,  numout, & 
     2086         &                    lvalid=llvalid ) 
     2087 
     2088      ! Dellocate the selection arrays 
     2089      DEALLOCATE( llvalid ) 
     2090 
     2091      ! ----------------------------------------------------------------------- 
     2092      ! Print information about what observations are left after qc 
     2093      ! ----------------------------------------------------------------------- 
     2094 
     2095      ! Update the total observation counter array 
     2096       
     2097      IF(lwp) THEN 
     2098         WRITE(numout,*) 
     2099         WRITE(numout,*) 'obs_pre_pco2 :' 
     2100         WRITE(numout,*) '~~~~~~~~~~~' 
     2101         WRITE(numout,*) 
     2102         WRITE(numout,*) ' pco2 data outside time domain                  = ', & 
     2103            &            iotdobsmpp 
     2104         WRITE(numout,*) ' Remaining pco2 data that failed grid search    = ', & 
     2105            &            igrdobsmpp 
     2106         WRITE(numout,*) ' Remaining pco2 data outside space domain       = ', & 
     2107            &            iosdsobsmpp 
     2108         WRITE(numout,*) ' Remaining pco2 data at land points             = ', & 
     2109            &            ilansobsmpp 
     2110         IF (ld_nea) THEN 
     2111            WRITE(numout,*) ' Remaining pco2 data near land points (removed) = ', & 
     2112               &            inlasobsmpp 
     2113         ELSE 
     2114            WRITE(numout,*) ' Remaining pco2 data near land points (kept)    = ', & 
     2115               &            inlasobsmpp 
     2116         ENDIF 
     2117         WRITE(numout,*) ' Remaining pco2 data near open boundary (removed) = ', & 
     2118           &            ibdysobsmpp 
     2119         WRITE(numout,*) ' pco2 data accepted                             = ', & 
     2120            &            pco2datqc%nsurfmpp 
     2121 
     2122         WRITE(numout,*) 
     2123         WRITE(numout,*) ' Number of observations per time step :' 
     2124         WRITE(numout,*) 
     2125         WRITE(numout,1997) 
     2126         WRITE(numout,1998) 
     2127      ENDIF 
     2128       
     2129      DO jobs = 1, pco2datqc%nsurf 
     2130         inrc = pco2datqc%mstp(jobs) + 2 - nit000 
     2131         pco2datqc%nsstp(inrc)  = pco2datqc%nsstp(inrc) + 1 
     2132      END DO 
     2133       
     2134      CALL obs_mpp_sum_integers( pco2datqc%nsstp, pco2datqc%nsstpmpp, & 
     2135         &                       nitend - nit000 + 2 ) 
     2136 
     2137      IF ( lwp ) THEN 
     2138         DO jstp = nit000 - 1, nitend 
     2139            inrc = jstp - nit000 + 2 
     2140            WRITE(numout,1999) jstp, pco2datqc%nsstpmpp(inrc) 
     2141         END DO 
     2142      ENDIF 
     2143 
     21441997  FORMAT(10X,'Time step',5X,'pco2 data') 
     21451998  FORMAT(10X,'---------',5X,'------------') 
     21461999  FORMAT(10X,I9,5X,I17) 
     2147       
     2148   END SUBROUTINE obs_pre_pco2 
     2149 
    13782150   SUBROUTINE obs_coo_tim( kcycle, & 
    13792151      &                    kyea0,   kmon0,   kday0,   khou0,   kmin0,     & 

  • branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90

    r5838 r7713  
    1111   !!   obs_wri_seaice: Write seaice observation related diagnostics 
    1212   !!   obs_wri_vel   : Write velocity observation diagnostics in NetCDF format 
     13   !!   obs_wri_logchl: Write logchl observation related diagnostics 
     14   !!   obs_wri_spm   : Write spm observation related diagnostics 
     15   !!   obs_wri_fco2  : Write fco2 observation related diagnostics 
     16   !!   obs_wri_pco2  : Write fco2 observation related diagnostics 
    1317   !!   obs_wri_stats : Print basic statistics on the data being written out 
    1418   !!---------------------------------------------------------------------- 
    … …  
    4549      &   obs_wri_seaice, & ! Write seaice observation related diagnostics 
    4650      &   obs_wri_vel, &    ! Write velocity observation related diagnostics 
     51      &   obs_wri_logchl, & ! Write logchl observation related diagnostics 
     52      &   obs_wri_spm, &    ! Write spm observation related diagnostics 
     53      &   obs_wri_fco2, &   ! Write fco2 observation related diagnostics 
     54      &   obs_wri_pco2, &   ! Write pco2 observation related diagnostics 
    4755      &   obswriinfo 
    4856    
    … …  
    930938   END SUBROUTINE obs_wri_vel 
    931939 
     940   SUBROUTINE obs_wri_logchl( cprefix, logchldata, padd, pext ) 
     941      !!----------------------------------------------------------------------- 
     942      !! 
     943      !!                     *** ROUTINE obs_wri_logchl  *** 
     944      !! 
     945      !! ** Purpose : Write logchl observation diagnostics 
     946      !!              related  
     947      !! 
     948      !! ** Method  : NetCDF 
     949      !!  
     950      !! ** Action  : 
     951      !! 
     952      !!----------------------------------------------------------------------- 
     953 
     954      !! * Modules used 
     955      IMPLICIT NONE 
     956 
     957      !! * Arguments 
     958      CHARACTER(LEN=*), INTENT(IN) :: cprefix       ! Prefix for output files 
     959      TYPE(obs_surf), INTENT(INOUT) :: logchldata   ! Full set of logchl 
     960      TYPE(obswriinfo), OPTIONAL :: padd            ! Additional info for each variable 
     961      TYPE(obswriinfo), OPTIONAL :: pext            ! Extra info 
     962 
     963      !! * Local declarations  
     964      TYPE(obfbdata) :: fbdata 
     965      CHARACTER(LEN=40) :: cfname             ! netCDF filename 
     966      CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_logchl' 
     967      INTEGER :: jo 
     968      INTEGER :: ja 
     969      INTEGER :: je 
     970      INTEGER :: nadd 
     971      INTEGER :: next 
     972 
     973      IF ( PRESENT( padd ) ) THEN 
     974         nadd = padd%inum 
     975      ELSE 
     976         nadd = 0 
     977      ENDIF 
     978 
     979      IF ( PRESENT( pext ) ) THEN 
     980         next = pext%inum 
     981      ELSE 
     982         next = 0 
     983      ENDIF 
     984 
     985      CALL init_obfbdata( fbdata ) 
     986 
     987      CALL alloc_obfbdata( fbdata, 1, logchldata%nsurf, 1, & 
     988         &                 1 + nadd, next, .TRUE. ) 
     989 
     990      fbdata%cname(1)      = 'LOGCHL' 
     991      fbdata%coblong(1)    = 'logchl concentration' 
     992      fbdata%cobunit(1)    = 'mg/m3' 
     993      DO je = 1, next 
     994         fbdata%cextname(je) = pext%cdname(je) 
     995         fbdata%cextlong(je) = pext%cdlong(je,1) 
     996         fbdata%cextunit(je) = pext%cdunit(je,1) 
     997      END DO 
     998      fbdata%caddname(1)   = 'Hx' 
     999      fbdata%caddlong(1,1) = 'Model interpolated LOGCHL' 
     1000      fbdata%caddunit(1,1) = 'mg/m3' 
     1001      fbdata%cgrid(1)      = 'T' 
     1002      DO ja = 1, nadd 
     1003         fbdata%caddname(1+ja) = padd%cdname(ja) 
     1004         fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
     1005         fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
     1006      END DO 
     1007 
     1008      WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cprefix), nproc 
     1009 
     1010      IF(lwp) THEN 
     1011         WRITE(numout,*) 
     1012         WRITE(numout,*)'obs_wri_logchl :' 
     1013         WRITE(numout,*)'~~~~~~~~~~~~~~~~' 
     1014         WRITE(numout,*)'Writing logchl feedback file : ',TRIM(cfname) 
     1015      ENDIF 
     1016 
     1017      ! Transform obs_prof data structure into obfbdata structure 
     1018      fbdata%cdjuldref = '19500101000000' 
     1019      DO jo = 1, logchldata%nsurf 
     1020         fbdata%plam(jo)      = logchldata%rlam(jo) 
     1021         fbdata%pphi(jo)      = logchldata%rphi(jo) 
     1022         WRITE(fbdata%cdtyp(jo),'(I4)') logchldata%ntyp(jo) 
     1023         fbdata%ivqc(jo,:)    = 0 
     1024         fbdata%ivqcf(:,jo,:) = 0 
     1025         IF ( logchldata%nqc(jo) > 10 ) THEN 
     1026            fbdata%ioqc(jo)    = 4 
     1027            fbdata%ioqcf(1,jo) = 0 
     1028            fbdata%ioqcf(2,jo) = logchldata%nqc(jo) - 10 
     1029         ELSE 
     1030            fbdata%ioqc(jo)    = MAX(logchldata%nqc(jo),1) 
     1031            fbdata%ioqcf(:,jo) = 0 
     1032         ENDIF 
     1033         fbdata%ipqc(jo)      = 0 
     1034         fbdata%ipqcf(:,jo)   = 0 
     1035         fbdata%itqc(jo)      = 0 
     1036         fbdata%itqcf(:,jo)   = 0 
     1037         fbdata%cdwmo(jo)     = '' 
     1038         fbdata%kindex(jo)    = logchldata%nsfil(jo) 
     1039         IF (ln_grid_global) THEN 
     1040            fbdata%iobsi(jo,1) = logchldata%mi(jo) 
     1041            fbdata%iobsj(jo,1) = logchldata%mj(jo) 
     1042         ELSE 
     1043            fbdata%iobsi(jo,1) = mig(logchldata%mi(jo)) 
     1044            fbdata%iobsj(jo,1) = mjg(logchldata%mj(jo)) 
     1045         ENDIF 
     1046         CALL greg2jul( 0, & 
     1047            &           logchldata%nmin(jo), & 
     1048            &           logchldata%nhou(jo), & 
     1049            &           logchldata%nday(jo), & 
     1050            &           logchldata%nmon(jo), & 
     1051            &           logchldata%nyea(jo), & 
     1052            &           fbdata%ptim(jo),   & 
     1053            &           krefdate = 19500101 ) 
     1054         fbdata%padd(1,jo,1,1) = logchldata%rmod(jo,1) 
     1055         fbdata%pob(1,jo,1)    = logchldata%robs(jo,1) 
     1056         fbdata%pdep(1,jo)     = 0.0 
     1057         fbdata%idqc(1,jo)     = 0 
     1058         fbdata%idqcf(:,1,jo)  = 0 
     1059         IF ( logchldata%nqc(jo) > 10 ) THEN 
     1060            fbdata%ivlqc(1,jo,1) = 4 
     1061            fbdata%ivlqcf(1,1,jo,1) = 0 
     1062            fbdata%ivlqcf(2,1,jo,1) = logchldata%nqc(jo) - 10 
     1063         ELSE 
     1064            fbdata%ivlqc(1,jo,1) = MAX(logchldata%nqc(jo),1) 
     1065            fbdata%ivlqcf(:,1,jo,1) = 0 
     1066         ENDIF 
     1067         fbdata%iobsk(1,jo,1)  = 0 
     1068         DO ja = 1, nadd 
     1069            fbdata%padd(1,jo,1+ja,1) = & 
     1070               & logchldata%rext(jo,padd%ipoint(ja)) 
     1071         END DO 
     1072         DO je = 1, next 
     1073            fbdata%pext(1,jo,je) = & 
     1074               & logchldata%rext(jo,pext%ipoint(je)) 
     1075         END DO 
     1076 
     1077      END DO 
     1078 
     1079      ! Write the obfbdata structure 
     1080      CALL write_obfbdata( cfname, fbdata ) 
     1081       
     1082      ! Output some basic statistics 
     1083      CALL obs_wri_stats( fbdata ) 
     1084 
     1085      CALL dealloc_obfbdata( fbdata ) 
     1086 
     1087   END SUBROUTINE obs_wri_logchl 
     1088 
     1089   SUBROUTINE obs_wri_spm( cprefix, spmdata, padd, pext ) 
     1090      !!----------------------------------------------------------------------- 
     1091      !! 
     1092      !!                     *** ROUTINE obs_wri_spm  *** 
     1093      !! 
     1094      !! ** Purpose : Write spm observation diagnostics 
     1095      !!              related  
     1096      !! 
     1097      !! ** Method  : NetCDF 
     1098      !!  
     1099      !! ** Action  : 
     1100      !! 
     1101      !!----------------------------------------------------------------------- 
     1102 
     1103      !! * Modules used 
     1104      IMPLICIT NONE 
     1105 
     1106      !! * Arguments 
     1107      CHARACTER(LEN=*), INTENT(IN) :: cprefix       ! Prefix for output files 
     1108      TYPE(obs_surf), INTENT(INOUT) :: spmdata   ! Full set of spm 
     1109      TYPE(obswriinfo), OPTIONAL :: padd            ! Additional info for each variable 
     1110      TYPE(obswriinfo), OPTIONAL :: pext            ! Extra info 
     1111 
     1112      !! * Local declarations  
     1113      TYPE(obfbdata) :: fbdata 
     1114      CHARACTER(LEN=40) :: cfname             ! netCDF filename 
     1115      CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_spm' 
     1116      INTEGER :: jo 
     1117      INTEGER :: ja 
     1118      INTEGER :: je 
     1119      INTEGER :: nadd 
     1120      INTEGER :: next 
     1121 
     1122      IF ( PRESENT( padd ) ) THEN 
     1123         nadd = padd%inum 
     1124      ELSE 
     1125         nadd = 0 
     1126      ENDIF 
     1127 
     1128      IF ( PRESENT( pext ) ) THEN 
     1129         next = pext%inum 
     1130      ELSE 
     1131         next = 0 
     1132      ENDIF 
     1133 
     1134      CALL init_obfbdata( fbdata ) 
     1135 
     1136      CALL alloc_obfbdata( fbdata, 1, spmdata%nsurf, 1, & 
     1137         &                 1 + nadd, next, .TRUE. ) 
     1138 
     1139      fbdata%cname(1)      = 'spm' 
     1140      fbdata%coblong(1)    = 'spm' 
     1141      fbdata%cobunit(1)    = 'g/m3' 
     1142      DO je = 1, next 
     1143         fbdata%cextname(je) = pext%cdname(je) 
     1144         fbdata%cextlong(je) = pext%cdlong(je,1) 
     1145         fbdata%cextunit(je) = pext%cdunit(je,1) 
     1146      END DO 
     1147      fbdata%caddname(1)   = 'Hx' 
     1148      fbdata%caddlong(1,1) = 'Model interpolated spm' 
     1149      fbdata%caddunit(1,1) = 'g/m3' 
     1150      fbdata%cgrid(1)      = 'T' 
     1151      DO ja = 1, nadd 
     1152         fbdata%caddname(1+ja) = padd%cdname(ja) 
     1153         fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
     1154         fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
     1155      END DO 
     1156 
     1157      WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cprefix), nproc 
     1158 
     1159      IF(lwp) THEN 
     1160         WRITE(numout,*) 
     1161         WRITE(numout,*)'obs_wri_spm :' 
     1162         WRITE(numout,*)'~~~~~~~~~~~~~~~~' 
     1163         WRITE(numout,*)'Writing spm feedback file : ',TRIM(cfname) 
     1164      ENDIF 
     1165 
     1166      ! Transform obs_prof data structure into obfbdata structure 
     1167      fbdata%cdjuldref = '19500101000000' 
     1168      DO jo = 1, spmdata%nsurf 
     1169         fbdata%plam(jo)      = spmdata%rlam(jo) 
     1170         fbdata%pphi(jo)      = spmdata%rphi(jo) 
     1171         WRITE(fbdata%cdtyp(jo),'(I4)') spmdata%ntyp(jo) 
     1172         fbdata%ivqc(jo,:)    = 0 
     1173         fbdata%ivqcf(:,jo,:) = 0 
     1174         IF ( spmdata%nqc(jo) > 10 ) THEN 
     1175            fbdata%ioqc(jo)    = 4 
     1176            fbdata%ioqcf(1,jo) = 0 
     1177            fbdata%ioqcf(2,jo) = spmdata%nqc(jo) - 10 
     1178         ELSE 
     1179            fbdata%ioqc(jo)    = MAX(spmdata%nqc(jo),1) 
     1180            fbdata%ioqcf(:,jo) = 0 
     1181         ENDIF 
     1182         fbdata%ipqc(jo)      = 0 
     1183         fbdata%ipqcf(:,jo)   = 0 
     1184         fbdata%itqc(jo)      = 0 
     1185         fbdata%itqcf(:,jo)   = 0 
     1186         fbdata%cdwmo(jo)     = '' 
     1187         fbdata%kindex(jo)    = spmdata%nsfil(jo) 
     1188         IF (ln_grid_global) THEN 
     1189            fbdata%iobsi(jo,1) = spmdata%mi(jo) 
     1190            fbdata%iobsj(jo,1) = spmdata%mj(jo) 
     1191         ELSE 
     1192            fbdata%iobsi(jo,1) = mig(spmdata%mi(jo)) 
     1193            fbdata%iobsj(jo,1) = mjg(spmdata%mj(jo)) 
     1194         ENDIF 
     1195         CALL greg2jul( 0, & 
     1196            &           spmdata%nmin(jo), & 
     1197            &           spmdata%nhou(jo), & 
     1198            &           spmdata%nday(jo), & 
     1199            &           spmdata%nmon(jo), & 
     1200            &           spmdata%nyea(jo), & 
     1201            &           fbdata%ptim(jo),   & 
     1202            &           krefdate = 19500101 ) 
     1203         fbdata%padd(1,jo,1,1) = spmdata%rmod(jo,1) 
     1204         fbdata%pob(1,jo,1)    = spmdata%robs(jo,1) 
     1205         fbdata%pdep(1,jo)     = 0.0 
     1206         fbdata%idqc(1,jo)     = 0 
     1207         fbdata%idqcf(:,1,jo)  = 0 
     1208         IF ( spmdata%nqc(jo) > 10 ) THEN 
     1209            fbdata%ivlqc(1,jo,1) = 4 
     1210            fbdata%ivlqcf(1,1,jo,1) = 0 
     1211            fbdata%ivlqcf(2,1,jo,1) = spmdata%nqc(jo) - 10 
     1212         ELSE 
     1213            fbdata%ivlqc(1,jo,1) = MAX(spmdata%nqc(jo),1) 
     1214            fbdata%ivlqcf(:,1,jo,1) = 0 
     1215         ENDIF 
     1216         fbdata%iobsk(1,jo,1)  = 0 
     1217         DO ja = 1, nadd 
     1218            fbdata%padd(1,jo,1+ja,1) = & 
     1219               & spmdata%rext(jo,padd%ipoint(ja)) 
     1220         END DO 
     1221         DO je = 1, next 
     1222            fbdata%pext(1,jo,je) = & 
     1223               & spmdata%rext(jo,pext%ipoint(je)) 
     1224         END DO 
     1225 
     1226      END DO 
     1227 
     1228      ! Write the obfbdata structure 
     1229      CALL write_obfbdata( cfname, fbdata ) 
     1230       
     1231      ! Output some basic statistics 
     1232      CALL obs_wri_stats( fbdata ) 
     1233 
     1234      CALL dealloc_obfbdata( fbdata ) 
     1235 
     1236   END SUBROUTINE obs_wri_spm 
     1237 
     1238   SUBROUTINE obs_wri_fco2( cprefix, fco2data, padd, pext ) 
     1239      !!----------------------------------------------------------------------- 
     1240      !! 
     1241      !!                     *** ROUTINE obs_wri_fco2  *** 
     1242      !! 
     1243      !! ** Purpose : Write fco2 observation diagnostics 
     1244      !!              related  
     1245      !! 
     1246      !! ** Method  : NetCDF 
     1247      !!  
     1248      !! ** Action  : 
     1249      !! 
     1250      !!----------------------------------------------------------------------- 
     1251 
     1252      !! * Modules used 
     1253      IMPLICIT NONE 
     1254 
     1255      !! * Arguments 
     1256      CHARACTER(LEN=*), INTENT(IN) :: cprefix       ! Prefix for output files 
     1257      TYPE(obs_surf), INTENT(INOUT) :: fco2data   ! Full set of fco2 
     1258      TYPE(obswriinfo), OPTIONAL :: padd            ! Additional info for each variable 
     1259      TYPE(obswriinfo), OPTIONAL :: pext            ! Extra info 
     1260 
     1261      !! * Local declarations  
     1262      TYPE(obfbdata) :: fbdata 
     1263      CHARACTER(LEN=40) :: cfname             ! netCDF filename 
     1264      CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_fco2' 
     1265      INTEGER :: jo 
     1266      INTEGER :: ja 
     1267      INTEGER :: je 
     1268      INTEGER :: nadd 
     1269      INTEGER :: next 
     1270 
     1271      IF ( PRESENT( padd ) ) THEN 
     1272         nadd = padd%inum 
     1273      ELSE 
     1274         nadd = 0 
     1275      ENDIF 
     1276 
     1277      IF ( PRESENT( pext ) ) THEN 
     1278         next = pext%inum 
     1279      ELSE 
     1280         next = 0 
     1281      ENDIF 
     1282 
     1283      CALL init_obfbdata( fbdata ) 
     1284 
     1285      CALL alloc_obfbdata( fbdata, 1, fco2data%nsurf, 1, & 
     1286         &                 1 + nadd, next, .TRUE. ) 
     1287 
     1288      fbdata%cname(1)      = 'fco2' 
     1289      fbdata%coblong(1)    = 'fco2' 
     1290      fbdata%cobunit(1)    = 'uatm' 
     1291      DO je = 1, next 
     1292         fbdata%cextname(je) = pext%cdname(je) 
     1293         fbdata%cextlong(je) = pext%cdlong(je,1) 
     1294         fbdata%cextunit(je) = pext%cdunit(je,1) 
     1295      END DO 
     1296      fbdata%caddname(1)   = 'Hx' 
     1297      fbdata%caddlong(1,1) = 'Model interpolated fco2' 
     1298      fbdata%caddunit(1,1) = 'uatm' 
     1299      fbdata%cgrid(1)      = 'T' 
     1300      DO ja = 1, nadd 
     1301         fbdata%caddname(1+ja) = padd%cdname(ja) 
     1302         fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
     1303         fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
     1304      END DO 
     1305 
     1306      WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cprefix), nproc 
     1307 
     1308      IF(lwp) THEN 
     1309         WRITE(numout,*) 
     1310         WRITE(numout,*)'obs_wri_fco2 :' 
     1311         WRITE(numout,*)'~~~~~~~~~~~~~~~~' 
     1312         WRITE(numout,*)'Writing fco2 feedback file : ',TRIM(cfname) 
     1313      ENDIF 
     1314 
     1315      ! Transform obs_prof data structure into obfbdata structure 
     1316      fbdata%cdjuldref = '19500101000000' 
     1317      DO jo = 1, fco2data%nsurf 
     1318         fbdata%plam(jo)      = fco2data%rlam(jo) 
     1319         fbdata%pphi(jo)      = fco2data%rphi(jo) 
     1320         WRITE(fbdata%cdtyp(jo),'(I4)') fco2data%ntyp(jo) 
     1321         fbdata%ivqc(jo,:)    = 0 
     1322         fbdata%ivqcf(:,jo,:) = 0 
     1323         IF ( fco2data%nqc(jo) > 10 ) THEN 
     1324            fbdata%ioqc(jo)    = 4 
     1325            fbdata%ioqcf(1,jo) = 0 
     1326            fbdata%ioqcf(2,jo) = fco2data%nqc(jo) - 10 
     1327         ELSE 
     1328            fbdata%ioqc(jo)    = MAX(fco2data%nqc(jo),1) 
     1329            fbdata%ioqcf(:,jo) = 0 
     1330         ENDIF 
     1331         fbdata%ipqc(jo)      = 0 
     1332         fbdata%ipqcf(:,jo)   = 0 
     1333         fbdata%itqc(jo)      = 0 
     1334         fbdata%itqcf(:,jo)   = 0 
     1335         fbdata%cdwmo(jo)     = '' 
     1336         fbdata%kindex(jo)    = fco2data%nsfil(jo) 
     1337         IF (ln_grid_global) THEN 
     1338            fbdata%iobsi(jo,1) = fco2data%mi(jo) 
     1339            fbdata%iobsj(jo,1) = fco2data%mj(jo) 
     1340         ELSE 
     1341            fbdata%iobsi(jo,1) = mig(fco2data%mi(jo)) 
     1342            fbdata%iobsj(jo,1) = mjg(fco2data%mj(jo)) 
     1343         ENDIF 
     1344         CALL greg2jul( 0, & 
     1345            &           fco2data%nmin(jo), & 
     1346            &           fco2data%nhou(jo), & 
     1347            &           fco2data%nday(jo), & 
     1348            &           fco2data%nmon(jo), & 
     1349            &           fco2data%nyea(jo), & 
     1350            &           fbdata%ptim(jo),   & 
     1351            &           krefdate = 19500101 ) 
     1352         fbdata%padd(1,jo,1,1) = fco2data%rmod(jo,1) 
     1353         fbdata%pob(1,jo,1)    = fco2data%robs(jo,1) 
     1354         fbdata%pdep(1,jo)     = 0.0 
     1355         fbdata%idqc(1,jo)     = 0 
     1356         fbdata%idqcf(:,1,jo)  = 0 
     1357         IF ( fco2data%nqc(jo) > 10 ) THEN 
     1358            fbdata%ivlqc(1,jo,1) = 4 
     1359            fbdata%ivlqcf(1,1,jo,1) = 0 
     1360            fbdata%ivlqcf(2,1,jo,1) = fco2data%nqc(jo) - 10 
     1361         ELSE 
     1362            fbdata%ivlqc(1,jo,1) = MAX(fco2data%nqc(jo),1) 
     1363            fbdata%ivlqcf(:,1,jo,1) = 0 
     1364         ENDIF 
     1365         fbdata%iobsk(1,jo,1)  = 0 
     1366         DO ja = 1, nadd 
     1367            fbdata%padd(1,jo,1+ja,1) = & 
     1368               & fco2data%rext(jo,padd%ipoint(ja)) 
     1369         END DO 
     1370         DO je = 1, next 
     1371            fbdata%pext(1,jo,je) = & 
     1372               & fco2data%rext(jo,pext%ipoint(je)) 
     1373         END DO 
     1374 
     1375      END DO 
     1376 
     1377      ! Write the obfbdata structure 
     1378      CALL write_obfbdata( cfname, fbdata ) 
     1379       
     1380      ! Output some basic statistics 
     1381      CALL obs_wri_stats( fbdata ) 
     1382 
     1383      CALL dealloc_obfbdata( fbdata ) 
     1384 
     1385   END SUBROUTINE obs_wri_fco2 
     1386 
     1387   SUBROUTINE obs_wri_pco2( cprefix, pco2data, padd, pext ) 
     1388      !!----------------------------------------------------------------------- 
     1389      !! 
     1390      !!                     *** ROUTINE obs_wri_pco2  *** 
     1391      !! 
     1392      !! ** Purpose : Write pco2 observation diagnostics 
     1393      !!              related  
     1394      !! 
     1395      !! ** Method  : NetCDF 
     1396      !!  
     1397      !! ** Action  : 
     1398      !! 
     1399      !!----------------------------------------------------------------------- 
     1400 
     1401      !! * Modules used 
     1402      IMPLICIT NONE 
     1403 
     1404      !! * Arguments 
     1405      CHARACTER(LEN=*), INTENT(IN) :: cprefix       ! Prefix for output files 
     1406      TYPE(obs_surf), INTENT(INOUT) :: pco2data   ! Full set of pco2 
     1407      TYPE(obswriinfo), OPTIONAL :: padd            ! Additional info for each variable 
     1408      TYPE(obswriinfo), OPTIONAL :: pext            ! Extra info 
     1409 
     1410      !! * Local declarations  
     1411      TYPE(obfbdata) :: fbdata 
     1412      CHARACTER(LEN=40) :: cfname             ! netCDF filename 
     1413      CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_pco2' 
     1414      INTEGER :: jo 
     1415      INTEGER :: ja 
     1416      INTEGER :: je 
     1417      INTEGER :: nadd 
     1418      INTEGER :: next 
     1419 
     1420      IF ( PRESENT( padd ) ) THEN 
     1421         nadd = padd%inum 
     1422      ELSE 
     1423         nadd = 0 
     1424      ENDIF 
     1425 
     1426      IF ( PRESENT( pext ) ) THEN 
     1427         next = pext%inum 
     1428      ELSE 
     1429         next = 0 
     1430      ENDIF 
     1431 
     1432      CALL init_obfbdata( fbdata ) 
     1433 
     1434      CALL alloc_obfbdata( fbdata, 1, pco2data%nsurf, 1, & 
     1435         &                 1 + nadd, next, .TRUE. ) 
     1436 
     1437      fbdata%cname(1)      = 'pco2' 
     1438      fbdata%coblong(1)    = 'pco2' 
     1439      fbdata%cobunit(1)    = 'uatm' 
     1440      DO je = 1, next 
     1441         fbdata%cextname(je) = pext%cdname(je) 
     1442         fbdata%cextlong(je) = pext%cdlong(je,1) 
     1443         fbdata%cextunit(je) = pext%cdunit(je,1) 
     1444      END DO 
     1445      fbdata%caddname(1)   = 'Hx' 
     1446      fbdata%caddlong(1,1) = 'Model interpolated pco2' 
     1447      fbdata%caddunit(1,1) = 'uatm' 
     1448      fbdata%cgrid(1)      = 'T' 
     1449      DO ja = 1, nadd 
     1450         fbdata%caddname(1+ja) = padd%cdname(ja) 
     1451         fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
     1452         fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
     1453      END DO 
     1454 
     1455      WRITE(cfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(cprefix), nproc 
     1456 
     1457      IF(lwp) THEN 
     1458         WRITE(numout,*) 
     1459         WRITE(numout,*)'obs_wri_pco2 :' 
     1460         WRITE(numout,*)'~~~~~~~~~~~~~~~~' 
     1461         WRITE(numout,*)'Writing pco2 feedback file : ',TRIM(cfname) 
     1462      ENDIF 
     1463 
     1464      ! Transform obs_prof data structure into obfbdata structure 
     1465      fbdata%cdjuldref = '19500101000000' 
     1466      DO jo = 1, pco2data%nsurf 
     1467         fbdata%plam(jo)      = pco2data%rlam(jo) 
     1468         fbdata%pphi(jo)      = pco2data%rphi(jo) 
     1469         WRITE(fbdata%cdtyp(jo),'(I4)') pco2data%ntyp(jo) 
     1470         fbdata%ivqc(jo,:)    = 0 
     1471         fbdata%ivqcf(:,jo,:) = 0 
     1472         IF ( pco2data%nqc(jo) > 10 ) THEN 
     1473            fbdata%ioqc(jo)    = 4 
     1474            fbdata%ioqcf(1,jo) = 0 
     1475            fbdata%ioqcf(2,jo) = pco2data%nqc(jo) - 10 
     1476         ELSE 
     1477            fbdata%ioqc(jo)    = MAX(pco2data%nqc(jo),1) 
     1478            fbdata%ioqcf(:,jo) = 0 
     1479         ENDIF 
     1480         fbdata%ipqc(jo)      = 0 
     1481         fbdata%ipqcf(:,jo)   = 0 
     1482         fbdata%itqc(jo)      = 0 
     1483         fbdata%itqcf(:,jo)   = 0 
     1484         fbdata%cdwmo(jo)     = '' 
     1485         fbdata%kindex(jo)    = pco2data%nsfil(jo) 
     1486         IF (ln_grid_global) THEN 
     1487            fbdata%iobsi(jo,1) = pco2data%mi(jo) 
     1488            fbdata%iobsj(jo,1) = pco2data%mj(jo) 
     1489         ELSE 
     1490            fbdata%iobsi(jo,1) = mig(pco2data%mi(jo)) 
     1491            fbdata%iobsj(jo,1) = mjg(pco2data%mj(jo)) 
     1492         ENDIF 
     1493         CALL greg2jul( 0, & 
     1494            &           pco2data%nmin(jo), & 
     1495            &           pco2data%nhou(jo), & 
     1496            &           pco2data%nday(jo), & 
     1497            &           pco2data%nmon(jo), & 
     1498            &           pco2data%nyea(jo), & 
     1499            &           fbdata%ptim(jo),   & 
     1500            &           krefdate = 19500101 ) 
     1501         fbdata%padd(1,jo,1,1) = pco2data%rmod(jo,1) 
     1502         fbdata%pob(1,jo,1)    = pco2data%robs(jo,1) 
     1503         fbdata%pdep(1,jo)     = 0.0 
     1504         fbdata%idqc(1,jo)     = 0 
     1505         fbdata%idqcf(:,1,jo)  = 0 
     1506         IF ( pco2data%nqc(jo) > 10 ) THEN 
     1507            fbdata%ivlqc(1,jo,1) = 4 
     1508            fbdata%ivlqcf(1,1,jo,1) = 0 
     1509            fbdata%ivlqcf(2,1,jo,1) = pco2data%nqc(jo) - 10 
     1510         ELSE 
     1511            fbdata%ivlqc(1,jo,1) = MAX(pco2data%nqc(jo),1) 
     1512            fbdata%ivlqcf(:,1,jo,1) = 0 
     1513         ENDIF 
     1514         fbdata%iobsk(1,jo,1)  = 0 
     1515         DO ja = 1, nadd 
     1516            fbdata%padd(1,jo,1+ja,1) = & 
     1517               & pco2data%rext(jo,padd%ipoint(ja)) 
     1518         END DO 
     1519         DO je = 1, next 
     1520            fbdata%pext(1,jo,je) = & 
     1521               & pco2data%rext(jo,pext%ipoint(je)) 
     1522         END DO 
     1523 
     1524      END DO 
     1525 
     1526      ! Write the obfbdata structure 
     1527      CALL write_obfbdata( cfname, fbdata ) 
     1528       
     1529      ! Output some basic statistics 
     1530      CALL obs_wri_stats( fbdata ) 
     1531 
     1532      CALL dealloc_obfbdata( fbdata ) 
     1533 
     1534   END SUBROUTINE obs_wri_pco2 
     1535 
    9321536   SUBROUTINE obs_wri_stats( fbdata ) 
    9331537      !!----------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.