Changeset 4746


Ignore:
Timestamp:
2014-08-15T11:26:21+02:00 (6 years ago)
Author:
jwhile
Message:

Updated to allow interpolation in s-coordinates

Location:
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/OBS
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90

    r4624 r4746  
    118118   !!---------------------------------------------------------------------- 
    119119 
     120   !! * Substitutions  
     121#  include "domzgr_substitute.h90" 
    120122CONTAINS 
    121123 
     
    207209      ! Read namelist parameters 
    208210      !----------------------------------------------------------------------- 
    209  
     211       
     212      !Initalise all values in namelist arrays 
     213      enactfiles(:) = '' 
     214      coriofiles(:) = '' 
     215      profbfiles(:) = '' 
     216      slafilesact(:) = '' 
     217      slafilespas(:) = '' 
     218      slafbfiles(:) = '' 
     219      sstfiles(:)   = '' 
     220      sstfbfiles(:) = '' 
     221      seaicefiles(:) = '' 
    210222      velcurfiles(:) = '' 
    211223      veladcpfiles(:) = '' 
     224      velavcurfiles(:) = '' 
     225      velhrcurfiles(:) = '' 
     226      velavadcpfiles(:) = '' 
     227      velhradcpfiles(:) = '' 
     228      velfbfiles(:) = '' 
     229      velcurfiles(:) = '' 
     230      veladcpfiles(:) = '' 
     231      endailyavtypes(:) = -1 
     232      endailyavtypes(1) = 820 
     233      ln_profb_ena(:) = .FALSE. 
     234      ln_profb_enatim(:) = .TRUE. 
     235      ln_velfb_av(:) = .FALSE. 
     236      ln_ignmis = .FALSE.       
    212237      CALL ini_date( dobsini ) 
    213238      CALL fin_date( dobsend ) 
     
    972997      !!        !  07-04  (G. Smith) Generalized surface operators 
    973998      !!        !  08-10  (M. Valdivieso) obs operator for velocity profiles 
     999      !!        !  14-08  (J. While) observation operator for profiles in  
     1000      !!                             generalised vertical coordinates 
    9741001      !!---------------------------------------------------------------------- 
    9751002      !! * Modules used 
     
    9771004         & rdt,           &                        
    9781005         & gdept_1d,       &              
     1006#if defined key_vvl  
     1007         & gdept_n,       &  
     1008#else  
     1009         & gdept_1d,         &  
     1010#endif                                         
    9791011         & tmask, umask, vmask                             
    9801012      USE phycst, ONLY : &              ! Physical constants 
     
    10341066      IF ( ln_t3d .OR. ln_s3d ) THEN 
    10351067         DO jprofset = 1, nprofsets 
    1036             IF ( ld_enact(jprofset) ) THEN 
    1037                CALL obs_pro_opt( prodatqc(jprofset),                     & 
    1038                   &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
    1039                   &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
    1040                   &              gdept_1d, tmask, n1dint, n2dint,        & 
    1041                   &              kdailyavtypes = endailyavtypes ) 
     1068            IF( ln_zco .OR. ln_zps ) THEN  
     1069               IF ( ld_enact(jprofset) ) THEN  
     1070                  CALL obs_pro_opt( prodatqc(jprofset),                     &  
     1071                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   &  
     1072                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   &  
     1073                     &              gdept_1d, tmask, n1dint, n2dint,         &  
     1074                     &              kdailyavtypes = endailyavtypes )  
     1075               ELSE  
     1076                  CALL obs_pro_opt( prodatqc(jprofset),                     &  
     1077                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   &  
     1078                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   &  
     1079                     &              gdept_1d, tmask, n1dint, n2dint               )  
     1080               ENDIF  
    10421081            ELSE 
    1043                CALL obs_pro_opt( prodatqc(jprofset),                     & 
    1044                   &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
    1045                   &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
    1046                   &              gdept_1d, tmask, n1dint, n2dint              ) 
     1082               IF ( ld_enact(jprofset) ) THEN  
     1083                  CALL obs_pro_sco_opt( prodatqc(jprofset),                 &  
     1084                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   &  
     1085                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   &  
     1086                     &              fsdept(:,:,:), tmask, n1dint, n2dint,   &  
     1087                     &              kdailyavtypes = endailyavtypes )  
     1088               ELSE  
     1089                  CALL obs_pro_sco_opt( prodatqc(jprofset),                 &  
     1090                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   &  
     1091                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   &  
     1092                     &              fsdept(:,:,:), tmask, n1dint, n2dint )  
     1093               ENDIF  
    10471094            ENDIF 
    10481095         END DO 
  • branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

    r4245 r4746  
    99   !!   obs_pro_opt :    Compute the model counterpart of temperature and 
    1010   !!                    salinity observations from profiles 
     11   !!   obs_pro_sco_opt: Compute the model counterpart of temperature and  
     12   !!                    salinity observations from profiles in generalised  
     13   !!                    vertical coordinates  
    1114   !!   obs_sla_opt :    Compute the model counterpart of sea level anomaly 
    1215   !!                    observations 
     
    3740   USE dom_oce,       ONLY : & 
    3841      & glamt, glamu, glamv, & 
    39       & gphit, gphiu, gphiv 
     42      & gphit, gphiu, gphiv, &  
     43#if defined key_vvl  
     44      & gdept_n  
     45#else  
     46      & gdept_0  
     47#endif   
    4048   USE lib_mpp,       ONLY : & 
    4149      & ctl_warn, ctl_stop 
    42  
     50   USE obs_grid,      ONLY : &  
     51      & obs_level_search      
     52       
    4353   IMPLICIT NONE 
    4454 
     
    4757 
    4858   PUBLIC obs_pro_opt, &  ! Compute the model counterpart of profile observations 
     59      &   obs_pro_sco_opt, &  ! Compute the model counterpart of profile observations  
     60                              ! in generalised vertical coordinates  
    4961      &   obs_sla_opt, &  ! Compute the model counterpart of SLA observations 
    5062      &   obs_sst_opt, &  ! Compute the model counterpart of SST observations 
     
    6173   !!---------------------------------------------------------------------- 
    6274 
     75   !! * Substitutions  
     76#  include "domzgr_substitute.h90"  
    6377CONTAINS 
    6478 
     
    449463   END SUBROUTINE obs_pro_opt 
    450464 
     465   SUBROUTINE obs_pro_sco_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &  
     466      &                    ptn, psn, pgdept, ptmask, k1dint, k2dint, &  
     467      &                    kdailyavtypes )  
     468      !!-----------------------------------------------------------------------  
     469      !!  
     470      !!                     ***  ROUTINE obs_pro_opt  ***  
     471      !!  
     472      !! ** Purpose : Compute the model counterpart of profiles  
     473      !!              data by interpolating from the model grid to the   
     474      !!              observation point. Generalised vertical coordinate version  
     475      !!  
     476      !! ** Method  : Linearly interpolate to each observation point using   
     477      !!              the model values at the corners of the surrounding grid box.  
     478      !!  
     479      !!          First, model values on the model grid are interpolated vertically to the  
     480      !!          Depths of the profile observations.  Two vertical interpolation schemes are  
     481      !!          available:  
     482      !!          - linear       (k1dint = 0)  
     483      !!          - Cubic spline (k1dint = 1)     
     484      !!  
     485      !!  
     486      !!         Secondly the interpolated values are interpolated horizontally to the   
     487      !!         obs (lon, lat) point.  
     488      !!         Several horizontal interpolation schemes are available:  
     489      !!        - distance-weighted (great circle) (k2dint = 0)  
     490      !!        - distance-weighted (small angle)  (k2dint = 1)  
     491      !!        - bilinear (geographical grid)     (k2dint = 2)  
     492      !!        - bilinear (quadrilateral grid)    (k2dint = 3)  
     493      !!        - polynomial (quadrilateral grid)  (k2dint = 4)  
     494      !!  
     495      !!    For the cubic spline the 2nd derivative of the interpolating   
     496      !!    polynomial is computed before entering the vertical interpolation   
     497      !!    routine.  
     498      !!  
     499      !!    For ENACT moored buoy data (e.g., TAO), the model equivalent is  
     500      !!    a daily mean model temperature field. So, we first compute  
     501      !!    the mean, then interpolate only at the end of the day.  
     502      !!  
     503      !!    This is the procedure to be used with generalised vertical model   
     504      !!    coordinates (ie s-coordinates. It is ~4x slower than the equivalent  
     505      !!    horizontal then vertical interpolation algorithm, but can deal with situations  
     506      !!    where the model levels are not flat.  
     507      !!    ONLY PERFORMED if ln_sco=.TRUE.   
     508      !!        
     509      !!    Note: the in situ temperature observations must be converted  
     510      !!    to potential temperature (the model variable) prior to  
     511      !!    assimilation.   
     512      !!??????????????????????????????????????????????????????????????  
     513      !!    INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???  
     514      !!??????????????????????????????????????????????????????????????  
     515      !!  
     516      !! ** Action  :  
     517      !!  
     518      !! History :  
     519      !!      ! 2014-08 (J. While) Adapted from obs_pro_opt to handel generalised  
     520      !!                           vertical coordinates 
     521      !!-----------------------------------------------------------------------  
     522    
     523      !! * Modules used  
     524      USE obs_profiles_def   ! Definition of storage space for profile obs.  
     525      USE dom_oce,  ONLY : &  
     526#if defined key_vvl   
     527      gdepw_n  
     528#else  
     529      gdepw_0  
     530#endif  
     531        
     532      IMPLICIT NONE  
     533  
     534      !! * Arguments  
     535      TYPE(obs_prof), INTENT(INOUT) :: prodatqc   ! Subset of profile data not failing screening  
     536      INTEGER, INTENT(IN) :: kt        ! Time step  
     537      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters  
     538      INTEGER, INTENT(IN) :: kpj  
     539      INTEGER, INTENT(IN) :: kpk  
     540      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step   
     541                                       !   (kit000-1 = restart time)  
     542      INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header)  
     543      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)  
     544      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                      
     545      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &  
     546         & ptn,    &    ! Model temperature field  
     547         & psn,    &    ! Model salinity field  
     548         & ptmask       ! Land-sea mask  
     549      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,jpj,kpk) :: &  
     550         & pgdept       ! Model array of depth levels  
     551      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &  
     552         & kdailyavtypes   ! Types for daily averages  
     553      !! * Local declarations  
     554      INTEGER ::   ji  
     555      INTEGER ::   jj  
     556      INTEGER ::   jk  
     557      INTEGER ::   iico, ijco  
     558      INTEGER ::   jobs  
     559      INTEGER ::   inrc  
     560      INTEGER ::   ipro  
     561      INTEGER ::   idayend  
     562      INTEGER ::   ista  
     563      INTEGER ::   iend  
     564      INTEGER ::   iobs  
     565      INTEGER ::   iin, ijn, ikn, ik   ! looping indices over interpolation nodes  
     566      INTEGER, DIMENSION(imaxavtypes) :: &  
     567         & idailyavtypes  
     568      REAL(KIND=wp) :: zlam  
     569      REAL(KIND=wp) :: zphi  
     570      REAL(KIND=wp) :: zdaystp  
     571      REAL(KIND=wp), DIMENSION(kpk) :: &  
     572         & zobsmask, &  
     573         & zobsk,    &  
     574         & zobs2k  
     575      REAL(KIND=wp), DIMENSION(2,2,1) :: &  
     576         & zweig, &  
     577         & l_zweig  
     578      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &  
     579         & zmask, &  
     580         & zintt, &  
     581         & zints, &  
     582         & zinmt, &  
     583         & zgdept,&  
     584         & zgdepw,&  
     585         & zinms  
     586      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &  
     587         & zglam, &  
     588         & zgphi  
     589      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &  
     590         & igrdi, &  
     591         & igrdj  
     592      INTEGER :: &  
     593         & inum_obs     
     594      REAL(KIND=wp), DIMENSION(1) :: zmsk_1        
     595      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner        
     596      INTEGER, ALLOCATABLE, DIMENSION(:) ::           v_indic     
     597  
     598      !------------------------------------------------------------------------  
     599      ! Local initialization   
     600      !------------------------------------------------------------------------  
     601      ! ... Record and data counters  
     602      inrc = kt - kit000 + 2  
     603      ipro = prodatqc%npstp(inrc)  
     604   
     605      ! Daily average types  
     606      IF ( PRESENT(kdailyavtypes) ) THEN  
     607         idailyavtypes(:) = kdailyavtypes(:)  
     608      ELSE  
     609         idailyavtypes(:) = -1  
     610      ENDIF  
     611  
     612      ! Initialize daily mean for first time-step  
     613      idayend = MOD( kt - kit000 + 1, kdaystp )  
     614  
     615      ! Added kt == 0 test to catch restart case   
     616      IF ( idayend == 1 .OR. kt == 0) THEN  
     617           
     618         IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt  
     619         DO jk = 1, jpk  
     620            DO jj = 1, jpj  
     621               DO ji = 1, jpi  
     622                  prodatqc%vdmean(ji,jj,jk,1) = 0.0  
     623                  prodatqc%vdmean(ji,jj,jk,2) = 0.0  
     624               END DO  
     625            END DO  
     626         END DO  
     627        
     628      ENDIF  
     629        
     630      DO jk = 1, jpk  
     631         DO jj = 1, jpj  
     632            DO ji = 1, jpi  
     633               ! Increment the temperature field for computing daily mean  
     634               prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &  
     635               &                        + ptn(ji,jj,jk)  
     636               ! Increment the salinity field for computing daily mean  
     637               prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &  
     638               &                        + psn(ji,jj,jk)  
     639            END DO  
     640         END DO  
     641      END DO  
     642     
     643      ! Compute the daily mean at the end of day  
     644      zdaystp = 1.0 / REAL( kdaystp )  
     645      IF ( idayend == 0 ) THEN  
     646         DO jk = 1, jpk  
     647            DO jj = 1, jpj  
     648               DO ji = 1, jpi  
     649                  prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &  
     650                  &                        * zdaystp  
     651                  prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &  
     652                  &                           * zdaystp  
     653               END DO  
     654            END DO  
     655         END DO  
     656      ENDIF  
     657  
     658      ! Get the data for interpolation  
     659      ALLOCATE( &  
     660      & igrdi(2,2,ipro),      &  
     661      & igrdj(2,2,ipro),      &  
     662      & zglam(2,2,ipro),      &  
     663      & zgphi(2,2,ipro),      &  
     664      & zmask(2,2,kpk,ipro),  &  
     665      & zintt(2,2,kpk,ipro),  &  
     666      & zints(2,2,kpk,ipro),  &  
     667      & zgdept(2,2,kpk,ipro), &  
     668      & zgdepw(2,2,kpk,ipro)  &  
     669      & )  
     670  
     671      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro  
     672         iobs = jobs - prodatqc%nprofup  
     673         igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1  
     674         igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1  
     675         igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1  
     676         igrdj(1,2,iobs) = prodatqc%mj(jobs,1)  
     677         igrdi(2,1,iobs) = prodatqc%mi(jobs,1)  
     678         igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1  
     679         igrdi(2,2,iobs) = prodatqc%mi(jobs,1)  
     680         igrdj(2,2,iobs) = prodatqc%mj(jobs,1)  
     681      END DO  
     682  
     683      CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam )  
     684      CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi )  
     685      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask )  
     686      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn,   zintt )  
     687      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn,   zints )  
     688      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, fsdept(:,:,:), &  
     689        &                     zgdept )  
     690      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, fsdepw(:,:,:), &  
     691        &                     zgdepw )  
     692  
     693      ! At the end of the day also get interpolated means  
     694      IF ( idayend == 0 ) THEN  
     695  
     696         ALLOCATE( &  
     697         & zinmt(2,2,kpk,ipro),  &  
     698         & zinms(2,2,kpk,ipro)   &  
     699         & )  
     700  
     701         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &  
     702         &                  prodatqc%vdmean(:,:,:,1), zinmt )  
     703         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &  
     704         &                  prodatqc%vdmean(:,:,:,2), zinms )  
     705  
     706      ENDIF  
     707        
     708      ! Return if no observations to process  
     709      ! Has to be done after comm commands to ensure processors  
     710      ! stay in sync  
     711      IF ( ipro == 0 ) RETURN  
     712  
     713      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro  
     714     
     715         iobs = jobs - prodatqc%nprofup  
     716     
     717         IF ( kt /= prodatqc%mstp(jobs) ) THEN  
     718              
     719            IF(lwp) THEN  
     720               WRITE(numout,*)  
     721               WRITE(numout,*) ' E R R O R : Observation',              &  
     722                  &            ' time step is not consistent with the', &  
     723                  &            ' model time step'  
     724               WRITE(numout,*) ' ========='  
     725               WRITE(numout,*)  
     726               WRITE(numout,*) ' Record  = ', jobs,                    &  
     727                  &            ' kt      = ', kt,                      &  
     728                  &            ' mstp    = ', prodatqc%mstp(jobs), &  
     729                  &            ' ntyp    = ', prodatqc%ntyp(jobs)  
     730            ENDIF  
     731            CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )  
     732         ENDIF  
     733           
     734         zlam = prodatqc%rlam(jobs)  
     735         zphi = prodatqc%rphi(jobs)  
     736           
     737         ! Horizontal weights  
     738         ! Only calculated once, for both T and S.  
     739         ! Masked values are calculated later.   
     740  
     741         IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. &  
     742         & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN  
     743  
     744            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     &  
     745            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &  
     746            &                   zmask(:,:,1,iobs), zweig, zmsk_1 )  
     747  
     748         ENDIF  
     749          
     750         ! IF zmsk_1 = 0; then ob is on land  
     751         IF (zmsk_1(1) < 0.1) THEN  
     752            WRITE(numout,*) 'WARNING (obs_oper) :- profile found within landmask'  
     753    
     754         ELSE   
     755              
     756            ! Temperature  
     757              
     758            IF ( prodatqc%npvend(jobs,1) > 0 ) THEN   
     759     
     760               zobsk(:) = obfillflt  
     761     
     762               IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN  
     763     
     764                  IF ( idayend == 0 )  THEN  
     765                    
     766                     ! Daily averaged moored buoy (MRB) data  
     767                    
     768                     ! vertically interpolate all 4 corners  
     769                     ista = prodatqc%npvsta(jobs,1)  
     770                     iend = prodatqc%npvend(jobs,1)  
     771                     inum_obs = iend - ista + 1  
     772                     ALLOCATE(interp_corner(2,2,inum_obs),v_indic(inum_obs))  
     773       
     774                     DO iin=1,2  
     775                        DO ijn=1,2  
     776                                        
     777                                        
     778            
     779                           IF ( k1dint == 1 ) THEN  
     780                              CALL obs_int_z1d_spl( kpk, &  
     781                              &     zinmt(iin,ijn,:,jobs), &  
     782                              &     zobs2k, zgdept(iin,ijn,:,jobs), &  
     783                              &     zmask(iin,ijn,:,jobs))  
     784                           ENDIF  
     785        
     786                           CALL obs_level_search(kpk, &  
     787                           &    zgdept(iin,ijn,:,jobs), &  
     788                           &    inum_obs, prodatqc%var(1)%vdep(ista:iend), &  
     789                           &    v_indic)  
     790                           CALL obs_int_z1d(kpk, v_indic, k1dint, inum_obs, &  
     791                           &    prodatqc%var(1)%vdep(ista:iend), &  
     792                           &    zinmt(iin,ijn,:,jobs), &  
     793                           &    zobs2k, interp_corner(iin,ijn,:), &  
     794                           &    zgdept(iin,ijn,:,jobs), &  
     795                           &    zmask(iin,ijn,:,jobs))  
     796        
     797                        ENDDO  
     798                     ENDDO  
     799                    
     800                    
     801                  ELSE  
     802                 
     803                     CALL ctl_stop( ' A nonzero' //     &  
     804                     &           ' number of profile T BUOY data should' // &  
     805                     &           ' only occur at the end of a given day' )  
     806     
     807                  ENDIF  
     808          
     809               ELSE   
     810                 
     811                  ! Point data  
     812      
     813                  ! vertically interpolate all 4 corners  
     814                  ista = prodatqc%npvsta(jobs,1)  
     815                  iend = prodatqc%npvend(jobs,1)  
     816                  inum_obs = iend - ista + 1  
     817                  ALLOCATE(interp_corner(2,2,inum_obs), v_indic(inum_obs))  
     818                  DO iin=1,2   
     819                     DO ijn=1,2  
     820                                     
     821                                     
     822                        IF ( k1dint == 1 ) THEN  
     823                           CALL obs_int_z1d_spl( kpk, &  
     824                           &    zintt(iin,ijn,:,jobs),&  
     825                           &    zobs2k, zgdept(iin,ijn,:,jobs), &  
     826                           &    zmask(iin,ijn,:,jobs))  
     827   
     828                        ENDIF  
     829        
     830                        CALL obs_level_search(kpk, &  
     831                         &        zgdept(iin,ijn,:,jobs),&  
     832                         &        inum_obs, prodatqc%var(1)%vdep(ista:iend), &  
     833                         &         v_indic)  
     834                        CALL obs_int_z1d(kpk, v_indic, k1dint, inum_obs,     &  
     835                         &          prodatqc%var(1)%vdep(ista:iend),     &  
     836                         &          zintt(iin,ijn,:,jobs),            &  
     837                         &          zobs2k,interp_corner(iin,ijn,:), &  
     838                         &          zgdept(iin,ijn,:,jobs),         &  
     839                         &          zmask(iin,ijn,:,jobs) )       
     840          
     841                     ENDDO  
     842                  ENDDO  
     843              
     844               ENDIF  
     845        
     846               !-------------------------------------------------------------  
     847               ! Compute the horizontal interpolation for every profile level  
     848               !-------------------------------------------------------------  
     849              
     850               DO ikn=1,inum_obs  
     851                  iend=ista+ikn-1  
     852    
     853                  ! This code forces the horizontal weights to be   
     854                  ! zero IF the observation is below the bottom of the   
     855                  ! corners of the interpolation nodes, Or if it is in   
     856                  ! the mask. This is important for observations are near   
     857                  ! steep bathymetry  
     858                  DO iin=1,2  
     859                     DO ijn=1,2  
     860      
     861                        depth_loop1: DO ik=kpk,2,-1  
     862                           IF(zmask(iin,ijn,ik-1,jobs ) > 0.9 )THEN    
     863                             
     864                              l_zweig(iin,ijn,1) = &   
     865                           &  zweig(iin,ijn,1) *   &  
     866                           &  MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,jobs) ) &  
     867                           &  - prodatqc%var(1)%vdep(iend)),0._wp)  
     868                             
     869                              EXIT depth_loop1  
     870                           ENDIF  
     871                        ENDDO depth_loop1  
     872      
     873                     ENDDO  
     874                  ENDDO  
     875    
     876                  CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), &  
     877                  &          prodatqc%var(1)%vmod(iend:iend) )  
     878  
     879               ENDDO  
     880  
     881  
     882               DEALLOCATE(interp_corner,v_indic)  
     883           
     884            ENDIF  
     885        
     886  
     887            ! Salinity   
     888           
     889            IF ( prodatqc%npvend(jobs,2) > 0 ) THEN   
     890     
     891               zobsk(:) = obfillflt  
     892     
     893               IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN  
     894     
     895                  IF ( idayend == 0 )  THEN  
     896                    
     897                     ! Daily averaged moored buoy (MRB) data  
     898                    
     899                     ! vertically interpolate all 4 corners  
     900                     ista = prodatqc%npvsta(jobs,2)  
     901                     iend = prodatqc%npvend(jobs,2)  
     902                     inum_obs = iend - ista + 1  
     903                     ALLOCATE(interp_corner(2,2,inum_obs),v_indic(inum_obs))  
     904       
     905                     DO iin=1,2  
     906                        DO ijn=1,2  
     907                                        
     908                                        
     909            
     910                           IF ( k1dint == 1 ) THEN  
     911                              CALL obs_int_z1d_spl( kpk, &  
     912                              &     zinms(iin,ijn,:,jobs), &  
     913                              &     zobs2k, zgdept(iin,ijn,:,jobs), &  
     914                              &     zmask(iin,ijn,:,jobs))  
     915                           ENDIF  
     916        
     917                           CALL obs_level_search(kpk, &  
     918                           &    zgdept(iin,ijn,:,jobs), &  
     919                           &    inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
     920                           &    v_indic)  
     921                           CALL obs_int_z1d(kpk, v_indic, k1dint, inum_obs, &  
     922                           &    prodatqc%var(2)%vdep(ista:iend), &  
     923                           &    zinms(iin,ijn,:,jobs), &  
     924                           &    zobs2k, interp_corner(iin,ijn,:), &  
     925                           &    zgdept(iin,ijn,:,jobs), &  
     926                           &    zmask(iin,ijn,:,jobs))  
     927        
     928                        ENDDO  
     929                     ENDDO  
     930                    
     931                    
     932                  ELSE  
     933                 
     934                     CALL ctl_stop( ' A nonzero' //     &  
     935                     &           ' number of profile T BUOY data should' // &  
     936                     &           ' only occur at the end of a given day' )  
     937     
     938                  ENDIF  
     939          
     940               ELSE   
     941                 
     942                  ! Point data  
     943      
     944                  ! vertically interpolate all 4 corners  
     945                  ista = prodatqc%npvsta(jobs,2)  
     946                  iend = prodatqc%npvend(jobs,2)  
     947                  inum_obs = iend - ista + 1  
     948                  ALLOCATE(interp_corner(2,2,inum_obs), v_indic(inum_obs))  
     949                    
     950                  DO iin=1,2      
     951                     DO ijn=1,2   
     952                                  
     953                                  
     954                        IF ( k1dint == 1 ) THEN  
     955                           CALL obs_int_z1d_spl( kpk, &  
     956                           &    zints(iin,ijn,:,jobs),&  
     957                           &    zobs2k, zgdept(iin,ijn,:,jobs), &  
     958                           &    zmask(iin,ijn,:,jobs))  
     959   
     960                        ENDIF  
     961        
     962                        CALL obs_level_search(kpk, &  
     963                         &        zgdept(iin,ijn,:,jobs),&  
     964                         &        inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
     965                         &         v_indic)  
     966                        CALL obs_int_z1d(kpk, v_indic, k1dint, inum_obs,     &  
     967                         &          prodatqc%var(2)%vdep(ista:iend),     &  
     968                         &          zints(iin,ijn,:,jobs),            &  
     969                         &          zobs2k,interp_corner(iin,ijn,:), &  
     970                         &          zgdept(iin,ijn,:,jobs),         &  
     971                         &          zmask(iin,ijn,:,jobs) )       
     972          
     973                     ENDDO  
     974                  ENDDO  
     975              
     976               ENDIF  
     977        
     978               !-------------------------------------------------------------  
     979               ! Compute the horizontal interpolation for every profile level  
     980               !-------------------------------------------------------------  
     981              
     982               DO ikn=1,inum_obs  
     983                  iend=ista+ikn-1  
     984    
     985                  ! This code forces the horizontal weights to be   
     986                  ! zero IF the observation is below the bottom of the   
     987                  ! corners of the interpolation nodes, Or if it is in   
     988                  ! the mask. This is important for observations are near   
     989                  ! steep bathymetry  
     990                  DO iin=1,2  
     991                     DO ijn=1,2  
     992      
     993                        depth_loop2: DO ik=kpk,2,-1  
     994                           IF(zmask(iin,ijn,ik-1,jobs ) > 0.9 )THEN    
     995                             
     996                              l_zweig(iin,ijn,1) = &   
     997                           &  zweig(iin,ijn,1) *   &  
     998                           &  MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,jobs) ) &  
     999                           &  - prodatqc%var(2)%vdep(iend)),0._wp)  
     1000                             
     1001                              EXIT depth_loop2  
     1002                           ENDIF  
     1003                        ENDDO depth_loop2  
     1004      
     1005                     ENDDO  
     1006                  ENDDO  
     1007    
     1008                  CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), &  
     1009                  &          prodatqc%var(2)%vmod(iend:iend) )  
     1010  
     1011               ENDDO  
     1012  
     1013  
     1014               DEALLOCATE(interp_corner,v_indic)  
     1015           
     1016            ENDIF  
     1017           
     1018         ENDIF  
     1019        
     1020      END DO  
     1021      
     1022      ! Deallocate the data for interpolation  
     1023      DEALLOCATE( &  
     1024         & igrdi, &  
     1025         & igrdj, &  
     1026         & zglam, &  
     1027         & zgphi, &  
     1028         & zmask, &  
     1029         & zintt, &  
     1030         & zints  &  
     1031         & )  
     1032      ! At the end of the day also get interpolated means  
     1033      IF ( idayend == 0 ) THEN  
     1034         DEALLOCATE( &  
     1035            & zinmt,  &  
     1036            & zinms   &  
     1037            & )  
     1038      ENDIF  
     1039     
     1040      prodatqc%nprofup = prodatqc%nprofup + ipro   
     1041        
     1042   END SUBROUTINE obs_pro_sco_opt  
     1043  
    4511044   SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, & 
    4521045      &                    psshn, psshmask, k2dint ) 
     
    12931886            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
    12941887               &                   zglamv(:,:,iobs), zgphiv(:,:,iobs), & 
    1295                &                   zvmask(:,:,:,iobs), zweigv, zobsmasku ) 
     1888               &                   zvmask(:,:,:,iobs), zweigv, zobsmaskv ) 
    12961889 
    12971890         ENDIF 
Note: See TracChangeset for help on using the changeset viewer.