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 13393 for branches – NEMO

Changeset 13393 for branches


Ignore:
Timestamp:
2020-08-13T16:44:28+02:00 (4 years ago)
Author:
mattmartin
Message:

Include surface velocity observation operator. See ticket ​https://code.metoffice.gov.uk/trac/utils/ticket/376.

Location:
branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/CONFIG/SHARED/namelist_ref

    r12610 r13393  
    12011201   ln_sit     = .false.             ! Logical switch for Sea Ice Thickness observations 
    12021202   ln_vel3d   = .false.             ! Logical switch for velocity observations 
    1203    ln_sss     = .false.             ! Logical swithc for SSS observations 
     1203   ln_sss     = .false.             ! Logical switch for SSS observations 
     1204   ln_ssv     = .false.             ! Logical switch for SSV (surface velocity) observations    
    12041205   ln_slchltot = .false.            ! Logical switch for surface total              log10(chlorophyll) obs 
    12051206   ln_slchldia = .false.            ! Logical switch for surface diatom             log10(chlorophyll) obs 
    … …  
    12401241   ln_sst_fp_indegs = .true. 
    12411242   ln_sss_fp_indegs = .true. 
     1243   ln_ssv_fp_indegs = .true.    
    12421244   ln_sic_fp_indegs = .true. 
    12431245   ln_sit_fp_indegs = .true. 
    … …  
    12501252   cn_velfbfiles = 'vel_01.nc'           ! Velocity feedback input observation file names 
    12511253   cn_sssfbfiles = 'sss_01.nc'           ! SSS feedback input observation file names 
     1254   cn_ssvfbfiles = 'ssv_01.nc'           ! SSV feedback input observation file names    
    12521255   cn_slchltotfbfiles = 'slchltot_01.nc' ! Surface total              log10(chlorophyll) obs file names 
    12531256   cn_slchldiafbfiles = 'slchldia_01.nc' ! Surface diatom             log10(chlorophyll) obs file names 
    … …  
    12861289   rn_sss_avglamscl = 0.                 ! E/W diameter of SSS observation footprint (metres/degrees) 
    12871290   rn_sss_avgphiscl = 0.                 ! N/S diameter of SSS observation footprint (metres/degrees) 
     1291   rn_ssv_avglamscl = 0.                 ! E/W diameter of SSV observation footprint (metres/degrees) 
     1292   rn_ssv_avgphiscl = 0.                 ! N/S diameter of SSV observation footprint (metres/degrees) 
    12881293   rn_sic_avglamscl = 0.                 ! E/W diameter of SIC observation footprint (metres/degrees) 
    12891294   rn_sic_avgphiscl = 0.                 ! N/S diameter of SIC observation footprint (metres/degrees) 
    … …  
    12951300   nn_2dint_sst = -1                     ! Horizontal interpolation method for SST 
    12961301   nn_2dint_sss = -1                     ! Horizontal interpolation method for SSS 
     1302   nn_2dint_ssv = -1                     ! Horizontal interpolation method for SSV    
    12971303   nn_2dint_sic = -1                     ! Horizontal interpolation method for SIC 
    12981304   nn_2dint_sit = -1                     ! Horizontal interpolation method for SIT 

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

    r12610 r13393  
    5151   LOGICAL :: ln_sstnight          !: Logical switch for night mean SST obs 
    5252   LOGICAL :: ln_default_fp_indegs !: T=> Default obs footprint size specified in degrees, F=> in metres 
    53    LOGICAL :: ln_sla_fp_indegs     !: T=>     SLA obs footprint size specified in degrees, F=> in metres 
    54    LOGICAL :: ln_sst_fp_indegs     !: T=>     SST obs footprint size specified in degrees, F=> in metres 
    55    LOGICAL :: ln_sss_fp_indegs     !: T=>     SSS obs footprint size specified in degrees, F=> in metres 
     53   LOGICAL :: ln_sla_fp_indegs     !: T=> SLA obs footprint size specified in degrees, F=> in metres 
     54   LOGICAL :: ln_sst_fp_indegs     !: T=> SST obs footprint size specified in degrees, F=> in metres 
     55   LOGICAL :: ln_sss_fp_indegs     !: T=> SSS obs footprint size specified in degrees, F=> in metres 
     56   LOGICAL :: ln_ssv_fp_indegs     !: T=> SSV obs footprint size specified in degrees, F=> in metres    
    5657   LOGICAL :: ln_sic_fp_indegs     !: T=> SIC obs footprint size specified in degrees, F=> in metres 
    5758   LOGICAL :: ln_sit_fp_indegs     !: T=> SIT obs footprint size specified in degrees, F=> in metres 
    … …  
    6768   REAL(wp) :: rn_sss_avglamscl     !: E/W diameter of SSS observation footprint 
    6869   REAL(wp) :: rn_sss_avgphiscl     !: N/S diameter of SSS observation footprint 
     70   REAL(wp) :: rn_ssv_avglamscl     !: E/W diameter of SSV observation footprint 
     71   REAL(wp) :: rn_ssv_avgphiscl     !: N/S diameter of SSV observation footprint 
    6972   REAL(wp) :: rn_sic_avglamscl     !: E/W diameter of SIC observation footprint 
    7073   REAL(wp) :: rn_sic_avgphiscl     !: N/S diameter of SIC observation footprint 
    … …  
    8083   INTEGER :: nn_2dint_sst     !: SST horizontal interpolation method (-1 = default) 
    8184   INTEGER :: nn_2dint_sss     !: SSS horizontal interpolation method (-1 = default) 
     85   INTEGER :: nn_2dint_ssv     !: SSV horizontal interpolation method (-1 = default)    
    8286   INTEGER :: nn_2dint_sic     !: SIC horizontal interpolation method (-1 = default) 
    8387   INTEGER :: nn_2dint_sit     !: SIT horizontal interpolation method (-1 = default) 
    … …  
    169173         & cn_velfbfiles,      & ! Velocity profile input filenames 
    170174         & cn_sssfbfiles,      & ! Sea surface salinity input filenames 
     175         & cn_ssvfbfiles,      & ! Sea surface velocity input filenames          
    171176         & cn_slchltotfbfiles, & ! Surface total              log10(chlorophyll) input filenames 
    172177         & cn_slchldiafbfiles, & ! Surface diatom             log10(chlorophyll) input filenames 
    … …  
    207212      LOGICAL :: ln_sit          ! Logical switch for sea ice thickness 
    208213      LOGICAL :: ln_sss          ! Logical switch for sea surface salinity obs 
     214      LOGICAL :: ln_ssv          ! Logical switch for sea surface velocity obs       
    209215      LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs 
    210216      LOGICAL :: ln_slchltot     ! Logical switch for surface total              log10(chlorophyll) obs 
    … …  
    261267      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 
    262268         & zmask                 ! Model land/sea mask associated with variables 
     269      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
     270         & zmask_surf            ! Surface model land/sea mask associated with variables 
    263271 
    264272 
    265273      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              & 
    266          &            ln_sst, ln_sic, ln_sit, ln_sss, ln_vel3d,               & 
     274         &            ln_sst, ln_sic, ln_sit, ln_sss,                 & 
     275         &            ln_ssv, ln_vel3d,                               & 
    267276         &            ln_slchltot, ln_slchldia, ln_slchlnon,          & 
    268277         &            ln_slchldin, ln_slchlmic, ln_slchlnan,          & 
    … …  
    280289         &            ln_time_mean_sla_bkg, ln_default_fp_indegs,     & 
    281290         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             & 
    282          &            ln_sss_fp_indegs, ln_sic_fp_indegs,             & 
    283          &            ln_sit_fp_indegs,                               & 
     291         &            ln_sss_fp_indegs, ln_ssv_fp_indegs,             & 
     292         &            ln_sic_fp_indegs, ln_sit_fp_indegs,             & 
    284293         &            cn_profbfiles, cn_slafbfiles,                   & 
    285294         &            cn_sstfbfiles, cn_sicfbfiles,                   & 
    286          &            cn_sitfbfiles,                                  & 
     295         &            cn_sitfbfiles, cn_ssvfbfiles,                   & 
    287296         &            cn_velfbfiles, cn_sssfbfiles,                   & 
    288297         &            cn_slchltotfbfiles, cn_slchldiafbfiles,         & 
    … …  
    305314         &            rn_sst_avglamscl, rn_sst_avgphiscl,             & 
    306315         &            rn_sss_avglamscl, rn_sss_avgphiscl,             & 
     316         &            rn_ssv_avglamscl, rn_ssv_avgphiscl,             &          
    307317         &            rn_sic_avglamscl, rn_sic_avgphiscl,             & 
    308318         &            rn_sit_avglamscl, rn_sit_avgphiscl,             & 
    309319         &            nn_1dint, nn_2dint_default,                     & 
    310320         &            nn_2dint_sla, nn_2dint_sst,                     & 
    311          &            nn_2dint_sss, nn_2dint_sic, nn_2dint_sit,       & 
     321         &            nn_2dint_sss, nn_2dint_ssv,                     & 
     322         &            nn_2dint_sic, nn_2dint_sit,                     & 
    312323         &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             & 
    313324         &            nn_profdavtypes 
    … …  
    325336      cn_velfbfiles(:)      = '' 
    326337      cn_sssfbfiles(:)      = '' 
     338      cn_ssvfbfiles(:)      = ''       
    327339      cn_slchltotfbfiles(:) = '' 
    328340      cn_slchldiafbfiles(:) = '' 
    … …  
    389401         WRITE(numout,*) '             Logical switch for velocity observations               ln_vel3d = ', ln_vel3d 
    390402         WRITE(numout,*) '             Logical switch for SSS observations                      ln_sss = ', ln_sss 
     403         WRITE(numout,*) '             Logical switch for SSV observations                      ln_ssv = ', ln_ssv          
    391404         WRITE(numout,*) '             Logical switch for surface total logchl obs         ln_slchltot = ', ln_slchltot 
    392405         WRITE(numout,*) '             Logical switch for surface diatom logchl obs        ln_slchldia = ', ln_slchldia 
    … …  
    424437         WRITE(numout,*) '             Type of horizontal interpolation method for SST    nn_2dint_sst = ', nn_2dint_sst 
    425438         WRITE(numout,*) '             Type of horizontal interpolation method for SSS    nn_2dint_sss = ', nn_2dint_sss 
     439         WRITE(numout,*) '             Type of horizontal interpolation method for SSV    nn_2dint_ssv = ', nn_2dint_ssv          
    426440         WRITE(numout,*) '             Type of horizontal interpolation method for SIC    nn_2dint_sic = ', nn_2dint_sic 
    427441         WRITE(numout,*) '             Type of horizontal interpolation method for SIT    nn_2dint_sit = ', nn_2dint_sit 
    … …  
    462476         &                  ln_pchltot,  ln_pno3,     ln_psi4,     ln_ppo4,     & 
    463477         &                  ln_pdic,     ln_palk,     ln_pph,      ln_po2 /) ) 
    464       nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sit, ln_sss,             & 
     478      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sit, ln_sss, ln_ssv,     & 
    465479         &                  ln_slchltot, ln_slchldia, ln_slchlnon, ln_slchldin, & 
    466480         &                  ln_slchlmic, ln_slchlnan, ln_slchlpic, ln_schltot,  & 
    … …  
    591605            cobstypessurf(jtype) = 'sss' 
    592606            clsurffiles(jtype,:) = cn_sssfbfiles 
     607         ENDIF 
     608         IF (ln_ssv) THEN 
     609            jtype = jtype + 1 
     610            cobstypessurf(jtype) = 'ssv' 
     611            clsurffiles(jtype,:) = cn_ssvfbfiles 
    593612         ENDIF 
    594613         IF (ln_slchltot) THEN 
    … …  
    721740               ztype_avgphiscl = rn_sss_avgphiscl 
    722741               ltype_fp_indegs = ln_sss_fp_indegs 
     742               ltype_night     = .FALSE. 
     743            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'ssv' ) THEN 
     744               IF ( nn_2dint_ssv == -1 ) THEN 
     745                  n2dint_type  = nn_2dint_default 
     746               ELSE 
     747                  n2dint_type  = nn_2dint_ssv 
     748               ENDIF 
     749               ztype_avglamscl = rn_ssv_avglamscl 
     750               ztype_avgphiscl = rn_ssv_avgphiscl 
     751               ltype_fp_indegs = ln_ssv_fp_indegs 
    723752               ltype_night     = .FALSE. 
    724753            ELSE 
    … …  
    903932               nvarssurf(jtype) = 1 
    904933               nextrsurf(jtype) = 2 
     934            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'ssv' ) THEN 
     935               nvarssurf(jtype) = 2 
     936               nextrsurf(jtype) = 0 
    905937            ELSE 
    906938               nvarssurf(jtype) = 1 
    … …  
    909941             
    910942            ALLOCATE( clvars( nvarssurf(jtype) ) ) 
    911  
     943            CALL wrk_alloc( jpi, jpj, nvarssurf(jtype), zglam ) 
     944            CALL wrk_alloc( jpi, jpj, nvarssurf(jtype), zgphi ) 
     945            CALL wrk_alloc( jpi, jpj, nvarssurf(jtype), zmask_surf ) 
     946 
     947            IF ( TRIM(cobstypessurf(jtype)) == 'ssv' ) THEN 
     948               zglam(:,:,1) = glamu(:,:) 
     949               zglam(:,:,2) = glamv(:,:) 
     950               zgphi(:,:,1) = gphiu(:,:) 
     951               zgphi(:,:,2) = gphiv(:,:) 
     952               zmask_surf(:,:,1) = umask(:,:,1) 
     953               zmask_surf(:,:,2) = vmask(:,:,1) 
     954            ELSE             
     955               DO jvar = 1, nvarssurf(jtype) 
     956                  zglam(:,:,jvar) = glamt(:,:) 
     957                  zgphi(:,:,jvar) = gphit(:,:) 
     958                  zmask_surf(:,:,jvar) = tmask(:,:,1)                                        
     959               END DO 
     960            ENDIF 
     961             
    912962            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
    913963               clvars(1) = 'SLA' 
    … …  
    921971            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN 
    922972               clvars(1) = 'SSS' 
     973            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'ssv' ) THEN 
     974               clvars(1) = 'UVEL' 
     975               clvars(2) = 'VVEL'            
    923976            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchltot' ) THEN 
    924977               clvars(1) = 'SLCHLTOT' 
    … …  
    9601013               &               llnightav(jtype), ltype_clim, ln_time_mean_sla_bkg, clvars ) 
    9611014 
    962             CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject, ln_seaicetypes ) 
     1015            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), & 
     1016               &               jpi, jpj, &             
     1017               &               zmask_surf, zglam, zgphi, & 
     1018               &               ln_nea, ln_bound_reject, ln_seaicetypes ) 
    9631019 
    9641020            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
    … …  
    9901046             
    9911047            DEALLOCATE( clvars ) 
     1048            CALL wrk_dealloc( jpi, jpj, nvarssurf(jtype), zglam ) 
     1049            CALL wrk_dealloc( jpi, jpj, nvarssurf(jtype), zgphi ) 
     1050            CALL wrk_dealloc( jpi, jpj, nvarssurf(jtype), zmask_surf ) 
    9921051 
    9931052         END DO 
    … …  
    10941153      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 
    10951154         & zprofmask               ! Mask associated with zprofvar 
    1096       REAL(wp), POINTER, DIMENSION(:,:) :: & 
     1155      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
    10971156         & zsurfvar, &             ! Model values equivalent to surface ob. 
    10981157         & zsurfclim, &            ! Climatology values for variables in a surface ob. 
    10991158         & zsurfmask               ! Mask associated with surface variable 
    11001159      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
    1101          & zglam,    &             ! Model longitudes for prof variables 
    1102          & zgphi                   ! Model latitudes for prof variables 
     1160         & zglam,    &             ! Model longitudes 
     1161         & zgphi                   ! Model latitudes 
    11031162      LOGICAL :: llog10            ! Perform log10 transform of variable 
    11041163#if defined key_fabm 
    … …  
    13461405 
    13471406      IF ( nsurftypes > 0 ) THEN 
    1348  
    1349          !Allocate local work arrays 
    1350          CALL wrk_alloc( jpi, jpj, zsurfvar ) 
    1351          CALL wrk_alloc( jpi, jpj, zsurfclim )          
    1352          CALL wrk_alloc( jpi, jpj, zsurfmask ) 
     1407          
    13531408#if defined key_fabm 
    13541409         CALL wrk_alloc( jpi, jpj, jpk, fabm_3d ) 
    … …  
    13571412         DO jtype = 1, nsurftypes 
    13581413 
     1414            !Allocate local work arrays 
     1415            CALL wrk_alloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfvar  ) 
     1416            CALL wrk_alloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfclim )          
     1417            CALL wrk_alloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfmask ) 
     1418            CALL wrk_alloc( jpi, jpj, surfdataqc(jtype)%nvar, zglam     ) 
     1419            CALL wrk_alloc( jpi, jpj, surfdataqc(jtype)%nvar, zgphi     ) 
     1420 
    13591421            !Defaults which might be changed 
    1360             zsurfmask(:,:) = tmask(:,:,1) 
    1361             zsurfclim(:,:) = 0._wp           
     1422            DO jvar = 1, surfdataqc(jtype)%nvar             
     1423               zsurfmask(:,:,jvar) = tmask(:,:,1) 
     1424               zsurfclim(:,:,jvar) = 0._wp 
     1425               zglam(:,:,jvar) = glamt(:,:) 
     1426               zgphi(:,:,jvar) = gphit(:,:)    
     1427            END DO              
    13621428            llog10 = .FALSE. 
    13631429 
    13641430            SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 
    13651431            CASE('sst') 
    1366                zsurfvar(:,:) = tsn(:,:,1,jp_tem) 
    1367                IF ( ln_output_clim ) zsurfclim(:,:) = tclim(:,:,1) 
     1432               zsurfvar(:,:,1) = tsn(:,:,1,jp_tem) 
     1433               IF ( ln_output_clim ) zsurfclim(:,:,1) = tclim(:,:,1) 
    13681434            CASE('sla') 
    1369                zsurfvar(:,:) = sshn(:,:) 
     1435               zsurfvar(:,:,1) = sshn(:,:) 
    13701436            CASE('sss') 
    1371                zsurfvar(:,:) = tsn(:,:,1,jp_sal) 
    1372                IF ( ln_output_clim ) zsurfclim(:,:) = sclim(:,:,1)               
     1437               zsurfvar(:,:,1) = tsn(:,:,1,jp_sal) 
     1438               IF ( ln_output_clim ) zsurfclim(:,:,1) = sclim(:,:,1)               
     1439            CASE('ssv') 
     1440               zsurfvar(:,:,1) = un(:,:,1) 
     1441               zsurfvar(:,:,2) = vn(:,:,1) 
     1442               zsurfmask(:,:,1) = umask(:,:,1) 
     1443               zsurfmask(:,:,2) = vmask(:,:,1)    
     1444               zglam(:,:,1) = glamu(:,:) 
     1445               zglam(:,:,2) = glamv(:,:)                
     1446               zgphi(:,:,1) = gphiu(:,:) 
     1447               zgphi(:,:,2) = gphiv(:,:)                   
    13731448            CASE('sic') 
    13741449               IF ( kstp == 0 ) THEN 
    … …  
    13811456               ELSE 
    13821457#if defined key_cice 
    1383                   zsurfvar(:,:) = fr_i(:,:) 
     1458                  zsurfvar(:,:,1) = fr_i(:,:) 
    13841459#elif defined key_lim2 || defined key_lim3 
    1385                   zsurfvar(:,:) = 1._wp - frld(:,:) 
     1460                  zsurfvar(:,:,1) = 1._wp - frld(:,:) 
    13861461#else 
    13871462               CALL ctl_stop( ' Trying to run sea-ice concentration observation operator', & 
    … …  
    14001475               ELSE        
    14011476#if defined key_cice 
    1402                   zsurfvar(:,:) = thick_i(:,:) 
     1477                  zsurfvar(:,:,1) = thick_i(:,:) 
    14031478#elif defined key_lim2 || defined key_lim3 
    14041479                  CALL ctl_stop( ' No sea-ice thickness observation operator defined for LIM model' ) 
    … …  
    14121487#if defined key_hadocc 
    14131488               ! Surface chlorophyll from HadOCC 
    1414                zsurfvar(:,:) = HADOCC_CHL(:,:,1) 
     1489               zsurfvar(:,:,1) = HADOCC_CHL(:,:,1) 
    14151490#elif defined key_medusa 
    14161491               ! Add non-diatom and diatom surface chlorophyll from MEDUSA 
    1417                zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 
     1492               zsurfvar(:,:,1) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 
    14181493#elif defined key_fabm 
    14191494               ! Add all surface chlorophyll groups from ERSEM 
    1420                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 
     1495               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 
    14211496                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
    14221497#else 
    … …  
    14321507#elif defined key_medusa 
    14331508               ! Diatom surface chlorophyll from MEDUSA 
    1434                zsurfvar(:,:) = trn(:,:,1,jpchd) 
     1509               zsurfvar(:,:,1) = trn(:,:,1,jpchd) 
    14351510#elif defined key_fabm 
    14361511               ! Diatom surface chlorophyll from ERSEM 
    1437                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) 
     1512               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) 
    14381513#else 
    14391514               CALL ctl_stop( ' Trying to run slchldia observation operator', & 
    … …  
    14481523#elif defined key_medusa 
    14491524               ! Non-diatom surface chlorophyll from MEDUSA 
    1450                zsurfvar(:,:) = trn(:,:,1,jpchn) 
     1525               zsurfvar(:,:,1) = trn(:,:,1,jpchn) 
    14511526#elif defined key_fabm 
    14521527               ! Add all non-diatom surface chlorophyll groups from ERSEM 
    1453                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 
     1528               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 
    14541529                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
    14551530#else 
    … …  
    14681543#elif defined key_fabm 
    14691544               ! Dinoflagellate surface chlorophyll from ERSEM 
    1470                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
     1545               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
    14711546#else 
    14721547               CALL ctl_stop( ' Trying to run slchldin observation operator', & 
    … …  
    14841559#elif defined key_fabm 
    14851560               ! Add diatom and dinoflagellate surface chlorophyll from ERSEM 
    1486                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
     1561               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
    14871562#else 
    14881563               CALL ctl_stop( ' Trying to run slchlmic observation operator', & 
    … …  
    15001575#elif defined key_fabm 
    15011576               ! Nanophytoplankton surface chlorophyll from ERSEM 
    1502                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) 
     1577               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) 
    15031578#else 
    15041579               CALL ctl_stop( ' Trying to run slchlnan observation operator', & 
    … …  
    15161591#elif defined key_fabm 
    15171592               ! Picophytoplankton surface chlorophyll from ERSEM 
    1518                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) 
     1593               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) 
    15191594#else 
    15201595               CALL ctl_stop( ' Trying to run slchlpic observation operator', & 
    … …  
    15261601#if defined key_hadocc 
    15271602               ! Surface chlorophyll from HadOCC 
    1528                zsurfvar(:,:) = HADOCC_CHL(:,:,1) 
     1603               zsurfvar(:,:,1) = HADOCC_CHL(:,:,1) 
    15291604#elif defined key_medusa 
    15301605               ! Add non-diatom and diatom surface chlorophyll from MEDUSA 
    1531                zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 
     1606               zsurfvar(:,:,1) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 
    15321607#elif defined key_fabm 
    15331608               ! Add all surface chlorophyll groups from ERSEM 
    1534                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 
    1535                   &            trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
     1609               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 
     1610                  &              trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
    15361611#else 
    15371612               CALL ctl_stop( ' Trying to run schltot observation operator', & 
    … …  
    15421617#if defined key_hadocc 
    15431618               ! Surface phytoplankton nitrogen from HadOCC multiplied by C:N ratio 
    1544                zsurfvar(:,:) = trn(:,:,1,jp_had_phy) * c2n_p 
     1619               zsurfvar(:,:,1) = trn(:,:,1,jp_had_phy) * c2n_p 
    15451620#elif defined key_medusa 
    15461621               ! Add non-diatom and diatom surface phytoplankton nitrogen from MEDUSA 
    15471622               ! multiplied by C:N ratio for each 
    1548                zsurfvar(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd) 
     1623               zsurfvar(:,:,1) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd) 
    15491624#elif defined key_fabm 
    15501625               ! Add all surface phytoplankton carbon groups from ERSEM 
    1551                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 
     1626               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 
    15521627                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c) 
    15531628#else 
    … …  
    15631638#elif defined key_medusa 
    15641639               ! Diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio 
    1565                zsurfvar(:,:) = trn(:,:,1,jpphd) * xthetapd 
     1640               zsurfvar(:,:,1) = trn(:,:,1,jpphd) * xthetapd 
    15661641#elif defined key_fabm 
    15671642               ! Diatom surface phytoplankton carbon from ERSEM 
    1568                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) 
     1643               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) 
    15691644#else 
    15701645               CALL ctl_stop( ' Trying to run slphydia observation operator', & 
    … …  
    15791654#elif defined key_medusa 
    15801655               ! Non-diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio 
    1581                zsurfvar(:,:) = trn(:,:,1,jpphn) * xthetapn 
     1656               zsurfvar(:,:,1) = trn(:,:,1,jpphn) * xthetapn 
    15821657#elif defined key_fabm 
    15831658               ! Add all non-diatom surface phytoplankton carbon groups from ERSEM 
    1584                zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 
    1585                   &            trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c) 
     1659               zsurfvar(:,:,1) = trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 
     1660                  &              trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c) 
    15861661#else 
    15871662               CALL ctl_stop( ' Trying to run slphynon observation operator', & 
    … …  
    15921667            CASE('sspm') 
    15931668#if defined key_spm 
    1594                zsurfvar(:,:) = 0.0 
     1669               zsurfvar(:,:,1) = 0.0 
    15951670               DO jn = 1, jp_spm 
    1596                   zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn)   ! sum SPM sizes 
     1671                  zsurfvar(:,:,1) = zsurfvar(:,:,1) + trn(:,:,1,jn)   ! sum SPM sizes 
    15971672               END DO 
    15981673#else 
    … …  
    16111686               ! light_xEPS diagnostic variable 
    16121687               fabm_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_xeps) 
    1613                zsurfvar(:,:) = fabm_3d(:,:,1) 
     1688               zsurfvar(:,:,1) = fabm_3d(:,:,1) 
    16141689#else 
    16151690               CALL ctl_stop( ' Trying to run skd490 observation operator', & 
    … …  
    16191694            CASE('sfco2') 
    16201695#if defined key_hadocc 
    1621                zsurfvar(:,:) = HADOCC_FCO2(:,:)    ! fCO2 from HadOCC 
     1696               zsurfvar(:,:,1) = HADOCC_FCO2(:,:)    ! fCO2 from HadOCC 
    16221697               IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. & 
    16231698                  & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN 
    16241699                  zsurfvar(:,:) = obfillflt 
    1625                   zsurfmask(:,:) = 0 
     1700                  zsurfmask(:,:,1) = 0 
    16261701                  CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', & 
    16271702                     &           ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' ) 
    16281703               ENDIF 
    16291704#elif defined key_medusa && defined key_roam 
    1630                zsurfvar(:,:) = f2_fco2w(:,:) 
     1705               zsurfvar(:,:,1) = f2_fco2w(:,:) 
    16311706#elif defined key_fabm 
    16321707               ! First, get pCO2 from FABM 
    16331708               fabm_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc) 
    1634                zsurfvar(:,:) = fabm_3d(:,:,1) 
     1709               zsurfvar(:,:,1) = fabm_3d(:,:,1) 
    16351710               ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of: 
    16361711               ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems 
    … …  
    16461721               ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal 
    16471722               ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway. 
    1648                zsurfvar(:,:) = zsurfvar(:,:) * EXP((-1636.75                                                          + & 
    1649                   &            12.0408      * (tsn(:,:,1,jp_tem)+rt0)                                                 - & 
    1650                   &            0.0327957    * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)                         + & 
    1651                   &            0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & 
    1652                   &            2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0)))                                        / & 
    1653                   &            (82.0578 * (tsn(:,:,1,jp_tem)+rt0))) 
     1723               zsurfvar(:,:,1) = zsurfvar(:,:,1) * EXP((-1636.75                                                          + & 
     1724                  &              12.0408      * (tsn(:,:,1,jp_tem)+rt0)                                                 - & 
     1725                  &              0.0327957    * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)                         + & 
     1726                  &              0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & 
     1727                  &              2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0)))                                        / & 
     1728                  &              (82.0578 * (tsn(:,:,1,jp_tem)+rt0))) 
    16541729#else 
    16551730               CALL ctl_stop( ' Trying to run sfco2 observation operator', & 
    … …  
    16591734            CASE('spco2') 
    16601735#if defined key_hadocc 
    1661                zsurfvar(:,:) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC 
     1736               zsurfvar(:,:,1) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC 
    16621737               IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. & 
    16631738                  & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN 
    1664                   zsurfvar(:,:) = obfillflt 
    1665                   zsurfmask(:,:) = 0 
     1739                  zsurfvar(:,:,1) = obfillflt 
     1740                  zsurfmask(:,:,1) = 0 
    16661741                  CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', & 
    16671742                     &           ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' ) 
    16681743               ENDIF 
    16691744#elif defined key_medusa && defined key_roam 
    1670                zsurfvar(:,:) = f2_pco2w(:,:) 
     1745               zsurfvar(:,:,1) = f2_pco2w(:,:) 
    16711746#elif defined key_fabm 
    16721747               fabm_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc) 
    1673                zsurfvar(:,:) = fabm_3d(:,:,1) 
     1748               zsurfvar(:,:,1) = fabm_3d(:,:,1) 
    16741749#else 
    16751750               CALL ctl_stop( ' Trying to run spco2 observation operator', & 
    … …  
    16861761               ! Take the log10 where we can, otherwise exclude 
    16871762               tiny = 1.0e-20 
    1688                WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt ) 
    1689                   zsurfvar(:,:)  = LOG10(zsurfvar(:,:)) 
     1763               WHERE(zsurfvar(:,:,1) > tiny .AND. zsurfvar(:,:,1) /= obfillflt ) 
     1764                  zsurfvar(:,:,1)  = LOG10(zsurfvar(:,:,1)) 
    16901765               ELSEWHERE 
    1691                   zsurfvar(:,:)  = obfillflt 
    1692                   zsurfmask(:,:) = 0 
     1766                  zsurfvar(:,:,1)  = obfillflt 
     1767                  zsurfmask(:,:,1) = 0 
    16931768               END WHERE 
    16941769            ENDIF 
    16951770 
    1696             IF ( TRIM(cobstypessurf(jtype)) == 'sla' .AND.                 & 
    1697                   &  ln_time_mean_sla_bkg ) THEN 
    1698                !Number of time-steps in meaning period 
    1699                imeanstp = NINT( ( MeanPeriodHours * 60. * 60. ) / rdt ) 
    1700                CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
    1701                   &               nit000, idaystp, zsurfvar,               & 
    1702                   &               zsurfclim, zsurfmask,                    & 
    1703                   &               n2dintsurf(jtype), llnightav(jtype),     & 
    1704                   &               ravglamscl(jtype), ravgphiscl(jtype),    & 
    1705                   &               lfpindegs(jtype), kmeanstp = imeanstp ) 
    1706  
    1707  
    1708             ELSE 
    1709                CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
    1710                   &               nit000, idaystp, zsurfvar,               & 
    1711                   &               zsurfclim, zsurfmask,                    & 
    1712                   &               n2dintsurf(jtype), llnightav(jtype),     & 
    1713                   &               ravglamscl(jtype), ravgphiscl(jtype),    & 
    1714                   &               lfpindegs(jtype) ) 
    1715             ENDIF 
    1716  
     1771            DO jvar = 1, surfdataqc(jtype)%nvar 
     1772 
     1773               IF ( TRIM(cobstypessurf(jtype)) == 'sla' .AND.                 & 
     1774                     &  ln_time_mean_sla_bkg ) THEN 
     1775                  !Number of time-steps in meaning period 
     1776                  imeanstp = NINT( ( MeanPeriodHours * 60. * 60. ) / rdt ) 
     1777                  CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
     1778                     &               nit000, idaystp, jvar,                   & 
     1779                     &               zsurfvar(:,:,jvar),                      & 
     1780                     &               zsurfclim(:,:,jvar),                     & 
     1781                     &               zsurfmask(:,:,jvar),                     & 
     1782                     &               zglam(:,:,jvar), zgphi(:,:,jvar),        &                      
     1783                     &               n2dintsurf(jtype), llnightav(jtype),     & 
     1784                     &               ravglamscl(jtype), ravgphiscl(jtype),    & 
     1785                     &               lfpindegs(jtype), kmeanstp = imeanstp ) 
     1786 
     1787               ELSE 
     1788                  CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
     1789                     &               nit000, idaystp, jvar,                   & 
     1790                     &               zsurfvar(:,:,jvar),                      & 
     1791                     &               zsurfclim(:,:,jvar),                     & 
     1792                     &               zsurfmask(:,:,jvar),                     & 
     1793                     &               zglam(:,:,jvar), zgphi(:,:,jvar),        &                                           
     1794                     &               n2dintsurf(jtype), llnightav(jtype),     & 
     1795                     &               ravglamscl(jtype), ravgphiscl(jtype),    & 
     1796                     &               lfpindegs(jtype) ) 
     1797               ENDIF 
     1798 
     1799            END DO 
     1800             
    17171801            ! Change label of data from FBD ("freeboard") to SIT ("Sea Ice 
    17181802            ! Thickness") 
    17191803            IF ( TRIM(surfdataqc(jtype)%cvars(1)) == 'FBD' ) THEN 
    1720                  surfdata(jtype)%cvars(1) = 'SIT' 
     1804                 surfdataqc(jtype)%cvars(1) = 'SIT' 
    17211805            ENDIF 
    17221806 
     1807            CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfvar  ) 
     1808            CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfclim )                   
     1809            CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfmask ) 
     1810            CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zglam     ) 
     1811            CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zgphi     )          
     1812 
    17231813         END DO 
    1724  
    1725          CALL wrk_dealloc( jpi, jpj, zsurfvar ) 
    1726          CALL wrk_dealloc( jpi, jpj, zsurfmask ) 
    17271814#if defined key_fabm 
    17281815         CALL wrk_dealloc( jpi, jpj, jpk, fabm_3d ) 
    … …  
    17801867                  & ) 
    17811868 
    1782                CALL obs_rotvel( profdataqc(jtype), nn_2dint_default, zu, zv ) 
     1869               CALL obs_rotvel_pro( profdataqc(jtype), nn_2dint_default, zu, zv ) 
    17831870 
    17841871               DO jo = 1, profdataqc(jtype)%nprof 
    … …  
    18131900 
    18141901         DO jtype = 1, nsurftypes 
     1902 
     1903            IF ( TRIM(cobstypessurf(jtype)) == 'vel' ) THEN 
     1904 
     1905               ! For velocity data, rotate the model velocities to N/S, E/W 
     1906               ! using the compressed data structure. 
     1907               ALLOCATE( & 
     1908                  & zu(surfdataqc(jtype)%nsurf), & 
     1909                  & zv(surfdataqc(jtype)%nsurf)  & 
     1910                  & ) 
     1911 
     1912               CALL obs_rotvel_surf( surfdataqc(jtype), nn_2dint_default, zu, zv ) 
     1913 
     1914               DO jo = 1, surfdataqc(jtype)%nsurf 
     1915                  surfdataqc(jtype)%rmod(jo,1) = zu(jo) 
     1916                  surfdataqc(jtype)%rmod(jo,2) = zv(jo) 
     1917               END DO 
     1918 
     1919               DEALLOCATE( zu ) 
     1920               DEALLOCATE( zv ) 
     1921 
     1922            END IF 
     1923 
    18151924 
    18161925            CALL obs_surf_decompress( surfdataqc(jtype), & 

  • branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

    r12610 r13393  
    544544 
    545545   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,            & 
    546       &                     kit000, kdaystp, psurf, pclim, psurfmask,   & 
     546      &                     kit000, kdaystp, kvar, & 
     547      &                     psurf, pclim, psurfmask, & 
     548      &                     plam, pphi, & 
    547549      &                     k2dint, ldnightav, plamscl, pphiscl, & 
    548550      &                     lindegrees, kmeanstp ) 
    … …  
    579581      !!      ! 15-02 (M. Martin) Combined routine for surface types 
    580582      !!      ! 17-03 (M. Martin) Added horizontal averaging options 
     583      !!      ! 20-08 (M. Martin) Added surface velocity options       
    581584      !!----------------------------------------------------------------------- 
    582585 
    … …  
    595598                                       !   (kit000-1 = restart time) 
    596599      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day 
     600      INTEGER, INTENT(IN) :: kvar      ! Number of variables in surfdataqc      
    597601      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header) 
    598602      INTEGER, INTENT(IN), OPTIONAL :: & 
    … …  
    603607         & pclim,  &                   ! Climatological surface field          
    604608         & psurfmask                   ! Land-sea mask 
     609      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
     610         & plam,   &                   ! Model longitudes for variable 
     611         & pphi                        ! Model latitudes for variable          
    605612      LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data 
    606613      REAL(KIND=wp), INTENT(IN) :: & 
    … …  
    670677      CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp ) 
    671678 
    672  
    673679      IF ( l_timemean ) THEN 
    674680         ! Initialize time mean for first timestep 
    … …  
    681687            DO jj = 1, jpj 
    682688               DO ji = 1, jpi 
    683                   surfdataqc%vdmean(ji,jj) = 0.0 
     689                  surfdataqc%vdmean(ji,jj,kvar) = 0.0 
    684690               END DO 
    685691            END DO 
    … …  
    690696         DO jj = 1, jpj 
    691697            DO ji = 1, jpi 
    692                surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 
    693                   &                        + psurf(ji,jj) 
     698               surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) & 
     699                  &                            + psurf(ji,jj) 
    694700            END DO 
    695701         END DO 
    … …  
    701707            DO jj = 1, jpj 
    702708               DO ji = 1, jpi 
    703                   surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 
    704                      &                       * zmeanstp 
     709                  surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) & 
     710                     &                            * zmeanstp 
    705711               END DO 
    706712            END DO 
    … …  
    729735            DO jj = 1, jpj 
    730736               DO ji = 1, jpi 
    731                   surfdataqc%vdmean(ji,jj) = 0.0 
     737                  surfdataqc%vdmean(ji,jj,kvar) = 0.0 
    732738                  zmeanday(ji,jj) = 0.0 
    733739                  icount_night(ji,jj) = 0 
    … …  
    743749            DO ji = 1, jpi 
    744750               ! Increment the temperature field for computing night mean and counter 
    745                surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj)  & 
    746                       &                    + psurf(ji,jj) * REAL( imask_night(ji,jj) ) 
     751               surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar)  & 
     752                      &                        + psurf(ji,jj) * REAL( imask_night(ji,jj) ) 
    747753               zmeanday(ji,jj)          = zmeanday(ji,jj) + psurf(ji,jj) 
    748754               icount_night(ji,jj)      = icount_night(ji,jj) + imask_night(ji,jj) 
    … …  
    758764                  ! Test if "no night" point 
    759765                  IF ( icount_night(ji,jj) > 0 ) THEN 
    760                      surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 
    761                        &                        / REAL( icount_night(ji,jj) ) 
     766                     surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) & 
     767                       &                            / REAL( icount_night(ji,jj) ) 
    762768                  ELSE 
    763769                     !At locations where there is no night (e.g. poles), 
    764770                     ! calculate daily mean instead of night-time mean. 
    765                      surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp 
     771                     surfdataqc%vdmean(ji,jj,kvar) = zmeanday(ji,jj) * zdaystp 
    766772                  ENDIF 
    767773               END DO 
    … …  
    772778 
    773779      ! Get the data for interpolation 
    774  
    775780      ALLOCATE( & 
    776          & zweig(imaxifp,imaxjfp,1),      & 
    777          & igrdi(imaxifp,imaxjfp,isurf), & 
    778          & igrdj(imaxifp,imaxjfp,isurf), & 
    779          & zglam(imaxifp,imaxjfp,isurf), & 
    780          & zgphi(imaxifp,imaxjfp,isurf), & 
    781          & zmask(imaxifp,imaxjfp,isurf), & 
    782          & zsurf(imaxifp,imaxjfp,isurf), & 
    783          & zsurftmp(imaxifp,imaxjfp,isurf),  & 
    784          & zglamf(imaxifp+1,imaxjfp+1,isurf), & 
    785          & zgphif(imaxifp+1,imaxjfp+1,isurf), & 
    786          & igrdip1(imaxifp+1,imaxjfp+1,isurf), & 
    787          & igrdjp1(imaxifp+1,imaxjfp+1,isurf) & 
     781         & zweig(imaxifp,imaxjfp,1),       & 
     782         & igrdi(imaxifp,imaxjfp,isurf),   & 
     783         & igrdj(imaxifp,imaxjfp,isurf),   & 
     784         & zglam(imaxifp,imaxjfp,isurf),   & 
     785         & zgphi(imaxifp,imaxjfp,isurf),   & 
     786         & zmask(imaxifp,imaxjfp,isurf),   & 
     787         & zsurf(imaxifp,imaxjfp,isurf),   & 
     788         & zsurftmp(imaxifp,imaxjfp,isurf) & 
    788789         & ) 
     790 
     791      IF ( k2dint > 4 ) THEN 
     792         ALLOCATE( & 
     793            & zglamf(imaxifp+1,imaxjfp+1,isurf),  & 
     794            & zgphif(imaxifp+1,imaxjfp+1,isurf),  & 
     795            & igrdip1(imaxifp+1,imaxjfp+1,isurf), & 
     796            & igrdjp1(imaxifp+1,imaxjfp+1,isurf)  & 
     797            & ) 
     798      ENDIF 
    789799 
    790800      IF ( surfdataqc%lclim ) ALLOCATE( zclim(imaxifp,imaxjfp,isurf) ) 
    … …  
    793803         iobs = jobs - surfdataqc%nsurfup 
    794804         DO ji = 0, imaxifp 
    795             imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1 
     805            imodi = surfdataqc%mi(jobs,kvar) - int(imaxifp/2) + ji - 1 
    796806             
    797807            !Deal with wrap around in longitude 
    … …  
    800810             
    801811            DO jj = 0, imaxjfp 
    802                imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1 
     812               imodj = surfdataqc%mj(jobs,kvar) - int(imaxjfp/2) + jj - 1 
    803813               !If model values are out of the domain to the north/south then 
    804814               !set them to be the edge of the domain 
    … …  
    806816               IF ( imodj > jpjglo ) imodj = jpjglo 
    807817 
    808                igrdip1(ji+1,jj+1,iobs) = imodi 
    809                igrdjp1(ji+1,jj+1,iobs) = imodj 
     818               IF ( k2dint > 4 ) THEN 
     819                  igrdip1(ji+1,jj+1,iobs) = imodi 
     820                  igrdjp1(ji+1,jj+1,iobs) = imodj 
     821               ENDIF 
    810822                
    811823               IF ( ji >= 1 .AND. jj >= 1 ) THEN 
    … …  
    819831 
    820832      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    821          &                  igrdi, igrdj, glamt, zglam ) 
     833         &                  igrdi, igrdj, plam, zglam ) 
    822834      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    823          &                  igrdi, igrdj, gphit, zgphi ) 
     835         &                  igrdi, igrdj, pphi, zgphi ) 
    824836      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    825837         &                  igrdi, igrdj, psurfmask, zmask ) 
    … …  
    831843            IF (lwp) WRITE(numout,*)' Interpolating the time mean values on time step: ',kt 
    832844            CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    833                &                  igrdi, igrdj, surfdataqc%vdmean(:,:), zsurfm ) 
     845               &                  igrdi, igrdj, surfdataqc%vdmean(:,:,kvar), zsurfm ) 
    834846         ENDIF 
    835847      ELSE 
    … …  
    858870 
    859871         CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, & 
    860             &               surfdataqc%vdmean(:,:), zsurfm ) 
     872            &                  surfdataqc%vdmean(:,:,kvar), zsurfm ) 
    861873 
    862874      ENDIF 
    … …  
    937949               surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2) 
    938950            ELSE 
    939                surfdataqc%rmod(jobs,1) = zext(1) 
     951               surfdataqc%rmod(jobs,kvar) = zext(1) 
    940952            ENDIF 
    941953 
    … …  
    985997         & zmask, & 
    986998         & zsurf, & 
    987          & zsurftmp, & 
    988          & zglamf, & 
    989          & zgphif, & 
    990          & igrdip1,& 
    991          & igrdjp1 & 
     999         & zsurftmp & 
    9921000         & ) 
    9931001 
     1002      IF ( k2dint > 4 ) THEN 
     1003         DEALLOCATE( &      
     1004            & zglamf, & 
     1005            & zgphif, & 
     1006            & igrdip1,& 
     1007            & igrdjp1 & 
     1008            & ) 
     1009      ENDIF 
     1010             
    9941011      IF ( surfdataqc%lclim ) DEALLOCATE( zclim ) 
    9951012 
    … …  
    10011018      ENDIF 
    10021019 
    1003       surfdataqc%nsurfup = surfdataqc%nsurfup + isurf 
     1020      IF ( kvar == surfdataqc%nvar ) THEN 
     1021         surfdataqc%nsurfup = surfdataqc%nsurfup + isurf 
     1022      ENDIF 
    10041023 
    10051024   END SUBROUTINE obs_surf_opt 

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

    r12610 r13393  
    5353CONTAINS 
    5454 
    55    SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject, & 
    56                             ld_seaicetypes, kqc_cutoff ) 
    57       !!---------------------------------------------------------------------- 
    58       !!                    ***  ROUTINE obs_pre_sla  *** 
     55   SUBROUTINE obs_pre_surf( surfdata, surfdataqc, & 
     56      &                     kpi, kpj, &    
     57      &                     zmask, pglam, pgphi, & 
     58      &                     ld_nea, ld_bound_reject, & 
     59      &                     ld_seaicetypes, kqc_cutoff ) 
     60      !!---------------------------------------------------------------------- 
     61      !!                    ***  ROUTINE obs_pre_surf  *** 
    5962      !! 
    6063      !! ** Purpose : First level check and screening of surface observations 
    … …  
    8285      TYPE(obs_surf), INTENT(INOUT) :: surfdata    ! Full set of surface data 
    8386      TYPE(obs_surf), INTENT(INOUT) :: surfdataqc  ! Subset of surface data not failing screening 
     87      INTEGER, INTENT(IN) :: kpi, kpj              ! Local domain sizes       
     88      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,surfdata%nvar) :: & 
     89         & zmask       
     90      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,surfdata%nvar) :: & 
     91         & pglam, & 
     92         & pgphi 
    8493      LOGICAL, INTENT(IN) :: ld_nea                ! Switch for rejecting observation near land 
    8594      LOGICAL, INTENT(IN) :: ld_bound_reject       ! Switch for rejecting obs near the boundary 
    … …  
    94103      INTEGER :: imin0 
    95104      INTEGER :: icycle       ! Current assimilation cycle 
    96                               ! Counters for observations that 
    97       INTEGER :: iotdobs      !  - outside time domain 
    98       INTEGER :: iosdsobs     !  - outside space domain 
    99       INTEGER :: ilansobs     !  - within a model land cell 
    100       INTEGER :: inlasobs     !  - close to land 
    101       INTEGER :: igrdobs      !  - fail the grid search 
    102       INTEGER :: ibdysobs     !  - close to open boundary 
    103                               ! Global counters for observations that 
    104       INTEGER :: iotdobsmpp     !  - outside time domain 
    105       INTEGER :: iosdsobsmpp    !  - outside space domain 
    106       INTEGER :: ilansobsmpp    !  - within a model land cell 
    107       INTEGER :: inlasobsmpp    !  - close to land 
    108       INTEGER :: igrdobsmpp     !  - fail the grid search 
    109       INTEGER :: ibdysobsmpp  !  - close to open boundary 
     105                                                        ! Counters for observations that are 
     106      INTEGER                           :: iotdobs      !  - outside time domain 
     107      INTEGER, DIMENSION(surfdata%nvar) :: iosdsobs     !  - outside space domain 
     108      INTEGER, DIMENSION(surfdata%nvar) :: ilansobs     !  - within a model land cell 
     109      INTEGER, DIMENSION(surfdata%nvar) :: inlasobs     !  - close to land 
     110      INTEGER, DIMENSION(surfdata%nvar) :: ibdysobs     !  - close to open boundary 
     111      INTEGER                           :: igrdobs      !  - fail the grid search       
     112                                                        ! Global counters for observations that 
     113      INTEGER                           :: iotdobsmpp   !  - outside time domain 
     114      INTEGER, DIMENSION(surfdata%nvar) :: iosdsobsmpp  !  - outside space domain 
     115      INTEGER, DIMENSION(surfdata%nvar) :: ilansobsmpp  !  - within a model land cell 
     116      INTEGER, DIMENSION(surfdata%nvar) :: inlasobsmpp  !  - close to land 
     117      INTEGER, DIMENSION(surfdata%nvar) :: ibdysobsmpp  !  - close to open boundary 
     118      INTEGER                           :: igrdobsmpp   !  - fail the grid search 
     119 
    110120      LOGICAL, DIMENSION(:), ALLOCATABLE :: &  
    111121         & llvalid            ! SLA data selection 
    112       INTEGER :: jobs         ! Obs. loop variable 
     122      INTEGER :: jobs         ! Obs. loop counter 
     123      INTEGER :: jvar         ! Variable loop counter 
    113124      INTEGER :: jstp         ! Time loop variable 
    114125      INTEGER :: inrc         ! Time index variable 
    115  
     126      CHARACTER(LEN=256) :: cout1  ! Diagnostic output line 
     127       
    116128      IF(lwp) WRITE(numout,*)'obs_pre_surf : Preparing the surface observations...' 
    117129      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    … …  
    130142      iotdobs  = 0 
    131143      igrdobs  = 0 
    132       iosdsobs = 0 
    133       ilansobs = 0 
    134       inlasobs = 0 
    135       ibdysobs = 0  
     144      iosdsobs(:) = 0 
     145      ilansobs(:) = 0 
     146      inlasobs(:) = 0 
     147      ibdysobs(:) = 0  
    136148 
    137149      ! Set QC cutoff to optional value if provided 
    … …  
    162174      ! Check for surface data failing the grid search 
    163175      ! ----------------------------------------------------------------------- 
    164  
    165       CALL obs_coo_grd( surfdata%nsurf,   surfdata%mi, surfdata%mj, & 
    166          &              surfdata%nqc,     igrdobs                         ) 
    167  
     176      DO jvar = 1, surfdata%nvar 
     177         CALL obs_coo_grd( surfdata%nsurf,   surfdata%mi(:,jvar), surfdata%mj(:,jvar), & 
     178            &              surfdata%nqc,     igrdobs ) 
     179      END DO 
     180       
    168181      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) 
    169182 
    … …  
    172185      ! ----------------------------------------------------------------------- 
    173186 
    174       CALL obs_coo_spc_2d( surfdata%nsurf,              & 
    175          &                 jpi,          jpj,          & 
    176          &                 surfdata%mi,   surfdata%mj,   &  
    177          &                 surfdata%rlam, surfdata%rphi, & 
    178          &                 glamt,        gphit,        & 
    179          &                 tmask(:,:,1), surfdata%nqc,  & 
    180          &                 iosdsobs,     ilansobs,     & 
    181          &                 inlasobs,     ld_nea,       & 
    182          &                 ibdysobs,     ld_bound_reject, & 
    183          &                 iqc_cutoff                     ) 
    184  
    185       CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 
    186       CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 
    187       CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 
    188       CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp ) 
    189  
     187      DO jvar = 1, surfdata%nvar 
     188         CALL obs_coo_spc_2d( surfdata%nsurf,                         & 
     189            &                 jpi,          jpj,                      & 
     190            &                 surfdata%mi(:,jvar), surfdata%mj(:,jvar), &  
     191            &                 surfdata%rlam, surfdata%rphi,           & 
     192            &                 pglam(:,:,jvar), pgphi(:,:,jvar),       & 
     193            &                 zmask(:,:,jvar), surfdata%nqc(:),       & 
     194            &                 iosdsobs(jvar),     ilansobs(jvar),     & 
     195            &                 inlasobs(jvar),     ld_nea,             & 
     196            &                 ibdysobs(jvar),     ld_bound_reject,    & 
     197            &                 iqc_cutoff                     ) 
     198         CALL obs_mpp_sum_integer( iosdsobs(jvar), iosdsobsmpp(jvar) ) 
     199         CALL obs_mpp_sum_integer( ilansobs(jvar), ilansobsmpp(jvar) ) 
     200         CALL obs_mpp_sum_integer( inlasobs(jvar), inlasobsmpp(jvar) ) 
     201         CALL obs_mpp_sum_integer( ibdysobs(jvar), ibdysobsmpp(jvar) ) 
     202      END DO 
     203            
    190204      ! ----------------------------------------------------------------------- 
    191205      ! Copy useful data from the surfdata data structure to 
    … …  
    216230       
    217231      IF(lwp) THEN 
     232 
     233         DO jvar = 1, surfdataqc%nvar        
     234            IF ( jvar == 1 ) THEN 
     235               cout1=TRIM(surfdataqc%cvars(1))                   
     236            ELSE 
     237               WRITE(cout1,'(A,A1,A)') TRIM(cout1), '/', TRIM(surfdataqc%cvars(jvar))             
     238            ENDIF 
     239         END DO 
     240                
    218241         WRITE(numout,*) 
    219          WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data outside time domain                  = ', & 
     242         WRITE(numout,*) ' '//TRIM(cout1)//' data outside time domain                  = ', & 
    220243            &            iotdobsmpp 
    221          WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data that failed grid search    = ', & 
     244         WRITE(numout,*) ' Remaining '//TRIM(cout1)//' data that failed grid search    = ', & 
    222245            &            igrdobsmpp 
    223          WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data outside space domain       = ', & 
    224             &            iosdsobsmpp 
    225          WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data at land points             = ', & 
    226             &            ilansobsmpp 
    227          IF (ld_nea) THEN 
    228             WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near land points (removed) = ', & 
    229                &            inlasobsmpp 
    230          ELSE 
    231             WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near land points (kept)    = ', & 
    232                &            inlasobsmpp 
    233          ENDIF 
    234          WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near open boundary (removed) = ', & 
    235             &            ibdysobsmpp   
    236          WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data accepted                             = ', & 
    237             &            surfdataqc%nsurfmpp 
     246 
     247         DO jvar = 1, surfdataqc%nvar             
     248            WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data outside space domain       = ', & 
     249                &            iosdsobsmpp(jvar) 
     250             WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data at land points             = ', & 
     251                &            ilansobsmpp(jvar) 
     252             IF (ld_nea) THEN 
     253                WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data near land points (removed) = ', & 
     254                   &            inlasobsmpp(jvar) 
     255             ELSE 
     256                WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data near land points (kept)    = ', & 
     257                   &            inlasobsmpp(jvar) 
     258             ENDIF      
     259             WRITE(numout,*) ' Remaining '//surfdataqc%cvars(jvar)//' data near open boundary (removed) = ', & 
     260                &            ibdysobsmpp(jvar) 
     261          END DO 
     262          WRITE(numout,*) ' '//TRIM(cout1)//' data accepted                             = ', & 
     263             &            surfdataqc%nsurfmpp 
    238264 
    239265         WRITE(numout,*) 
    240266         WRITE(numout,*) ' Number of observations per time step :' 
    241267         WRITE(numout,*) 
    242          WRITE(numout,'(10X,A,10X,A)')'Time step',surfdataqc%cvars(1) 
     268         WRITE(numout,'(10X,A,10X,A)')'Time step',TRIM(cout1) 
    243269         WRITE(numout,'(10X,A,5X,A)')'---------','-----------------' 
    244270         CALL FLUSH(numout) 
    … …  
    445471 
    446472      IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 
    447          CALL obs_uv_rej( profdata, iuvchku, iuvchkv, iqc_cutoff ) 
     473         CALL obs_uv_rej_pro( profdata, iuvchku, iuvchkv, iqc_cutoff ) 
    448474         CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp ) 
    449475         CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp ) 
    … …  
    14571483   END SUBROUTINE obs_pro_rej 
    14581484 
    1459    SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff ) 
    1460       !!---------------------------------------------------------------------- 
    1461       !!                    ***  ROUTINE obs_uv_rej *** 
     1485   SUBROUTINE obs_uv_rej_pro( profdata, knumu, knumv, kqc_cutoff ) 
     1486      !!---------------------------------------------------------------------- 
     1487      !!                    ***  ROUTINE obs_uv_rej_pro *** 
    14621488      !! 
    14631489      !! ** Purpose : Reject u if v is rejected and vice versa 
    … …  
    15131539      END DO 
    15141540 
    1515    END SUBROUTINE obs_uv_rej 
     1541   END SUBROUTINE obs_uv_rej_pro 
    15161542 
    15171543END MODULE obs_prep 

  • branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90

    r12610 r13393  
    8585      INTEGER :: jj 
    8686      INTEGER :: jk 
     87      INTEGER :: jvar 
    8788      INTEGER :: iflag 
    8889      INTEGER :: inobf 
    … …  
    107108         & ityp, & 
    108109         & itypmpp 
    109       INTEGER, DIMENSION(:), ALLOCATABLE :: & 
     110      INTEGER, DIMENSION(:,:), ALLOCATABLE :: & 
    110111         & iobsi,    & 
    111112         & iobsj,    & 
    112          & iproc,    & 
     113         & iproc 
     114      INTEGER, DIMENSION(:), ALLOCATABLE :: &          
    113115         & iindx,    & 
    114116         & ifileidx, & 
    … …  
    126128      TYPE(obfbdata), POINTER, DIMENSION(:) :: & 
    127129         & inpfiles 
    128  
     130      CHARACTER(LEN=256) :: cout1  ! Diagnostic output line 
     131       
    129132      ! Local initialization 
    130133      iobs  = 0 
    … …  
    278281 
    279282            IF ( inpfiles(jj)%nobs > 0 ) THEN 
    280                inpfiles(jj)%iproc = -1 
    281                inpfiles(jj)%iobsi = -1 
    282                inpfiles(jj)%iobsj = -1 
     283               inpfiles(jj)%iproc(:,:) = -1 
     284               inpfiles(jj)%iobsi(:,:) = -1 
     285               inpfiles(jj)%iobsj(:,:) = -1 
    283286            ENDIF 
    284287 
    … …  
    295298               END DO 
    296299            ENDIF    
    297  
     300             
    298301            inowin = 0 
    299302            DO ji = 1, inpfiles(jj)%nobs 
    … …  
    305308            ALLOCATE( zlam(inowin)  ) 
    306309            ALLOCATE( zphi(inowin)  ) 
    307             ALLOCATE( iobsi(inowin) ) 
    308             ALLOCATE( iobsj(inowin) ) 
    309             ALLOCATE( iproc(inowin) ) 
     310            ALLOCATE( iobsi(inowin,kvars) ) 
     311            ALLOCATE( iobsj(inowin,kvars) ) 
     312            ALLOCATE( iproc(inowin,kvars) ) 
    310313            inowin = 0 
    311314            DO ji = 1, inpfiles(jj)%nobs 
    … …  
    318321            END DO 
    319322 
    320             CALL obs_grid_search( inowin, zlam, zphi, iobsi, iobsj, iproc, 'T' ) 
     323            ! Assume anything other than velocity is on T grid 
     324            IF ( TRIM( inpfiles(jj)%cname(1) ) == 'UVEL' ) THEN 
     325               CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,1), iobsj(:,1), & 
     326                  &                  iproc(:,1), 'U' ) 
     327               CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,2), iobsj(:,2), & 
     328                  &                  iproc(:,2), 'V' ) 
     329            ELSE 
     330               CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,1), iobsj(:,1), & 
     331                  &                  iproc(:,1), 'T' ) 
     332               IF ( kvars > 1 ) THEN 
     333                  DO jvar = 2, kvars 
     334                     iobsi(:,jvar) = iobsi(:,1) 
     335                     iobsj(:,jvar) = iobsj(:,1) 
     336                     iproc(:,jvar) = iproc(:,1) 
     337                  END DO 
     338               ENDIF 
     339            ENDIF 
    321340 
    322341            inowin = 0 
    … …  
    325344                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
    326345                  inowin = inowin + 1 
    327                   inpfiles(jj)%iproc(ji,1) = iproc(inowin) 
    328                   inpfiles(jj)%iobsi(ji,1) = iobsi(inowin) 
    329                   inpfiles(jj)%iobsj(ji,1) = iobsj(inowin) 
     346                  DO jvar = 1, kvars 
     347                     inpfiles(jj)%iproc(ji,jvar) = iproc(inowin,jvar) 
     348                     inpfiles(jj)%iobsi(ji,jvar) = iobsi(inowin,jvar) 
     349                     inpfiles(jj)%iobsj(ji,jvar) = iobsj(inowin,jvar) 
     350                  END DO 
    330351               ENDIF 
    331352            END DO 
    … …  
    448469 
    449470               ! Coordinate search parameters 
    450                surfdata%mi  (iobs) = inpfiles(jj)%iobsi(ji,1) 
    451                surfdata%mj  (iobs) = inpfiles(jj)%iobsj(ji,1) 
    452  
     471               DO jvar = 1, kvars 
     472                  surfdata%mi(iobs,jvar) = inpfiles(jj)%iobsi(ji,jvar) 
     473                  surfdata%mj(iobs,jvar) = inpfiles(jj)%iobsj(ji,jvar) 
     474               END DO 
     475                
    453476               ! WMO number 
    454477               surfdata%cwmo(iobs) = inpfiles(jj)%cdwmo(ji) 
    … …  
    480503 
    481504               ! Observed value 
    482                surfdata%robs(iobs,1) = inpfiles(jj)%pob(1,ji,1) 
     505               DO jvar = 1, kvars                
     506                  surfdata%robs(iobs,jvar) = inpfiles(jj)%pob(1,ji,jvar) 
     507               END DO 
    483508               IF ( TRIM(surfdata%cvars(1)) == 'FBD' ) THEN 
    484509                   surfdata%rext(iobs,1) = inpfiles(jj)%pob(1,ji,1) 
    … …  
    488513               ! Model and MDT is set to fbrmdi unless read from file 
    489514               IF ( ldmod ) THEN 
    490                   surfdata%rmod(iobs,1) = inpfiles(jj)%padd(1,ji,1,1) 
     515                  DO jvar = 1, kvars                               
     516                     surfdata%rmod(iobs,jvar) = inpfiles(jj)%padd(1,ji,1,jvar) 
     517                  END DO 
    491518                  IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) THEN 
    492519                     surfdata%rext(iobs,1) = inpfiles(jj)%padd(1,ji,2,1) 
    … …  
    494521                  ENDIF 
    495522                ELSE 
    496                   surfdata%rmod(iobs,1) = fbrmdi 
     523                  DO jvar = 1, kvars                 
     524                     surfdata%rmod(iobs,jvar) = fbrmdi 
     525                  END DO 
    497526                  IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) surfdata%rext(iobs,:) = fbrmdi 
    498527               ENDIF 
    499528                
    500529               ! Initialise climatology if set 
    501                IF ( surfdata%lclim ) surfdata%rclm(iobs,1) = fbrmdi 
     530               IF ( surfdata%lclim ) THEN 
     531                  DO jvar = 1, kvars 
     532                     surfdata%rclm(iobs,jvar) = fbrmdi 
     533                  END DO 
     534               ENDIF 
    502535                
    503536               ! STD (obs error standard deviation) read from file and passed through obs operator 
    … …  
    521554      !----------------------------------------------------------------------- 
    522555      IF (lwp) THEN 
    523  
     556         DO jvar = 1, surfdata%nvar        
     557            IF ( jvar == 1 ) THEN 
     558               cout1=TRIM(surfdata%cvars(1))                   
     559            ELSE 
     560               WRITE(cout1,'(A,A1,A)') TRIM(cout1), '/', TRIM(surfdata%cvars(jvar))             
     561            ENDIF 
     562         END DO 
     563  
    524564         WRITE(numout,*) 
    525          WRITE(numout,'(1X,A)')TRIM( surfdata%cvars(1) )//' data' 
     565         WRITE(numout,'(1X,A)')TRIM( cout1 )//' data' 
    526566         WRITE(numout,'(1X,A)')'--------------' 
    527567         DO jj = 1,8 
    … …  
    533573            & '---------------------------------------------------------------' 
    534574         WRITE(numout,'(1X,A,I8)') & 
    535             & 'Total data for variable '//TRIM( surfdata%cvars(1) )// & 
     575            & 'Total data for variable '//TRIM( cout1 )// & 
    536576            & '           = ', iobsmpp 
    537577         WRITE(numout,'(1X,A)') & 

  • branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_rot_vel.F90

    r7992 r13393  
    66 
    77   !!---------------------------------------------------------------------- 
    8    !!   obs_rotvel : Rotate velocity data into N-S,E-W directorions 
     8   !!   obs_rotvel_pro :  Rotate profile velocity data into N-S,E-W directions 
     9   !!   obs_rotvel_surf : Rotate surface velocity data into N-S,E-W directions    
    910   !!---------------------------------------------------------------------- 
    1011   !! * Modules used    
    … …  
    1718   USE obs_utils                ! For error handling 
    1819   USE obs_profiles_def         ! Profile definitions 
     20   USE obs_surf_def             ! Surface definitions    
    1921   USE obs_inter_h2d            ! Horizontal interpolation 
    2022   USE obs_inter_sup            ! MPP support routines for interpolation 
    … …  
    2729   PRIVATE 
    2830 
    29    PUBLIC obs_rotvel            ! Rotate the observations 
    30  
     31   PUBLIC obs_rotvel_pro, &     ! Rotate the profile velocity observations 
     32      &   obs_rotvel_surf       ! Rotate the surface velocity observations 
     33       
    3134   !!---------------------------------------------------------------------- 
    3235   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    … …  
    3740CONTAINS 
    3841 
    39    SUBROUTINE obs_rotvel( profdata, k2dint, pu, pv ) 
     42   SUBROUTINE obs_rotvel_pro( profdata, k2dint, pu, pv ) 
    4043      !!--------------------------------------------------------------------- 
    4144      !! 
    … …  
    228231      CALL wrk_dealloc(jpi,jpj,zsingu,zcosgu,zsingv,zcosgv)  
    229232 
    230    END SUBROUTINE obs_rotvel 
     233   END SUBROUTINE obs_rotvel_pro 
     234 
     235   SUBROUTINE obs_rotvel_surf( surfdata, k2dint, pu, pv ) 
     236      !!--------------------------------------------------------------------- 
     237      !! 
     238      !!                   *** ROUTINE obs_rotvel_surf *** 
     239      !! 
     240      !! ** Purpose : Rotate surface velocity data into N-S,E-W directorions 
     241      !! 
     242      !! ** Method  : Interpolation of geo2ocean coefficients on U,V grid 
     243      !!              to observation point followed by a similar computations 
     244      !!              as in geo2ocean. 
     245      !! 
     246      !! ** Action  : Review if there is a better way to do this. 
     247      !! 
     248      !! References :  
     249      !! 
     250      !! History :   
     251      !!      ! :  2009-02 (K. Mogensen) : New routine 
     252      !!---------------------------------------------------------------------- 
     253      !! * Modules used 
     254      !! * Arguments 
     255      TYPE(obs_surf), INTENT(INOUT) :: surfdata    ! Surface data to be read 
     256      INTEGER, INTENT(IN) :: k2dint     ! Horizontal interpolation methed 
     257      REAL(wp), DIMENSION(*) :: & 
     258         & pu, & 
     259         & pv 
     260      !! * Local declarations 
     261      REAL(wp), DIMENSION(2,2,1) :: zweig 
     262      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
     263         & zmasku, & 
     264         & zmaskv, & 
     265         & zcoslu, & 
     266         & zsinlu, & 
     267         & zcoslv, & 
     268         & zsinlv, & 
     269         & zglamu, & 
     270         & zgphiu, & 
     271         & zglamv, & 
     272         & zgphiv 
     273      REAL(wp), DIMENSION(1) :: & 
     274         & zsinu, & 
     275         & zcosu, & 
     276         & zsinv, & 
     277         & zcosv 
     278      REAL(wp) :: zsin 
     279      REAL(wp) :: zcos 
     280      REAL(wp), DIMENSION(1) :: zobsmask 
     281      REAL(wp), POINTER, DIMENSION(:,:) :: zsingu,zcosgu,zsingv,zcosgv 
     282      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
     283         & igrdiu, & 
     284         & igrdju, & 
     285         & igrdiv, & 
     286         & igrdjv 
     287      INTEGER :: ji 
     288      INTEGER :: jk 
     289 
     290      CALL wrk_alloc(jpi,jpj,zsingu,zcosgu,zsingv,zcosgv)  
     291 
     292      !----------------------------------------------------------------------- 
     293      ! Allocate data for message parsing and interpolation 
     294      !----------------------------------------------------------------------- 
     295 
     296      ALLOCATE( & 
     297         & igrdiu(2,2,surfdata%nsurf), & 
     298         & igrdju(2,2,surfdata%nsurf), & 
     299         & zglamu(2,2,surfdata%nsurf), & 
     300         & zgphiu(2,2,surfdata%nsurf), & 
     301         & zmasku(2,2,surfdata%nsurf), & 
     302         & zcoslu(2,2,surfdata%nsurf), & 
     303         & zsinlu(2,2,surfdata%nsurf), & 
     304         & igrdiv(2,2,surfdata%nsurf), & 
     305         & igrdjv(2,2,surfdata%nsurf), & 
     306         & zglamv(2,2,surfdata%nsurf), & 
     307         & zgphiv(2,2,surfdata%nsurf), & 
     308         & zmaskv(2,2,surfdata%nsurf), & 
     309         & zcoslv(2,2,surfdata%nsurf), & 
     310         & zsinlv(2,2,surfdata%nsurf)  & 
     311         & ) 
     312 
     313      !----------------------------------------------------------------------- 
     314      ! Receive the angles on the U and V grids. 
     315      !----------------------------------------------------------------------- 
     316 
     317      CALL obs_rot( zsingu, zcosgu, zsingv, zcosgv ) 
     318 
     319      DO ji = 1, surfdata%nsurf 
     320         igrdiu(1,1,ji) = surfdata%mi(ji,1)-1 
     321         igrdju(1,1,ji) = surfdata%mj(ji,1)-1 
     322         igrdiu(1,2,ji) = surfdata%mi(ji,1)-1 
     323         igrdju(1,2,ji) = surfdata%mj(ji,1) 
     324         igrdiu(2,1,ji) = surfdata%mi(ji,1) 
     325         igrdju(2,1,ji) = surfdata%mj(ji,1)-1 
     326         igrdiu(2,2,ji) = surfdata%mi(ji,1) 
     327         igrdju(2,2,ji) = surfdata%mj(ji,1) 
     328         igrdiv(1,1,ji) = surfdata%mi(ji,2)-1 
     329         igrdjv(1,1,ji) = surfdata%mj(ji,2)-1 
     330         igrdiv(1,2,ji) = surfdata%mi(ji,2)-1 
     331         igrdjv(1,2,ji) = surfdata%mj(ji,2) 
     332         igrdiv(2,1,ji) = surfdata%mi(ji,2) 
     333         igrdjv(2,1,ji) = surfdata%mj(ji,2)-1 
     334         igrdiv(2,2,ji) = surfdata%mi(ji,2) 
     335         igrdjv(2,2,ji) = surfdata%mj(ji,2) 
     336      END DO 
     337 
     338      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, & 
     339         &                  glamu, zglamu ) 
     340      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, & 
     341         &                  gphiu, zgphiu ) 
     342      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, & 
     343         &                  umask(:,:,1), zmasku ) 
     344      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, & 
     345         &                  zsingu, zsinlu ) 
     346      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, & 
     347         &                  zcosgu, zcoslu ) 
     348      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, & 
     349         &                  glamv, zglamv ) 
     350      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, & 
     351         &                  gphiv, zgphiv ) 
     352      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, & 
     353         &                  vmask(:,:,1), zmaskv ) 
     354      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, & 
     355         &                  zsingv, zsinlv ) 
     356      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, & 
     357         &                  zcosgv, zcoslv ) 
     358 
     359      DO ji = 1, surfdata%nsurf 
     360             
     361         CALL obs_int_h2d_init( 1, 1, k2dint, & 
     362            &                   surfdata%rlam(ji), surfdata%rphi(ji), & 
     363            &                   zglamu(:,:,ji), zgphiu(:,:,ji), & 
     364            &                   zmasku(:,:,ji), zweig, zobsmask ) 
     365          
     366         CALL obs_int_h2d( 1, 1, zweig, zsinlu(:,:,ji),  zsinu ) 
     367 
     368         CALL obs_int_h2d( 1, 1, zweig, zcoslu(:,:,ji),  zcosu ) 
     369 
     370         CALL obs_int_h2d_init( 1, 1, k2dint, & 
     371            &                   surfdata%rlam(ji), surfdata%rphi(ji), & 
     372            &                   zglamv(:,:,ji), zgphiv(:,:,ji), & 
     373            &                   zmaskv(:,:,ji), zweig, zobsmask ) 
     374          
     375         CALL obs_int_h2d( 1, 1, zweig, zsinlv(:,:,ji),  zsinv ) 
     376 
     377         CALL obs_int_h2d( 1, 1, zweig, zcoslv(:,:,ji),  zcosv ) 
     378 
     379         ! Assume that the angle at observation point is the  
     380         ! mean of u and v cosines/sines 
     381 
     382         zcos = 0.5_wp * ( zcosu(1) + zcosv(1) ) 
     383         zsin = 0.5_wp * ( zsinu(1) + zsinv(1) ) 
     384 
     385         IF ( ( surfdata%rmod(ji,1) /= fbrmdi ) .AND. & 
     386            & ( surfdata%rmod(ji,2) /= fbrmdi ) ) THEN 
     387            pu(ji) = surfdata%rmod(ji,1) * zcos - & 
     388               &     surfdata%rmod(ji,2) * zsin 
     389            pv(ji) = surfdata%rmod(ji,2) * zcos + & 
     390               &     surfdata%rmod(ji,1) * zsin 
     391         ELSE 
     392            pu(ji) = fbrmdi 
     393            pv(ji) = fbrmdi 
     394         ENDIF 
     395 
     396 
     397      END DO 
     398       
     399      DEALLOCATE( & 
     400         & igrdiu, & 
     401         & igrdju, & 
     402         & zglamu, & 
     403         & zgphiu, & 
     404         & zmasku, & 
     405         & zcoslu, & 
     406         & zsinlu, & 
     407         & igrdiv, & 
     408         & igrdjv, & 
     409         & zglamv, & 
     410         & zgphiv, & 
     411         & zmaskv, & 
     412         & zcoslv, & 
     413         & zsinlv  & 
     414         & ) 
     415 
     416      CALL wrk_dealloc(jpi,jpj,zsingu,zcosgu,zsingv,zcosgv)  
     417 
     418   END SUBROUTINE obs_rotvel_surf 
    231419 
    232420END MODULE obs_rot_vel 

  • branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90

    r11468 r13393  
    5757 
    5858      INTEGER, POINTER, DIMENSION(:) :: & 
    59          & mi,   &        !: i-th grid coord. for interpolating to surface observation 
    60          & mj,   &        !: j-th grid coord. for interpolating to surface observation 
    6159         & mt,   &        !: time record number for gridded data 
    6260         & nsidx,&        !: Surface observation number 
    … …  
    7169         & ntyp           !: Type of surface observation product 
    7270 
     71      INTEGER, POINTER, DIMENSION(:,:) :: & 
     72         & mi,   &        !: i-th grid coord. for interpolating to surface observation 
     73         & mj             !: j-th grid coord. for interpolating to surface observation 
     74 
     75 
    7376      CHARACTER(len=8), POINTER, DIMENSION(:) :: & 
    7477         & cvars          !: Variable names 
    … …  
    9295         & rext           !: Extra fields interpolated to observation points 
    9396 
    94       REAL(KIND=wp), POINTER, DIMENSION(:,:) :: & 
     97      REAL(KIND=wp), POINTER, DIMENSION(:,:,:) :: & 
    9598         & vdmean         !: Time averaged of model field 
    9699 
    … …  
    176179 
    177180      ALLOCATE( & 
    178          & surf%mi(ksurf),      & 
    179          & surf%mj(ksurf),      & 
    180181         & surf%mt(ksurf),      & 
    181182         & surf%nsidx(ksurf),   & 
    … …  
    201202 
    202203      ALLOCATE( &  
     204         & surf%mi(ksurf,kvar),   & 
     205         & surf%mj(ksurf,kvar),   & 
    203206         & surf%robs(ksurf,kvar), & 
    204          & surf%rmod(ksurf,kvar) & 
     207         & surf%rmod(ksurf,kvar)  & 
    205208         & )    
    206209 
    … …  
    230233 
    231234      ALLOCATE( & 
    232          & surf%vdmean(kpi,kpj) & 
     235         & surf%vdmean(kpi,kpj,kvar) & 
    233236         & ) 
    234237 
    … …  
    405408            insurf = insurf + 1 
    406409 
    407             newsurf%mi(insurf)    = surf%mi(ji) 
    408             newsurf%mj(insurf)    = surf%mj(ji) 
     410            newsurf%mi(insurf,:)  = surf%mi(ji,:) 
     411            newsurf%mj(insurf,:)  = surf%mj(ji,:) 
    409412            newsurf%mt(insurf)    = surf%mt(ji) 
    410413            newsurf%nsidx(insurf) = surf%nsidx(ji) 
    … …  
    496499         jj=surf%nsind(ji) 
    497500 
    498          oldsurf%mi(jj)    = surf%mi(ji) 
    499          oldsurf%mj(jj)    = surf%mj(ji) 
     501         oldsurf%mi(jj,:)  = surf%mi(ji,:) 
     502         oldsurf%mj(jj,:)  = surf%mj(ji,:) 
    500503         oldsurf%mt(jj)    = surf%mt(ji) 
    501504         oldsurf%nsidx(jj) = surf%nsidx(ji) 

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

    r12610 r13393  
    437437      CHARACTER(LEN=10) :: clfiletype 
    438438      CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf' 
    439       CHARACTER(LEN=ilenlong) :: cllongname  ! Long name of variable 
    440       CHARACTER(LEN=ilenunit) :: clunits     ! Units of variable 
    441       CHARACTER(LEN=ilengrid) :: clgrid      ! Grid of variable 
     439      CHARACTER(LEN=ilenlong), DIMENSION(surfdata%nvar) :: cllongname  ! Long name of variable 
     440      CHARACTER(LEN=ilenunit), DIMENSION(surfdata%nvar) :: clunits     ! Units of variable 
     441      CHARACTER(LEN=ilengrid), DIMENSION(surfdata%nvar) :: clgrid      ! Grid of variable 
    442442      INTEGER :: jo 
    443443      INTEGER :: ja 
    444444      INTEGER :: je 
     445      INTEGER :: jv 
    445446      INTEGER :: iadd 
    446447      INTEGER :: iext 
    … …  
    522523      CASE('SST') 
    523524 
    524          clfiletype = 'sstfb' 
    525          cllongname = 'Sea surface temperature' 
    526          clunits    = 'Degree centigrade' 
    527          clgrid     = 'T' 
     525         clfiletype    = 'sstfb' 
     526         cllongname(1) = 'Sea surface temperature' 
     527         clunits(1)    = 'Degree centigrade' 
     528         clgrid(1)     = 'T' 
    528529          
    529530      CASE('ICECONC') 
    530531 
    531          clfiletype = 'sicfb' 
    532          cllongname = 'Sea ice concentration' 
    533          clunits    = 'Fraction' 
    534          clgrid     = 'T' 
     532         clfiletype    = 'sicfb' 
     533         cllongname(1) = 'Sea ice concentration' 
     534         clunits(1)    = 'Fraction' 
     535         clgrid(1)     = 'T' 
    535536 
    536537      CASE('SIT') 
    537538 
    538          clfiletype = 'sitfb' 
    539          cllongname = 'Sea ice thickness' 
    540          clunits    = 'm' 
    541          clgrid     = 'T' 
     539         clfiletype    = 'sitfb' 
     540         cllongname(1) = 'Sea ice thickness' 
     541         clunits(1)    = 'm' 
     542         clgrid(1)     = 'T' 
    542543 
    543544      CASE('SSS') 
    544545 
    545          clfiletype = 'sssfb' 
    546          cllongname = 'Sea surface salinity' 
    547          clunits    = 'psu' 
    548          clgrid     = 'T' 
     546         clfiletype    = 'sssfb' 
     547         cllongname(1) = 'Sea surface salinity' 
     548         clunits(1)    = 'psu' 
     549         clgrid(1)     = 'T' 
     550 
     551      CASE('UVEL') 
     552 
     553         clfiletype    = 'ssvfb' 
     554         cllongname(1) = 'Eastward sea surface velocity' 
     555         clunits(1)    = 'm s-1' 
     556         clgrid(1)     = 'U' 
     557         cllongname(2) = 'Northward sea surface velocity' 
     558         clunits(2)    = 'm s-1' 
     559         clgrid(2)     = 'V' 
    549560          
    550561      CASE('SLCHLTOT') 
    551562 
    552          clfiletype = 'slchltotfb' 
    553          cllongname = 'Surface total log10(chlorophyll)' 
    554          clunits    = 'log10(mg/m3)' 
    555          clgrid     = 'T' 
     563         clfiletype    = 'slchltotfb' 
     564         cllongname(1) = 'Surface total log10(chlorophyll)' 
     565         clunits(1)    = 'log10(mg/m3)' 
     566         clgrid(1)     = 'T' 
    556567 
    557568      CASE('SLCHLDIA') 
    558569 
    559          clfiletype = 'slchldiafb' 
    560          cllongname = 'Surface diatom log10(chlorophyll)' 
    561          clunits    = 'log10(mg/m3)' 
    562          clgrid     = 'T' 
     570         clfiletype    = 'slchldiafb' 
     571         cllongname(1) = 'Surface diatom log10(chlorophyll)' 
     572         clunits(1)    = 'log10(mg/m3)' 
     573         clgrid(1)     = 'T' 
    563574 
    564575      CASE('SLCHLNON') 
    565576 
    566          clfiletype = 'slchlnonfb' 
    567          cllongname = 'Surface non-diatom log10(chlorophyll)' 
    568          clunits    = 'log10(mg/m3)' 
    569          clgrid     = 'T' 
     577         clfiletype    = 'slchlnonfb' 
     578         cllongname(1) = 'Surface non-diatom log10(chlorophyll)' 
     579         clunits(1)    = 'log10(mg/m3)' 
     580         clgrid(1)     = 'T' 
    570581 
    571582      CASE('SLCHLDIN') 
    572583 
    573          clfiletype = 'slchldinfb' 
    574          cllongname = 'Surface dinoflagellate log10(chlorophyll)' 
    575          clunits    = 'log10(mg/m3)' 
    576          clgrid     = 'T' 
     584         clfiletype    = 'slchldinfb' 
     585         cllongname(1) = 'Surface dinoflagellate log10(chlorophyll)' 
     586         clunits(1)    = 'log10(mg/m3)' 
     587         clgrid(1)     = 'T' 
    577588 
    578589      CASE('SLCHLMIC') 
    579590 
    580          clfiletype = 'slchlmicfb' 
    581          cllongname = 'Surface microphytoplankton log10(chlorophyll)' 
    582          clunits    = 'log10(mg/m3)' 
    583          clgrid     = 'T' 
     591         clfiletype    = 'slchlmicfb' 
     592         cllongname(1) = 'Surface microphytoplankton log10(chlorophyll)' 
     593         clunits(1)    = 'log10(mg/m3)' 
     594         clgrid(1)     = 'T' 
    584595 
    585596      CASE('SLCHLNAN') 
    586597 
    587          clfiletype = 'slchlnanfb' 
    588          cllongname = 'Surface nanophytoplankton log10(chlorophyll)' 
    589          clunits    = 'log10(mg/m3)' 
    590          clgrid     = 'T' 
     598         clfiletype    = 'slchlnanfb' 
     599         cllongname(1) = 'Surface nanophytoplankton log10(chlorophyll)' 
     600         clunits(1)    = 'log10(mg/m3)' 
     601         clgrid(1)     = 'T' 
    591602 
    592603      CASE('SLCHLPIC') 
    593604 
    594          clfiletype = 'slchlpicfb' 
    595          cllongname = 'Surface picophytoplankton log10(chlorophyll)' 
    596          clunits    = 'log10(mg/m3)' 
    597          clgrid     = 'T' 
     605         clfiletype    = 'slchlpicfb' 
     606         cllongname(1) = 'Surface picophytoplankton log10(chlorophyll)' 
     607         clunits(1)    = 'log10(mg/m3)' 
     608         clgrid(1)     = 'T' 
    598609 
    599610      CASE('SCHLTOT') 
    … …  
    606617      CASE('SLPHYTOT') 
    607618 
    608          clfiletype = 'slphytotfb' 
    609          cllongname = 'Surface total log10(phytoplankton carbon)' 
    610          clunits    = 'log10(mmolC/m3)' 
    611          clgrid     = 'T' 
     619         clfiletype    = 'slphytotfb' 
     620         cllongname(1) = 'Surface total log10(phytoplankton carbon)' 
     621         clunits(1)    = 'log10(mmolC/m3)' 
     622         clgrid(1)     = 'T' 
    612623 
    613624      CASE('SLPHYDIA') 
    614625 
    615          clfiletype = 'slphydiafb' 
    616          cllongname = 'Surface diatom log10(phytoplankton carbon)' 
    617          clunits    = 'log10(mmolC/m3)' 
    618          clgrid     = 'T' 
     626         clfiletype    = 'slphydiafb' 
     627         cllongname(1) = 'Surface diatom log10(phytoplankton carbon)' 
     628         clunits(1)    = 'log10(mmolC/m3)' 
     629         clgrid(1)     = 'T' 
    619630 
    620631      CASE('SLPHYNON') 
    621632 
    622          clfiletype = 'slphynonfb' 
    623          cllongname = 'Surface non-diatom log10(phytoplankton carbon)' 
    624          clunits    = 'log10(mmolC/m3)' 
    625          clgrid     = 'T' 
     633         clfiletype    = 'slphynonfb' 
     634         cllongname(1) = 'Surface non-diatom log10(phytoplankton carbon)' 
     635         clunits(1)    = 'log10(mmolC/m3)' 
     636         clgrid(1)     = 'T' 
    626637 
    627638      CASE('SSPM') 
    628639 
    629          clfiletype = 'sspmfb' 
    630          cllongname = 'Surface suspended particulate matter' 
    631          clunits    = 'g/m3' 
    632          clgrid     = 'T' 
     640         clfiletype    = 'sspmfb' 
     641         cllongname(1) = 'Surface suspended particulate matter' 
     642         clunits(1)    = 'g/m3' 
     643         clgrid(1)     = 'T' 
    633644 
    634645      CASE('SKD490') 
    635646 
    636          clfiletype = 'skd490fb' 
    637          cllongname = 'Surface attenuation coefficient of downwelling radiation at 490 nm' 
    638          clunits    = 'm-1' 
    639          clgrid     = 'T' 
     647         clfiletype    = 'skd490fb' 
     648         cllongname(1) = 'Surface attenuation coefficient of downwelling radiation at 490 nm' 
     649         clunits(1)    = 'm-1' 
     650         clgrid(1)     = 'T' 
    640651 
    641652      CASE('SFCO2') 
    642653 
    643          clfiletype = 'sfco2fb' 
    644          cllongname = 'Surface fugacity of carbon dioxide' 
    645          clunits    = 'uatm' 
    646          clgrid     = 'T' 
     654         clfiletype    = 'sfco2fb' 
     655         cllongname(1) = 'Surface fugacity of carbon dioxide' 
     656         clunits(1)    = 'uatm' 
     657         clgrid(1)     = 'T' 
    647658 
    648659      CASE('SPCO2') 
    649660 
    650          clfiletype = 'spco2fb' 
    651          cllongname = 'Surface partial pressure of carbon dioxide' 
    652          clunits    = 'uatm' 
    653          clgrid     = 'T' 
     661         clfiletype    = 'spco2fb' 
     662         cllongname(1) = 'Surface partial pressure of carbon dioxide' 
     663         clunits(1)    = 'uatm' 
     664         clgrid(1)     = 'T' 
    654665 
    655666      CASE DEFAULT 
    … …  
    664675      IF ( TRIM(surfdata%cvars(1)) /= 'SLA' ) THEN 
    665676       
    666          CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
     677         CALL alloc_obfbdata( fbdata, surfdata%nvar, surfdata%nsurf, 1, & 
    667678            &                 1 + iadd_std + iadd_clm + iadd, iext, .TRUE. ) 
    668679 
    669          fbdata%cname(1)      = surfdata%cvars(1) 
    670          fbdata%coblong(1)    = cllongname 
    671          fbdata%cobunit(1)    = clunits 
     680         DO jv = 1, surfdata%nvar 
     681            fbdata%cname(jv)      = surfdata%cvars(jv) 
     682            fbdata%coblong(jv)    = cllongname(jv) 
     683            fbdata%cobunit(jv)    = clunits(jv) 
     684         END DO 
    672685         DO je = 1, iext 
    673686            fbdata%cextname(je) = pext%cdname(je) 
    … …  
    675688            fbdata%cextunit(je) = pext%cdunit(je,1) 
    676689         END DO 
    677          IF ( TRIM(surfdata%cvars(1)) == 'ICECONC' ) THEN 
    678             fbdata%caddlong(1,1) = 'Model interpolated ICE' 
    679          ELSE 
    680             fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(surfdata%cvars(1)) 
    681          ENDIF 
    682          fbdata%caddunit(1,1) = clunits 
    683          fbdata%cgrid(1)      = clgrid 
     690         DO jv = 1, surfdata%nvar          
     691            IF ( TRIM(surfdata%cvars(1)) == 'ICECONC' ) THEN 
     692               fbdata%caddlong(1,jv) = 'Model interpolated ICE' 
     693            ELSE 
     694               fbdata%caddlong(1,jv) = 'Model interpolated ' // TRIM(surfdata%cvars(jv)) 
     695            ENDIF 
     696            fbdata%caddunit(1,jv) = clunits(jv) 
     697            fbdata%cgrid(jv)      = clgrid(jv) 
     698         END DO             
    684699         DO ja = 1, iadd 
    685700            fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm+ja) = padd%cdname(ja) 
    686             fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm+ja,1) = padd%cdlong(ja,1) 
    687             fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm+ja,1) = padd%cdunit(ja,1) 
    688          END DO 
    689  
    690       ENDIF 
    691        
     701            DO jv = 1, surfdata%nvar                      
     702               fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = padd%cdlong(ja,jv) 
     703               fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = padd%cdunit(ja,jv) 
     704            END DO 
     705         END DO 
     706      ENDIF 
     707 
    692708      fbdata%caddname(1)   = 'Hx' 
    693709      IF ( indx_std /= -1 ) THEN 
    694710         fbdata%caddname(1+iadd_mdt+iadd_std)   = surfdata%cext(indx_std) 
    695          fbdata%caddlong(1+iadd_mdt+iadd_std,1) = 'Obs error standard deviation' 
    696          fbdata%caddunit(1+iadd_mdt+iadd_std,1) = fbdata%cobunit(1) 
     711         DO jv = 1, surfdata%nvar                         
     712            fbdata%caddlong(1+iadd_mdt+iadd_std,1) = 'Obs error standard deviation' 
     713            fbdata%caddunit(1+iadd_mdt+iadd_std,1) = fbdata%cobunit(1) 
     714         END DO 
    697715      ENDIF 
    698716       
    699717      IF ( surfdata%lclim ) THEN 
    700718         fbdata%caddname(1+iadd_mdt+iadd_std+iadd_clm)   = 'CLM' 
    701          fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm,1) = 'Climatology' 
    702          fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm,1) = fbdata%cobunit(1) 
    703       ENDIF 
    704        
    705        
     719         DO jv = 1, surfdata%nvar                                  
     720            fbdata%caddlong(1+iadd_mdt+iadd_std+iadd_clm,jv) = 'Climatology' 
     721            fbdata%caddunit(1+iadd_mdt+iadd_std+iadd_clm,jv) = fbdata%cobunit(1) 
     722         END DO 
     723      ENDIF 
     724 
    706725      WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc 
    707726 
    … …  
    710729         WRITE(numout,*)'obs_wri_surf :' 
    711730         WRITE(numout,*)'~~~~~~~~~~~~~' 
    712          WRITE(numout,*)'Writing '//TRIM(surfdata%cvars(1))//' feedback file : ',TRIM(clfname) 
     731         WRITE(numout,*)'Writing '//TRIM(clfiletype)//' feedback file : ',TRIM(clfname) 
    713732      ENDIF 
    714733 
    … …  
    736755         fbdata%kindex(jo)    = surfdata%nsfil(jo) 
    737756         IF (ln_grid_global) THEN 
    738             fbdata%iobsi(jo,1) = surfdata%mi(jo) 
    739             fbdata%iobsj(jo,1) = surfdata%mj(jo) 
     757            DO jv = 1, surfdata%nvar 
     758               fbdata%iobsi(jo,jv) = surfdata%mi(jo,jv) 
     759               fbdata%iobsj(jo,jv) = surfdata%mj(jo,jv) 
     760            END DO 
    740761         ELSE 
    741             fbdata%iobsi(jo,1) = mig(surfdata%mi(jo)) 
    742             fbdata%iobsj(jo,1) = mjg(surfdata%mj(jo)) 
     762            DO jv = 1, surfdata%nvar          
     763               fbdata%iobsi(jo,jv) = mig(surfdata%mi(jo,jv)) 
     764               fbdata%iobsj(jo,jv) = mjg(surfdata%mj(jo,jv)) 
     765            END DO                
    743766         ENDIF 
    744767         CALL greg2jul( 0, & 
    … …  
    750773            &           fbdata%ptim(jo),   & 
    751774            &           krefdate = 19500101 ) 
    752                      
    753          fbdata%pob(1,jo,1)    = surfdata%robs(jo,1)  
     775 
    754776         fbdata%pdep(1,jo)     = 0.0 
    755777         fbdata%idqc(1,jo)     = 0 
    756778         fbdata%idqcf(:,1,jo)  = 0 
    757          IF ( surfdata%nqc(jo) > 255 ) THEN 
    758             fbdata%ivqc(jo,1)       = 4 
    759             fbdata%ivlqc(1,jo,1)    = 4 
    760             fbdata%ivlqcf(1,1,jo,1) = 0 
    761             fbdata%ivlqcf(2,1,jo,1) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111') 
    762          ELSE 
    763             fbdata%ivqc(jo,1)       = surfdata%nqc(jo) 
    764             fbdata%ivlqc(1,jo,1)    = surfdata%nqc(jo) 
    765             fbdata%ivlqcf(:,1,jo,1) = 0 
    766          ENDIF 
    767          fbdata%iobsk(1,jo,1)  = 0 
    768   
    769          ! Additional variables. 
    770          ! Hx is always the first additional variable 
    771          fbdata%padd(1,jo,1,1) = surfdata%rmod(jo,1) 
    772          ! MDT is output as an additional variable if SLA obs type 
    773          IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) THEN 
    774             fbdata%padd(1,jo,1+iadd_mdt,1) = surfdata%rext(jo,1) 
    775          ENDIF     
    776          ! STD is output as an additional variable if available 
    777          IF ( indx_std /= -1 ) THEN 
    778             fbdata%padd(1,jo,1+iadd_mdt+iadd_std,1) = surfdata%rext(jo,indx_std) 
    779          ENDIF 
    780          ! CLM is output as an additional variable if available 
    781          IF ( surfdata%lclim ) THEN 
    782             fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm,1) = surfdata%rclm(jo,1) 
    783          ENDIF 
    784          ! Then other additional variables are output 
    785          DO ja = 1, iadd 
    786             fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm+ja,1) = & 
    787                & surfdata%rext(jo,padd%ipoint(ja)) 
    788          END DO 
     779         DO jv = 1, surfdata%nvar  
     780            fbdata%pob(1,jo,jv)    = surfdata%robs(jo,jv) 
    789781          
     782            IF ( surfdata%nqc(jo) > 255 ) THEN 
     783               fbdata%ivqc(jo,jv)       = 4 
     784               fbdata%ivlqc(1,jo,jv)    = 4 
     785               fbdata%ivlqcf(1,1,jo,jv) = 0 
     786               fbdata%ivlqcf(2,1,jo,jv) = IAND(surfdata%nqc(jo),b'0000 0000 1111 1111') 
     787            ELSE 
     788               fbdata%ivqc(jo,jv)       = surfdata%nqc(jo) 
     789               fbdata%ivlqc(1,jo,jv)    = surfdata%nqc(jo) 
     790               fbdata%ivlqcf(:,1,jo,jv) = 0 
     791            ENDIF 
     792            fbdata%iobsk(1,jo,jv)  = 0 
     793 
     794           
     795            ! Additional variables. 
     796            ! Hx is always the first additional variable 
     797            fbdata%padd(1,jo,1,jv) = surfdata%rmod(jo,jv) 
     798            ! MDT is output as an additional variable if SLA obs type 
     799            IF ( TRIM(surfdata%cvars(jv)) == 'SLA' ) THEN 
     800               fbdata%padd(1,jo,1+iadd_mdt,jv) = surfdata%rext(jo,1) 
     801            ENDIF     
     802            ! STD is output as an additional variable if available 
     803            IF ( indx_std /= -1 ) THEN 
     804               fbdata%padd(1,jo,1+iadd_mdt+iadd_std,jv) = surfdata%rext(jo,indx_std) 
     805            ENDIF 
     806            ! CLM is output as an additional variable if available 
     807            IF ( surfdata%lclim ) THEN 
     808               fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm,jv) = surfdata%rclm(jo,jv) 
     809            ENDIF 
     810            ! Then other additional variables are output 
     811            DO ja = 1, iadd 
     812               fbdata%padd(1,jo,1+iadd_mdt+iadd_std+iadd_clm+ja,jv) = & 
     813                  & surfdata%rext(jo,padd%ipoint(ja)) 
     814            END DO 
     815         END DO          
    790816         ! Extra variables 
    791817         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%pext(1,jo,1) = surfdata%rext(jo,2)          
Note: See TracChangeset for help on using the changeset viewer.