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

Ignore:
Timestamp:
2018-01-05T14:29:29+01:00 (7 years ago)
Author:
dford
Message:

Initial implementation of 3D biogeochemistry observation operator.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

    r8653 r9186  
    231231                  DO ji = 1, jpi 
    232232                     prodatqc%vdmean(ji,jj,jk,1) = 0.0 
    233                      prodatqc%vdmean(ji,jj,jk,2) = 0.0 
     233                     IF ( prodatqc%nvar == 2 ) THEN 
     234                        prodatqc%vdmean(ji,jj,jk,2) = 0.0 
     235                     ENDIF 
    234236                  END DO 
    235237               END DO 
     
    243245                  prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 
    244246                     &                        + pvar1(ji,jj,jk) 
    245                   ! Increment field 2 for computing daily mean 
    246                   prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
    247                      &                        + pvar2(ji,jj,jk) 
     247                  IF ( prodatqc%nvar == 2 ) THEN 
     248                     ! Increment field 2 for computing daily mean 
     249                     prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
     250                        &                        + pvar2(ji,jj,jk) 
     251                  ENDIF 
    248252               END DO 
    249253            END DO 
     
    260264                     prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 
    261265                        &                        * zdaystp 
    262                      prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
    263                         &                        * zdaystp 
     266                     IF ( prodatqc%nvar == 2 ) THEN 
     267                        prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
     268                           &                        * zdaystp 
     269                     ENDIF 
    264270                  END DO 
    265271               END DO 
     
    297303         igrdi1(2,2,iobs) = prodatqc%mi(jobs,1) 
    298304         igrdj1(2,2,iobs) = prodatqc%mj(jobs,1) 
    299          igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1 
    300          igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1 
    301          igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1 
    302          igrdj2(1,2,iobs) = prodatqc%mj(jobs,2) 
    303          igrdi2(2,1,iobs) = prodatqc%mi(jobs,2) 
    304          igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1 
    305          igrdi2(2,2,iobs) = prodatqc%mi(jobs,2) 
    306          igrdj2(2,2,iobs) = prodatqc%mj(jobs,2) 
     305         IF ( prodatqc%nvar == 2 ) THEN 
     306            igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1 
     307            igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1 
     308            igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1 
     309            igrdj2(1,2,iobs) = prodatqc%mj(jobs,2) 
     310            igrdi2(2,1,iobs) = prodatqc%mi(jobs,2) 
     311            igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1 
     312            igrdi2(2,2,iobs) = prodatqc%mi(jobs,2) 
     313            igrdj2(2,2,iobs) = prodatqc%mj(jobs,2) 
     314         ENDIF 
    307315      END DO 
    308316 
     
    316324      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pvar1,   zint1 ) 
    317325       
    318       CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, plam2, zglam2 ) 
    319       CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, pphi2, zgphi2 ) 
    320       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pmask2, zmask2 ) 
    321       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2,   zint2 ) 
     326      IF ( prodatqc%nvar == 2 ) THEN 
     327         CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, plam2, zglam2 ) 
     328         CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, pphi2, zgphi2 ) 
     329         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pmask2, zmask2 ) 
     330         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2,   zint2 ) 
     331      ENDIF 
    322332 
    323333      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept )  
     
    334344         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, & 
    335345            &                  prodatqc%vdmean(:,:,:,1), zinm1 ) 
    336          CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, & 
    337             &                  prodatqc%vdmean(:,:,:,2), zinm2 ) 
     346         IF ( prodatqc%nvar == 2 ) THEN 
     347            CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, & 
     348               &                  prodatqc%vdmean(:,:,:,2), zinm2 ) 
     349         ENDIF 
    338350 
    339351      ENDIF 
     
    378390         ENDIF 
    379391 
    380          IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    381  
    382             CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     & 
    383                &                   zglam2(:,:,iobs), zgphi2(:,:,iobs), & 
    384                &                   zmask2(:,:,1,iobs), zweig2, zmsk_2 ) 
    385   
     392         IF ( prodatqc%nvar == 2 ) THEN 
     393            IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
     394 
     395               CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     & 
     396                  &                   zglam2(:,:,iobs), zgphi2(:,:,iobs), & 
     397                  &                   zmask2(:,:,1,iobs), zweig2, zmsk_2 ) 
     398 
     399            ENDIF 
    386400         ENDIF 
    387401 
     
    512526         ENDIF  
    513527 
    514          ! For the second variable 
    515          IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    516  
    517             zobsk(:) = obfillflt 
    518  
    519             IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 
    520  
    521                IF ( idayend == 0 )  THEN 
    522                   ! Daily averaged data 
     528         IF ( prodatqc%nvar == 2 ) THEN 
     529            ! For the second variable 
     530            IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
     531 
     532               zobsk(:) = obfillflt 
     533 
     534               IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 
     535 
     536                  IF ( idayend == 0 )  THEN 
     537                     ! Daily averaged data 
     538 
     539                     ! vertically interpolate all 4 corners  
     540                     ista = prodatqc%npvsta(jobs,2)  
     541                     iend = prodatqc%npvend(jobs,2)  
     542                     inum_obs = iend - ista + 1  
     543                     ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
     544 
     545                     DO iin=1,2  
     546                        DO ijn=1,2  
     547 
     548                           IF ( k1dint == 1 ) THEN  
     549                              CALL obs_int_z1d_spl( kpk, &  
     550                                 &     zinm2(iin,ijn,:,iobs), &  
     551                                 &     zobs2k, zgdept(iin,ijn,:,iobs), &  
     552                                 &     zmask2(iin,ijn,:,iobs))  
     553                           ENDIF  
     554 
     555                           CALL obs_level_search(kpk, &  
     556                              &    zgdept(iin,ijn,:,iobs), &  
     557                              &    inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
     558                              &    iv_indic)  
     559 
     560                           CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
     561                              &    prodatqc%var(2)%vdep(ista:iend), &  
     562                              &    zinm2(iin,ijn,:,iobs), &  
     563                              &    zobs2k, interp_corner(iin,ijn,:), &  
     564                              &    zgdept(iin,ijn,:,iobs), &  
     565                              &    zmask2(iin,ijn,:,iobs))  
     566 
     567                        ENDDO  
     568                     ENDDO  
     569 
     570                  ENDIF !idayend 
     571 
     572               ELSE    
     573 
     574                  ! Point data  
    523575 
    524576                  ! vertically interpolate all 4 corners  
     
    526578                  iend = prodatqc%npvend(jobs,2)  
    527579                  inum_obs = iend - ista + 1  
    528                   ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
    529  
     580                  ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
     581                  DO iin=1,2   
     582                     DO ijn=1,2  
     583 
     584                        IF ( k1dint == 1 ) THEN  
     585                           CALL obs_int_z1d_spl( kpk, &  
     586                              &    zint2(iin,ijn,:,iobs),&  
     587                              &    zobs2k, zgdept(iin,ijn,:,iobs), &  
     588                              &    zmask2(iin,ijn,:,iobs))  
     589 
     590                        ENDIF  
     591 
     592                        CALL obs_level_search(kpk, &  
     593                            &        zgdept(iin,ijn,:,iobs),&  
     594                            &        inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
     595                            &        iv_indic)  
     596 
     597                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
     598                            &          prodatqc%var(2)%vdep(ista:iend),     &  
     599                            &          zint2(iin,ijn,:,iobs),            &  
     600                            &          zobs2k,interp_corner(iin,ijn,:), &  
     601                            &          zgdept(iin,ijn,:,iobs),         &  
     602                            &          zmask2(iin,ijn,:,iobs) )       
     603 
     604                     ENDDO  
     605                  ENDDO  
     606 
     607               ENDIF  
     608 
     609               !-------------------------------------------------------------  
     610               ! Compute the horizontal interpolation for every profile level  
     611               !-------------------------------------------------------------  
     612 
     613               DO ikn=1,inum_obs  
     614                  iend=ista+ikn-1 
     615 
     616                  zweig(:,:,1) = 0._wp  
     617 
     618                  ! This code forces the horizontal weights to be   
     619                  ! zero IF the observation is below the bottom of the   
     620                  ! corners of the interpolation nodes, Or if it is in   
     621                  ! the mask. This is important for observations near   
     622                  ! steep bathymetry  
    530623                  DO iin=1,2  
    531624                     DO ijn=1,2  
    532625 
    533                         IF ( k1dint == 1 ) THEN  
    534                            CALL obs_int_z1d_spl( kpk, &  
    535                               &     zinm2(iin,ijn,:,iobs), &  
    536                               &     zobs2k, zgdept(iin,ijn,:,iobs), &  
    537                               &     zmask2(iin,ijn,:,iobs))  
    538                         ENDIF  
    539         
    540                         CALL obs_level_search(kpk, &  
    541                            &    zgdept(iin,ijn,:,iobs), &  
    542                            &    inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
    543                            &    iv_indic)  
    544  
    545                         CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
    546                            &    prodatqc%var(2)%vdep(ista:iend), &  
    547                            &    zinm2(iin,ijn,:,iobs), &  
    548                            &    zobs2k, interp_corner(iin,ijn,:), &  
    549                            &    zgdept(iin,ijn,:,iobs), &  
    550                            &    zmask2(iin,ijn,:,iobs))  
    551         
     626                        depth_loop2: DO ik=kpk,2,-1  
     627                           IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
     628 
     629                              zweig(iin,ijn,1) = &   
     630                                 & zweig2(iin,ijn,1) * &  
     631                                 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
     632                                 &  - prodatqc%var(2)%vdep(iend)),0._wp)  
     633 
     634                              EXIT depth_loop2  
     635 
     636                           ENDIF  
     637 
     638                        ENDDO depth_loop2  
     639 
    552640                     ENDDO  
    553641                  ENDDO  
    554642 
    555                ENDIF !idayend 
    556  
    557             ELSE    
    558  
    559                ! Point data  
    560       
    561                ! vertically interpolate all 4 corners  
    562                ista = prodatqc%npvsta(jobs,2)  
    563                iend = prodatqc%npvend(jobs,2)  
    564                inum_obs = iend - ista + 1  
    565                ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
    566                DO iin=1,2   
    567                   DO ijn=1,2  
    568                      
    569                      IF ( k1dint == 1 ) THEN  
    570                         CALL obs_int_z1d_spl( kpk, &  
    571                            &    zint2(iin,ijn,:,iobs),&  
    572                            &    zobs2k, zgdept(iin,ijn,:,iobs), &  
    573                            &    zmask2(iin,ijn,:,iobs))  
    574    
    575                      ENDIF  
    576         
    577                      CALL obs_level_search(kpk, &  
    578                          &        zgdept(iin,ijn,:,iobs),&  
    579                          &        inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
    580                          &        iv_indic)  
    581  
    582                      CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
    583                          &          prodatqc%var(2)%vdep(ista:iend),     &  
    584                          &          zint2(iin,ijn,:,iobs),            &  
    585                          &          zobs2k,interp_corner(iin,ijn,:), &  
    586                          &          zgdept(iin,ijn,:,iobs),         &  
    587                          &          zmask2(iin,ijn,:,iobs) )       
    588           
    589                   ENDDO  
     643                  CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &  
     644                     &              prodatqc%var(2)%vmod(iend:iend) )  
     645 
     646                     ! Set QC flag for any observations found below the bottom 
     647                     ! needed as the check here is more strict than that in obs_prep 
     648                  IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4 
     649 
    590650               ENDDO  
    591               
     651 
     652               DEALLOCATE(interp_corner,iv_indic)  
     653 
    592654            ENDIF  
    593  
    594             !-------------------------------------------------------------  
    595             ! Compute the horizontal interpolation for every profile level  
    596             !-------------------------------------------------------------  
    597               
    598             DO ikn=1,inum_obs  
    599                iend=ista+ikn-1 
    600                    
    601                zweig(:,:,1) = 0._wp  
    602     
    603                ! This code forces the horizontal weights to be   
    604                ! zero IF the observation is below the bottom of the   
    605                ! corners of the interpolation nodes, Or if it is in   
    606                ! the mask. This is important for observations near   
    607                ! steep bathymetry  
    608                DO iin=1,2  
    609                   DO ijn=1,2  
    610       
    611                      depth_loop2: DO ik=kpk,2,-1  
    612                         IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
    613                              
    614                            zweig(iin,ijn,1) = &   
    615                               & zweig2(iin,ijn,1) * &  
    616                               & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
    617                               &  - prodatqc%var(2)%vdep(iend)),0._wp)  
    618                              
    619                            EXIT depth_loop2  
    620  
    621                         ENDIF  
    622  
    623                      ENDDO depth_loop2  
    624       
    625                   ENDDO  
    626                ENDDO  
    627     
    628                CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &  
    629                   &              prodatqc%var(2)%vmod(iend:iend) )  
    630  
    631                   ! Set QC flag for any observations found below the bottom 
    632                   ! needed as the check here is more strict than that in obs_prep 
    633                IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4 
    634   
    635             ENDDO  
    636   
    637             DEALLOCATE(interp_corner,iv_indic)  
    638            
    639655         ENDIF  
    640656 
Note: See TracChangeset for help on using the changeset viewer.