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 7773 for branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 – NEMO

Ignore:
Timestamp:
2017-03-09T13:52:43+01:00 (7 years ago)
Author:
mattmartin
Message:

Committing updates after doing the following:

  • merging the branch dev_r4650_general_vert_coord_obsoper@7763 into this branch
  • updating it so that the following OBS changes were implemented correctly on top of the simplification changes:
    • generalised vertical coordinate for profile obs. This was done so that is now the default option.
    • sst bias correction implemented with the new simplified obs code.
    • included the biogeochemical obs types int he new simplified obs code.
    • included the changes to exclude obs in the boundary for limited area models
    • included other changes for the efficiency of the obs operator to remove global arrays.
File:
1 edited

Legend:

Unmodified
Added
Removed

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

    • Property svn:keywords deleted
    r5704 r7773  
    4949   !!---------------------------------------------------------------------- 
    5050 
     51   !! * Substitutions  
     52#  include "domzgr_substitute.h90"  
    5153CONTAINS 
    5254 
    5355   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk,          & 
    5456      &                     kit000, kdaystp,                      & 
    55       &                     pvar1, pvar2, pgdept, pmask1, pmask2, & 
     57      &                     pvar1, pvar2, pgdept, pgdepw, 
     58      &                     pmask1, pmask2, & 
    5659      &                     plam1, plam2, pphi1, pphi2,           & 
    5760      &                     k1dint, k2dint, kdailyavtypes ) 
    … …  
    104107      !!      ! 07-03 (K. Mogensen) General handling of profiles 
    105108      !!      ! 15-02 (M. Martin) Combined routine for all profile types 
     109      !!      ! 17-02 (M. Martin) Include generalised vertical coordinate changes 
    106110      !!----------------------------------------------------------------------- 
    107111 
    … …  
    133137         & pphi1,    &               ! Model latitudes for variable 1 
    134138         & pphi2                     ! Model latitudes for variable 2 
    135       REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 
    136          & pgdept                    ! Model array of depth levels 
     139      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &  
     140         & pgdept, &                 ! Model array of depth T levels  
     141         & pgdepw                    ! Model array of depth W levels  
    137142      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
    138143         & kdailyavtypes             ! Types for daily averages 
    … …  
    164169         & zobsk,    & 
    165170         & zobs2k 
    166       REAL(KIND=wp), DIMENSION(2,2,kpk) :: & 
     171      REAL(KIND=wp), DIMENSION(2,2,1) :: & 
    167172         & zweig1, & 
    168          & zweig2 
     173         & zweig2, & 
     174         & zweig 
    169175      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 
    170176         & zmask1, & 
    171177         & zmask2, & 
    172          & zint1, & 
    173          & zint2, & 
    174          & zinm1, & 
    175          & zinm2 
     178         & zint1,  & 
     179         & zint2,  & 
     180         & zinm1,  & 
     181         & zinm2,  & 
     182         & zgdept, &  
     183         & zgdepw 
    176184      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    177185         & zglam1, & 
    … …  
    179187         & zgphi1, & 
    180188         & zgphi2 
     189      REAL(KIND=wp), DIMENSION(1) :: zmsk_1, zmsk_2    
     190      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner        
     191 
    181192      LOGICAL :: ld_dailyav 
    182193 
    … …  
    259270         & zmask1(2,2,kpk,ipro),  & 
    260271         & zmask2(2,2,kpk,ipro),  & 
    261          & zint1(2,2,kpk,ipro),  & 
    262          & zint2(2,2,kpk,ipro)   & 
     272         & zint1(2,2,kpk,ipro),   & 
     273         & zint2(2,2,kpk,ipro),   & 
     274         & zgdept(2,2,kpk,ipro),  &  
     275         & zgdepw(2,2,kpk,ipro)   &  
    263276         & ) 
    264277 
    … …  
    283296      END DO 
    284297 
     298      ! Initialise depth arrays 
     299      zgdept(:,:,:,:) = 0.0 
     300      zgdepw(:,:,:,:) = 0.0 
     301 
    285302      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 ) 
    286303      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 ) 
    … …  
    293310      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2,   zint2 ) 
    294311 
     312      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept )  
     313      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdepw, zgdepw )  
     314 
    295315      ! At the end of the day also get interpolated means 
    296316      IF ( ld_dailyav .AND. idayend == 0 ) THEN 
    … …  
    307327 
    308328      ENDIF 
     329 
     330      ! Return if no observations to process  
     331      ! Has to be done after comm commands to ensure processors  
     332      ! stay in sync  
     333      IF ( ipro == 0 ) RETURN  
    309334 
    310335      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
    … …  
    332357         zphi = prodatqc%rphi(jobs) 
    333358 
    334          ! Horizontal weights and vertical mask 
    335  
     359         ! Horizontal weights  
     360         ! Masked values are calculated later.   
    336361         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 
    337362 
    338             CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
     363            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     & 
    339364               &                   zglam1(:,:,iobs), zgphi1(:,:,iobs), & 
    340                &                   zmask1(:,:,:,iobs), zweig1, zobsmask1 ) 
     365               &                   zmask1(:,:,1,iobs), zweig1, zmsk_1 ) 
    341366 
    342367         ENDIF 
    … …  
    344369         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    345370 
    346             CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
     371            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     & 
    347372               &                   zglam2(:,:,iobs), zgphi2(:,:,iobs), & 
    348                &                   zmask2(:,:,:,iobs), zweig2, zobsmask2 ) 
     373               &                   zmask2(:,:,1,iobs), zweig2, zmsk_2 ) 
    349374  
    350375         ENDIF 
    … …  
    358383               IF ( idayend == 0 )  THEN 
    359384                  ! Daily averaged data 
    360                   CALL obs_int_h2d( kpk, kpk,      & 
    361                      &              zweig1, zinm1(:,:,:,iobs), zobsk ) 
    362  
    363                ENDIF 
    364  
    365             ELSE  
    366  
    367                ! Point data 
    368                CALL obs_int_h2d( kpk, kpk,      & 
    369                   &              zweig1, zint1(:,:,:,iobs), zobsk ) 
    370  
    371             ENDIF 
    372  
    373             !------------------------------------------------------------- 
    374             ! Compute vertical second-derivative of the interpolating  
    375             ! polynomial at obs points 
    376             !------------------------------------------------------------- 
    377  
    378             IF ( k1dint == 1 ) THEN 
    379                CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   & 
    380                   &                  pgdept, zobsmask1 ) 
    381             ENDIF 
    382  
    383             !----------------------------------------------------------------- 
    384             !  Vertical interpolation to the observation point 
    385             !----------------------------------------------------------------- 
    386             ista = prodatqc%npvsta(jobs,1) 
    387             iend = prodatqc%npvend(jobs,1) 
    388             CALL obs_int_z1d( kpk,                & 
    389                & prodatqc%var(1)%mvk(ista:iend),  & 
    390                & k1dint, iend - ista + 1,         & 
    391                & prodatqc%var(1)%vdep(ista:iend), & 
    392                & zobsk, zobs2k,                   & 
    393                & prodatqc%var(1)%vmod(ista:iend), & 
    394                & pgdept, zobsmask1 ) 
    395  
    396          ENDIF 
    397  
     385 
     386                  ! vertically interpolate all 4 corners  
     387                  ista = prodatqc%npvsta(jobs,1)  
     388                  iend = prodatqc%npvend(jobs,1)  
     389                  inum_obs = iend - ista + 1  
     390                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
     391 
     392                  DO iin=1,2  
     393                     DO ijn=1,2  
     394 
     395                        IF ( k1dint == 1 ) THEN  
     396                           CALL obs_int_z1d_spl( kpk, &  
     397                              &     zinm1(iin,ijn,:,iobs), &  
     398                              &     zobs2k, zgdept(iin,ijn,:,iobs), &  
     399                              &     zmask1(iin,ijn,:,iobs))  
     400                        ENDIF  
     401        
     402                        CALL obs_level_search(kpk, &  
     403                           &    zgdept(iin,ijn,:,iobs), &  
     404                           &    inum_obs, prodatqc%var(1)%vdep(ista:iend), &  
     405                           &    iv_indic)  
     406 
     407                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
     408                           &    prodatqc%var(1)%vdep(ista:iend), &  
     409                           &    zinm1(iin,ijn,:,iobs), &  
     410                           &    zobs2k, interp_corner(iin,ijn,:), &  
     411                           &    zgdept(iin,ijn,:,iobs), &  
     412                           &    zmask1(iin,ijn,:,iobs))  
     413        
     414                     ENDDO  
     415                  ENDDO  
     416 
     417               ENDIF !idayend 
     418 
     419            ELSE    
     420 
     421               ! Point data  
     422      
     423               ! vertically interpolate all 4 corners  
     424               ista = prodatqc%npvsta(jobs,1)  
     425               iend = prodatqc%npvend(jobs,1)  
     426               inum_obs = iend - ista + 1  
     427               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
     428               DO iin=1,2   
     429                  DO ijn=1,2  
     430                     
     431                     IF ( k1dint == 1 ) THEN  
     432                        CALL obs_int_z1d_spl( kpk, &  
     433                           &    zint1(iin,ijn,:,iobs),&  
     434                           &    zobs2k, zgdept(iin,ijn,:,iobs), &  
     435                           &    zmask1(iin,ijn,:,iobs))  
     436   
     437                     ENDIF  
     438        
     439                     CALL obs_level_search(kpk, &  
     440                         &        zgdept(iin,ijn,:,iobs),&  
     441                         &        inum_obs, prodatqc%var(1)%vdep(ista:iend), &  
     442                         &        iv_indic)  
     443 
     444                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
     445                         &          prodatqc%var(1)%vdep(ista:iend),     &  
     446                         &          zint1(iin,ijn,:,iobs),            &  
     447                         &          zobs2k,interp_corner(iin,ijn,:), &  
     448                         &          zgdept(iin,ijn,:,iobs),         &  
     449                         &          zmask1(iin,ijn,:,iobs) )       
     450          
     451                  ENDDO  
     452               ENDDO  
     453              
     454            ENDIF  
     455 
     456            !-------------------------------------------------------------  
     457            ! Compute the horizontal interpolation for every profile level  
     458            !-------------------------------------------------------------  
     459              
     460            DO ikn=1,inum_obs  
     461               iend=ista+ikn-1 
     462                   
     463               zweig(:,:,1) = 0._wp  
     464    
     465               ! This code forces the horizontal weights to be   
     466               ! zero IF the observation is below the bottom of the   
     467               ! corners of the interpolation nodes, Or if it is in   
     468               ! the mask. This is important for observations near   
     469               ! steep bathymetry  
     470               DO iin=1,2  
     471                  DO ijn=1,2  
     472      
     473                     depth_loop1: DO ik=kpk,2,-1  
     474                        IF(zmask1(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
     475                             
     476                           zweig(iin,ijn,1) = &   
     477                              & zweig1(iin,ijn,1) * &  
     478                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
     479                              &  - prodatqc%var(1)%vdep(iend)),0._wp)  
     480                             
     481                           EXIT depth_loop1  
     482 
     483                        ENDIF  
     484 
     485                     ENDDO depth_loop1  
     486      
     487                  ENDDO  
     488               ENDDO  
     489    
     490               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &  
     491                  &              prodatqc%var(1)%vmod(iend:iend) )  
     492 
     493                  ! Set QC flag for any observations found below the bottom 
     494                  ! needed as the check here is more strict than that in obs_prep 
     495               IF (sum(zweig) == 0.0_wp) prodatqc%var(1)%nvqc(iend:iend)=4 
     496  
     497            ENDDO  
     498  
     499            DEALLOCATE(interp_corner,iv_indic)  
     500           
     501         ENDIF  
     502 
     503         ! For the second variable 
    398504         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    399505 
    … …  
    403509 
    404510               IF ( idayend == 0 )  THEN 
    405  
    406511                  ! Daily averaged data 
    407                   CALL obs_int_h2d( kpk, kpk,      & 
    408                      &              zweig2, zinm2(:,:,:,iobs), zobsk ) 
    409  
    410                ENDIF 
    411  
    412             ELSE 
    413  
    414                ! Point data 
    415                CALL obs_int_h2d( kpk, kpk,      & 
    416                   &              zweig2, zint2(:,:,:,iobs), zobsk ) 
    417  
    418             ENDIF 
    419  
    420  
    421             !------------------------------------------------------------- 
    422             ! Compute vertical second-derivative of the interpolating  
    423             ! polynomial at obs points 
    424             !------------------------------------------------------------- 
    425  
    426             IF ( k1dint == 1 ) THEN 
    427                CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 
    428                   &                  pgdept, zobsmask2 ) 
    429             ENDIF 
    430  
    431             !---------------------------------------------------------------- 
    432             !  Vertical interpolation to the observation point 
    433             !---------------------------------------------------------------- 
    434             ista = prodatqc%npvsta(jobs,2) 
    435             iend = prodatqc%npvend(jobs,2) 
    436             CALL obs_int_z1d( kpk, & 
    437                & prodatqc%var(2)%mvk(ista:iend),& 
    438                & k1dint, iend - ista + 1, & 
    439                & prodatqc%var(2)%vdep(ista:iend),& 
    440                & zobsk, zobs2k, & 
    441                & prodatqc%var(2)%vmod(ista:iend),& 
    442                & pgdept, zobsmask2 ) 
    443  
    444          ENDIF 
    445  
    446       END DO 
     512 
     513                  ! vertically interpolate all 4 corners  
     514                  ista = prodatqc%npvsta(jobs,2)  
     515                  iend = prodatqc%npvend(jobs,2)  
     516                  inum_obs = iend - ista + 1  
     517                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
     518 
     519                  DO iin=1,2  
     520                     DO ijn=1,2  
     521 
     522                        IF ( k1dint == 1 ) THEN  
     523                           CALL obs_int_z1d_spl( kpk, &  
     524                              &     zinm2(iin,ijn,:,iobs), &  
     525                              &     zobs2k, zgdept(iin,ijn,:,iobs), &  
     526                              &     zmask2(iin,ijn,:,iobs))  
     527                        ENDIF  
     528        
     529                        CALL obs_level_search(kpk, &  
     530                           &    zgdept(iin,ijn,:,iobs), &  
     531                           &    inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
     532                           &    iv_indic)  
     533 
     534                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
     535                           &    prodatqc%var(2)%vdep(ista:iend), &  
     536                           &    zinm2(iin,ijn,:,iobs), &  
     537                           &    zobs2k, interp_corner(iin,ijn,:), &  
     538                           &    zgdept(iin,ijn,:,iobs), &  
     539                           &    zmask2(iin,ijn,:,iobs))  
     540        
     541                     ENDDO  
     542                  ENDDO  
     543 
     544               ENDIF !idayend 
     545 
     546            ELSE    
     547 
     548               ! Point data  
     549      
     550               ! vertically interpolate all 4 corners  
     551               ista = prodatqc%npvsta(jobs,2)  
     552               iend = prodatqc%npvend(jobs,2)  
     553               inum_obs = iend - ista + 1  
     554               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
     555               DO iin=1,2   
     556                  DO ijn=1,2  
     557                     
     558                     IF ( k1dint == 1 ) THEN  
     559                        CALL obs_int_z1d_spl( kpk, &  
     560                           &    zint2(iin,ijn,:,iobs),&  
     561                           &    zobs2k, zgdept(iin,ijn,:,iobs), &  
     562                           &    zmask2(iin,ijn,:,iobs))  
     563   
     564                     ENDIF  
     565        
     566                     CALL obs_level_search(kpk, &  
     567                         &        zgdept(iin,ijn,:,iobs),&  
     568                         &        inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
     569                         &        iv_indic)  
     570 
     571                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
     572                         &          prodatqc%var(2)%vdep(ista:iend),     &  
     573                         &          zint2(iin,ijn,:,iobs),            &  
     574                         &          zobs2k,interp_corner(iin,ijn,:), &  
     575                         &          zgdept(iin,ijn,:,iobs),         &  
     576                         &          zmask2(iin,ijn,:,iobs) )       
     577          
     578                  ENDDO  
     579               ENDDO  
     580              
     581            ENDIF  
     582 
     583            !-------------------------------------------------------------  
     584            ! Compute the horizontal interpolation for every profile level  
     585            !-------------------------------------------------------------  
     586              
     587            DO ikn=1,inum_obs  
     588               iend=ista+ikn-1 
     589                   
     590               zweig(:,:,1) = 0._wp  
     591    
     592               ! This code forces the horizontal weights to be   
     593               ! zero IF the observation is below the bottom of the   
     594               ! corners of the interpolation nodes, Or if it is in   
     595               ! the mask. This is important for observations near   
     596               ! steep bathymetry  
     597               DO iin=1,2  
     598                  DO ijn=1,2  
     599      
     600                     depth_loop2: DO ik=kpk,2,-1  
     601                        IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
     602                             
     603                           zweig(iin,ijn,1) = &   
     604                              & zweig2(iin,ijn,1) * &  
     605                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
     606                              &  - prodatqc%var(2)%vdep(iend)),0._wp)  
     607                             
     608                           EXIT depth_loop2  
     609 
     610                        ENDIF  
     611 
     612                     ENDDO depth_loop2  
     613      
     614                  ENDDO  
     615               ENDDO  
     616    
     617               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &  
     618                  &              prodatqc%var(2)%vmod(iend:iend) )  
     619 
     620                  ! Set QC flag for any observations found below the bottom 
     621                  ! needed as the check here is more strict than that in obs_prep 
     622               IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4 
     623  
     624            ENDDO  
     625  
     626            DEALLOCATE(interp_corner,iv_indic)  
     627           
     628         ENDIF  
    447629 
    448630      ! Deallocate the data for interpolation 
    … …  
    459641         & zmask2, & 
    460642         & zint1,  & 
    461          & zint2   & 
     643         & zint2,  & 
     644         & zgdept, & 
     645         & zgdepw  & 
    462646         & ) 
    463647 
Note: See TracChangeset for help on using the changeset viewer.