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

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

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

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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/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 
Note: See TracChangeset for help on using the changeset viewer.