Changeset 11449


Ignore:
Timestamp:
2019-08-19T10:54:47+02:00 (22 months ago)
Author:
mattmartin
Message:

Committed first version of changes to output climatology values at obs locations in the fdbk files.

Location:
branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90

    r11235 r11449  
    3333   USE mpp_map                  ! MPP mapping 
    3434   USE lib_mpp                  ! For ctl_warn/stop 
     35   USE dtatsd                   ! For climatological temperature & salinity 
     36   USE tradmp                   ! "" 
    3537 
    3638   IMPLICIT NONE 
     
    4850   LOGICAL :: ln_diaobs            !: Logical switch for the obs operator 
    4951   LOGICAL :: ln_sstnight          !: Logical switch for night mean SST obs 
     52   LOGICAL :: ln_output_clim       !: Logical switch for interpolating and outputting T/S climatology 
    5053   LOGICAL :: ln_default_fp_indegs !: T=> Default obs footprint size specified in degrees, F=> in metres 
    5154   LOGICAL :: ln_sla_fp_indegs     !: T=>     SLA obs footprint size specified in degrees, F=> in metres 
     
    254257         &            ln_grid_global, ln_grid_search_lookup,          & 
    255258         &            ln_ignmis, ln_s_at_t, ln_bound_reject,          & 
    256          &            ln_sstnight, ln_default_fp_indegs,              & 
     259         &            ln_sstnight,  ln_output_clim,                   & 
     260         &            ln_default_fp_indegs,                           & 
    257261         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             & 
    258262         &            ln_sss_fp_indegs, ln_sic_fp_indegs,             & 
     
    415419         WRITE(numout,*) '             Daily average types                             nn_profdavtypes = ', nn_profdavtypes 
    416420         WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight 
     421         WRITE(numout,*) '             Logical switch for writing climat. at ob locs    ln_output_clim = ', ln_output_clim 
    417422      ENDIF 
    418423      !----------------------------------------------------------------------- 
     
    438443         RETURN 
    439444      ENDIF 
     445 
     446      IF ( ln_output_clim .AND. ( .NOT. ln_tradmp ) ) THEN 
     447         IF(lwp) WRITE(numout,cform_war) 
     448         IF(lwp) WRITE(numout,*) ' ln_output_clim is true, but ln_tradmp is false', & 
     449            &                    ' so climatological T/S not available and will not be output' 
     450         nwarn = nwarn + 1 
     451         ln_output_clim = .FALSE. 
     452      ENDIF 
     453      
    440454 
    441455      IF(lwp) WRITE(numout,*) '          Number of profile obs types: ',nproftypes 
     
    926940      REAL(wp) :: tiny             ! small number 
    927941      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 
    928          & zprofvar                ! Model values for variables in a prof ob 
     942         & zprofvar, &             ! Model values for variables in a prof ob 
     943         & zprofclim               ! Climatology values for variables in a prof ob 
    929944      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 
    930945         & zprofmask               ! Mask associated with zprofvar 
    931946      REAL(wp), POINTER, DIMENSION(:,:) :: & 
    932947         & zsurfvar, &             ! Model values equivalent to surface ob. 
     948         & zsurfclim, &            ! Climatology values for variables in a surface ob. 
    933949         & zsurfmask               ! Mask associated with surface variable 
    934950      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
     
    940956         & pco2_3d                 ! 3D pCO2 from FABM 
    941957#endif 
    942  
     958      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zts_dta  
     959       
    943960      IF(lwp) THEN 
    944961         WRITE(numout,*) 
     
    950967      idaystp = NINT( rday / rdt ) 
    951968 
     969      ! Get the climatological T & S fields on this time step 
     970      IF ( ln_output_clim ) CALL dta_tsd( kstp, zts_dta ) 
     971 
    952972      !----------------------------------------------------------------------- 
    953973      ! Call the profile and surface observation operators 
     
    963983            CALL wrk_alloc( jpi, jpj,      profdataqc(jtype)%nvar, zglam     ) 
    964984            CALL wrk_alloc( jpi, jpj,      profdataqc(jtype)%nvar, zgphi     ) 
    965              
     985            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofclim )     
     986                               
    966987            ! Defaults which might change 
    967988            DO jvar = 1, profdataqc(jtype)%nvar 
     
    969990               zglam(:,:,jvar)       = glamt(:,:) 
    970991               zgphi(:,:,jvar)       = gphit(:,:) 
     992               zprofclim(:,:,:,jvar) = 0._wp 
    971993            END DO 
    972994 
     
    976998               zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem) 
    977999               zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal) 
    978  
     1000               IF ( ln_output_clim ) THEN                
     1001                  zprofclim(:,:,:,1) = zts_dta(:,:,:,jp_tem) 
     1002                  zprofclim(:,:,:,2) = zts_dta(:,:,:,jp_sal) 
     1003               ENDIF 
     1004                
    9791005            CASE('vel') 
    9801006               zprofvar(:,:,:,1) = un(:,:,:) 
     
    11551181                  &               nit000, idaystp, jvar,                   & 
    11561182                  &               zprofvar(:,:,:,jvar),                    & 
     1183                  &               zprofclim(:,:,:,jvar),                   & 
    11571184                  &               fsdept(:,:,:), fsdepw(:,:,:),            &  
    11581185                  &               zprofmask(:,:,:,jvar),                   & 
     
    11661193            CALL wrk_dealloc( jpi, jpj,      profdataqc(jtype)%nvar, zglam     ) 
    11671194            CALL wrk_dealloc( jpi, jpj,      profdataqc(jtype)%nvar, zgphi     ) 
     1195            CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofclim  )             
    11681196 
    11691197         END DO 
     
    11751203         !Allocate local work arrays 
    11761204         CALL wrk_alloc( jpi, jpj, zsurfvar ) 
     1205         CALL wrk_alloc( jpi, jpj, zsurfclim )          
    11771206         CALL wrk_alloc( jpi, jpj, zsurfmask ) 
    11781207#if defined key_fabm 
     
    11841213            !Defaults which might be changed 
    11851214            zsurfmask(:,:) = tmask(:,:,1) 
     1215            zsurfclim(:,:) = 0._wp           
    11861216            llog10 = .FALSE. 
    11871217 
     
    11891219            CASE('sst') 
    11901220               zsurfvar(:,:) = tsn(:,:,1,jp_tem) 
     1221               IF ( ln_output_clim ) zsurfclim(:,:) = zts_dta(:,:,1,jp_tem) 
    11911222            CASE('sla') 
    1192                zsurfvar(:,:) = sshn(:,:) 
     1223               zsurfvar(:,:)  = sshn(:,:) 
    11931224            CASE('sss') 
    11941225               zsurfvar(:,:) = tsn(:,:,1,jp_sal) 
     1226               IF ( ln_output_clim ) zsurfclim(:,:) = zts_dta(:,:,1,jp_sal)                
    11951227            CASE('sic') 
    11961228               IF ( kstp == 0 ) THEN 
     
    14851517 
    14861518            CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
    1487                &               nit000, idaystp, zsurfvar, zsurfmask,    & 
     1519               &               nit000, idaystp, zsurfvar,               & 
     1520               &               zsurfclim, zsurfmask,                    & 
    14881521               &               n2dintsurf(jtype), llnightav(jtype),     & 
    1489                &               ravglamscl(jtype), ravgphiscl(jtype),     & 
     1522               &               ravglamscl(jtype), ravgphiscl(jtype),    & 
    14901523               &               lfpindegs(jtype) ) 
    14911524 
  • branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

    r9306 r11449  
    6262   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, & 
    6363      &                     kit000, kdaystp, kvar,       & 
    64       &                     pvar, pgdept, pgdepw,        & 
    65       &                     pmask,                       &   
     64      &                     pvar, pclim,                 & 
     65      &                     pgdept, pgdepw, pmask,       &   
    6666      &                     plam, pphi,                  & 
    6767      &                     k1dint, k2dint, kdailyavtypes ) 
     
    137137      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 
    138138         & pvar,   &                 ! Model field for variable 
     139         & pclim,  &                 ! Climatology field for variable          
    139140         & pmask                     ! Land-sea mask for variable 
    140141      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
     
    172173      REAL(KIND=wp), DIMENSION(kpk) :: & 
    173174         & zobsk,    & 
    174          & zobs2k 
     175         & zobs2k,   & 
     176         & zclm2k          
    175177      REAL(KIND=wp), DIMENSION(2,2,1) :: & 
    176178         & zweig1, & 
     
    178180      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 
    179181         & zmask,  & 
     182         & zclim,  &          
    180183         & zint,   & 
    181184         & zinm,   & 
     
    189192 
    190193      LOGICAL :: ld_dailyav 
     194      LOGICAL :: ld_clim 
    191195 
    192196      !------------------------------------------------------------------------ 
     
    196200      inrc = kt - kit000 + 2 
    197201      ipro = prodatqc%npstp(inrc) 
     202 
     203      ! Check if climatology is available and set flag 
     204      IF ( SUM( pclim(:,:,:) ) == 0. ) THEN 
     205         ld_clim = .FALSE. 
     206      ELSE 
     207         ld_clim = .TRUE.       
     208      ENDIF 
    198209 
    199210      ! Daily average types 
     
    262273         & ) 
    263274 
     275      IF ( ld_clim ) ALLOCATE( zclim(2,2,kpk,ipro) ) 
     276 
    264277      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
    265278         iobs = jobs - prodatqc%nprofup 
     
    286299      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw, zgdepw )  
    287300 
     301      IF ( ld_clim ) THEN 
     302         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pclim, zclim )             
     303      ENDIF  
     304       
    288305      ! At the end of the day also get interpolated means 
    289306      IF ( ld_dailyav .AND. idayend == 0 ) THEN 
     
    349366                  inum_obs = iend - ista + 1  
    350367                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
    351  
     368                  IF ( ld_clim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) ) 
     369                   
    352370                  DO iin=1,2  
    353371                     DO ijn=1,2  
     
    358376                              &     zobs2k, zgdept(iin,ijn,:,iobs), &  
    359377                              &     zmask(iin,ijn,:,iobs))  
     378 
     379                           IF ( ld_clim ) THEN 
     380                              CALL obs_int_z1d_spl( kpk, &  
     381                                 &     zclim(iin,ijn,:,iobs), &  
     382                                 &     zclm2k, zgdept(iin,ijn,:,iobs), &  
     383                                 &     zmask(iin,ijn,:,iobs))  
     384                           ENDIF 
     385 
    360386                        ENDIF  
    361387        
     
    371397                           &    zgdept(iin,ijn,:,iobs), &  
    372398                           &    zmask(iin,ijn,:,iobs))  
    373         
     399 
     400                        IF ( ld_clim ) THEN 
     401                           CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
     402                              &    prodatqc%var(kvar)%vdep(ista:iend), &  
     403                              &    zclim(iin,ijn,:,iobs), &  
     404                              &    zclm2k, interp_corner_clim(iin,ijn,:), &  
     405                              &    zgdept(iin,ijn,:,iobs), &  
     406                              &    zmask(iin,ijn,:,iobs))  
     407                        ENDIF 
     408                         
    374409                     ENDDO  
    375410                  ENDDO  
     
    385420               iend = prodatqc%npvend(jobs,kvar)  
    386421               inum_obs = iend - ista + 1  
    387                ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
     422               ALLOCATE( interp_corner(2,2,inum_obs),      & 
     423                  &      iv_indic(inum_obs) )  
     424               IF ( ld_clim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) )                   
    388425               DO iin=1,2   
    389426                  DO ijn=1,2  
     
    394431                           &    zobs2k, zgdept(iin,ijn,:,iobs), &  
    395432                           &    zmask(iin,ijn,:,iobs))  
     433 
     434                        IF ( ld_clim ) THEN 
     435                           CALL obs_int_z1d_spl( kpk, &  
     436                              &    zclim(iin,ijn,:,iobs),&  
     437                              &    zclm2k, zgdept(iin,ijn,:,iobs), &  
     438                              &    zmask(iin,ijn,:,iobs))  
     439                        ENDIF 
    396440   
    397441                     ENDIF  
     
    408452                         &          zgdept(iin,ijn,:,iobs),         &  
    409453                         &          zmask(iin,ijn,:,iobs) )       
     454 
     455                     IF ( ld_clim ) THEN 
     456                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
     457                            &          prodatqc%var(kvar)%vdep(ista:iend),     &  
     458                            &          zclim(iin,ijn,:,iobs),            &  
     459                            &          zclm2k,interp_corner_clim(iin,ijn,:), &  
     460                            &          zgdept(iin,ijn,:,iobs),         &  
     461                            &          zmask(iin,ijn,:,iobs) )    
     462                     ENDIF    
    410463          
    411464                  ENDDO  
     
    451504                  &              prodatqc%var(kvar)%vmod(iend:iend) )  
    452505 
     506               IF ( ld_clim ) THEN 
     507                  CALL obs_int_h2d( 1, 1, zweig, interp_corner_clim(:,:,ikn), &  
     508                     &              prodatqc%var(kvar)%vclm(iend:iend) ) 
     509               ENDIF 
     510 
    453511                  ! Set QC flag for any observations found below the bottom 
    454512                  ! needed as the check here is more strict than that in obs_prep 
     
    458516  
    459517            DEALLOCATE(interp_corner,iv_indic)  
    460            
     518            IF ( ld_clim ) DEALLOCATE( interp_corner_clim )          
     519              
    461520         ENDIF 
    462521 
     
    475534         & ) 
    476535 
     536      IF ( ld_clim ) DEALLOCATE( zclim ) 
     537       
    477538      ! At the end of the day also get interpolated means 
    478539      IF ( ld_dailyav .AND. idayend == 0 ) THEN 
     
    487548 
    488549   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,            & 
    489       &                     kit000, kdaystp, psurf, psurfmask,   & 
     550      &                     kit000, kdaystp, psurf, pclim, psurfmask,   & 
    490551      &                     k2dint, ldnightav, plamscl, pphiscl, & 
    491552      &                     lindegrees ) 
     
    541602      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
    542603         & psurf,  &                   ! Model surface field 
     604         & pclim,  &                   ! Climatological surface field          
    543605         & psurfmask                   ! Land-sea mask 
    544606      LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data 
     
    569631      REAL(wp) :: zlam 
    570632      REAL(wp) :: zphi 
    571       REAL(wp), DIMENSION(1) :: zext, zobsmask 
     633      REAL(wp), DIMENSION(1) :: zext, zobsmask, zclm 
    572634      REAL(wp) :: zdaystp 
    573635      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
     
    577639         & zsurfm, & 
    578640         & zsurftmp, & 
     641         & zclim,  & 
    579642         & zglam,  & 
    580643         & zgphi,  & 
     
    586649         & zouttmp, & 
    587650         & zmeanday    ! to compute model sst in region of 24h daylight (pole) 
    588  
     651          
     652      LOGICAL :: ld_clim ! T => climatological data is available 
    589653      !------------------------------------------------------------------------ 
    590654      ! Local initialization  
     
    593657      inrc = kt - kit000 + 2 
    594658      isurf = surfdataqc%nsstp(inrc) 
     659       
     660      ! Check if climatological information is available 
     661      IF ( SUM(pclim(:,:)) == 0._wp ) THEN 
     662        ld_clim = .FALSE. 
     663      ELSE 
     664        ld_clim = .TRUE.       
     665      ENDIF 
    595666 
    596667      ! Work out the maximum footprint size for the  
     
    679750         & ) 
    680751 
     752      IF ( ld_clim ) ALLOCATE( zclim(imaxifp,imaxjfp,isurf) ) 
     753 
    681754      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 
    682755         iobs = jobs - surfdataqc%nsurfup 
     
    720793         &                  igrdip1, igrdjp1, gphif, zgphif ) 
    721794 
     795      IF ( ld_clim ) THEN  
     796         CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
     797            &                  igrdi, igrdj, pclim, zclim ) 
     798      ENDIF 
     799 
    722800      ! At the end of the day get interpolated means 
    723801      IF ( idayend == 0 .AND. ldnightav ) THEN 
     
    775853            CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 
    776854 
     855            IF ( ld_clim ) THEN   
     856               CALL obs_int_h2d( 1, 1, zweig, zclim(:,:,iobs), zclm ) 
     857            ENDIF 
     858 
     859 
    777860         ELSE 
    778861 
     
    788871               &              zweig, zsurftmp(:,:,iobs),  zext ) 
    789872 
     873            IF ( ld_clim ) THEN   
     874               CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 
     875                  &              zweig, zsurftmp(:,:,iobs),  zclm ) 
     876            ENDIF 
     877 
    790878         ENDIF 
    791879 
     
    797885            surfdataqc%rmod(jobs,1) = zext(1) 
    798886         ENDIF 
     887          
     888         IF ( ld_clim ) surfdataqc%rclm(jobs,1) = zclm(1) 
    799889          
    800890         IF ( zext(1) == obfillflt ) THEN 
     
    821911         & ) 
    822912 
     913      IF ( ld_clim ) DEALLOCATE( zclim ) 
     914 
    823915      ! At the end of the day also deallocate night-time mean array 
    824916      IF ( idayend == 0 .AND. ldnightav ) THEN 
  • branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_profiles_def.F90

    r7992 r11449  
    7272         & vdep,  &       !: Depth coordinate of profile data 
    7373         & vobs,  &       !: Profile data 
    74          & vmod           !: Model counterpart of the profile data vector 
    75  
     74         & vmod,  &       !: Model counterpart of the profile data vector 
     75         & vclm,  &       !: Climatological counterpart of the profile data vector 
     76          
    7677      REAL(KIND=wp), POINTER, DIMENSION(:,:) :: & 
    7778         & vext           !: Extra variables 
     
    492493         & prof%var(kvar)%vobs(kobs),      & 
    493494         & prof%var(kvar)%vmod(kobs),      & 
     495         & prof%var(kvar)%vclm(kobs),      &          
    494496         & prof%var(kvar)%nvind(kobs)      & 
    495497         & ) 
     
    530532         & prof%var(kvar)%vobs,   & 
    531533         & prof%var(kvar)%vmod,   & 
     534         & prof%var(kvar)%vclm,   &          
    532535         & prof%var(kvar)%nvind,  & 
    533536         & prof%var(kvar)%idqcf,  & 
     
    741744                     newprof%var(jvar)%vmod(invpro(jvar))   = & 
    742745                        &                           prof%var(jvar)%vmod(jj) 
     746                     newprof%var(jvar)%vclm(invpro(jvar))   = & 
     747                        &                           prof%var(jvar)%vclm(jj) 
    743748                     DO jext = 1, prof%next 
    744749                        newprof%var(jvar)%vext(invpro(jvar),jext) = & 
     
    864869               oldprof%var(jvar)%vobs(jl)   = prof%var(jvar)%vobs(jj) 
    865870               oldprof%var(jvar)%vmod(jl)   = prof%var(jvar)%vmod(jj) 
     871               oldprof%var(jvar)%vclm(jl)   = prof%var(jvar)%vclm(jj)                
    866872               oldprof%var(jvar)%idqcf(:,jl) = prof%var(jvar)%idqcf(:,jj) 
    867873               oldprof%var(jvar)%nvqcf(:,jl) = prof%var(jvar)%nvqcf(:,jj) 
  • branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_readmdt.F90

    r7992 r11449  
    3737   REAL(wp), PUBLIC :: rn_mdtcorr   = 1.61_wp  ! User specified MDT correction 
    3838   REAL(wp), PUBLIC :: rn_mdtcutoff = 65.0_wp  ! MDT cutoff for computed correction 
    39  
     39   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  mdt_fld    !: Mean Dynamic Topography on the model grid 
     40    
    4041   !!---------------------------------------------------------------------- 
    4142   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     
    9293      CALL iom_get ( nummdt, jpdom_data, 'sossheig', z_mdt(:,:), 1 ) 
    9394      CALL iom_close(nummdt)                 ! Close the file 
     95       
     96      ! Set the MDT gridded field to be the one read in directly with no adjustments 
     97      ALLOCATE( mdt_fld(jpi,jpj) ) 
     98      mdt_fld(:,:) = z_mdt(:,:) 
    9499       
    95100      ! Read in the fill value 
  • branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90

    r9308 r11449  
    8484      REAL(KIND=wp), POINTER, DIMENSION(:,:) :: & 
    8585         & robs, &        !: Surface observation  
    86          & rmod           !: Model counterpart of the surface observation vector 
    87  
     86         & rmod, &        !: Model counterpart of the surface observation vector 
     87         & rclm           !: Climatological counterpart of the surface observation vector 
     88          
    8889      REAL(KIND=wp), POINTER, DIMENSION(:,:) :: & 
    8990         & rext           !: Extra fields interpolated to observation points 
     
    197198      ALLOCATE( &  
    198199         & surf%robs(ksurf,kvar), & 
    199          & surf%rmod(ksurf,kvar)  & 
     200         & surf%rmod(ksurf,kvar), & 
     201         & surf%rclm(ksurf,kvar)  &          
    200202         & )    
    201203 
     
    290292      DEALLOCATE( &  
    291293         & surf%robs,    & 
    292          & surf%rmod     & 
     294         & surf%rmod,    & 
     295         & surf%rclm     &          
    293296         & ) 
    294297 
     
    418421               newsurf%robs(insurf,jk)  = surf%robs(ji,jk) 
    419422               newsurf%rmod(insurf,jk)  = surf%rmod(ji,jk) 
     423               newsurf%rclm(insurf,jk)  = surf%rclm(ji,jk)                
    420424                
    421425            END DO 
     
    514518            oldsurf%robs(jj,jk)  = surf%robs(ji,jk) 
    515519            oldsurf%rmod(jj,jk)  = surf%rmod(ji,jk) 
     520            oldsurf%rclm(jj,jk)  = surf%rclm(ji,jk)             
    516521 
    517522         END DO 
  • branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90

    r9308 r11449  
    2929   USE obs_mpp              ! MPP support routines for observation diagnostics 
    3030   USE lib_mpp        ! MPP routines 
     31   USE diaobs, ONLY: & 
     32      & ln_output_clim 
    3133 
    3234   IMPLICIT NONE 
     
    9597      INTEGER :: je 
    9698      INTEGER :: iadd 
     99      INTEGER :: iadd_exp  ! expected additional variables 
    97100      INTEGER :: iext 
    98101      REAL(wp) :: zpres 
    99102 
     103 
     104      ! Set up number of additional variables to be ouput: 
     105      ! Hx, CLIM, ... 
     106      iadd_exp = 1   ! Hx 
     107      IF ( ln_output_clim ) iadd_exp = iadd_exp + 1 
     108       
    100109      IF ( PRESENT( padd ) ) THEN 
    101110         iadd = padd%inum 
     
    123132         clfiletype='profb' 
    124133         CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, & 
    125             &                 1 + iadd, 1 + iext, .TRUE. ) 
     134            &                 iadd_exp + iadd, 1 + iext, .TRUE. ) 
    126135         fbdata%cname(1)      = profdata%cvars(1) 
    127136         fbdata%cname(2)      = profdata%cvars(2) 
     
    130139         fbdata%cobunit(1)    = 'Degrees centigrade' 
    131140         fbdata%cobunit(2)    = 'PSU' 
     141          
    132142         fbdata%cextname(1)   = 'TEMP' 
    133143         fbdata%cextlong(1)   = 'Insitu temperature' 
    134144         fbdata%cextunit(1)   = 'Degrees centigrade' 
     145          
    135146         fbdata%caddlong(1,1) = 'Model interpolated potential temperature' 
    136147         fbdata%caddlong(1,2) = 'Model interpolated practical salinity' 
    137148         fbdata%caddunit(1,1) = 'Degrees centigrade' 
    138149         fbdata%caddunit(1,2) = 'PSU' 
     150         IF ( ln_output_clim ) THEN 
     151            fbdata%caddlong(2,1) = 'Climatology interpolated potential temperature' 
     152            fbdata%caddlong(2,2) = 'Climatology interpolated practical salinity' 
     153            fbdata%caddunit(2,1) = 'Degrees centigrade' 
     154            fbdata%caddunit(2,2) = 'PSU' 
     155         ENDIF 
    139156         fbdata%cgrid(:)      = 'T' 
    140157         DO je = 1, iext 
     
    144161         END DO 
    145162         DO ja = 1, iadd 
    146             fbdata%caddname(1+ja) = padd%cdname(ja) 
     163            fbdata%caddname(iadd_exp+ja) = padd%cdname(ja) 
    147164            DO jvar = 1, 2 
    148                fbdata%caddlong(1+ja,jvar) = padd%cdlong(ja,jvar) 
    149                fbdata%caddunit(1+ja,jvar) = padd%cdunit(ja,jvar) 
     165               fbdata%caddlong(iadd_exp+ja,jvar) = padd%cdlong(ja,jvar) 
     166               fbdata%caddunit(iadd_exp+ja,jvar) = padd%cdunit(ja,jvar) 
    150167            END DO 
    151168         END DO 
     
    154171 
    155172         clfiletype='velfb' 
    156          CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, 1, 0, .TRUE. ) 
     173         CALL alloc_obfbdata( fbdata, 2, profdata%nprof, ilevel, & 
     174            &                 iadd_exp + iadd, 0, .TRUE. ) 
    157175         fbdata%cname(1)      = profdata%cvars(1) 
    158176         fbdata%cname(2)      = profdata%cvars(2) 
     
    170188         fbdata%caddunit(1,1) = 'm/s' 
    171189         fbdata%caddunit(1,2) = 'm/s' 
     190         IF ( ln_output_clim ) THEN 
     191            fbdata%caddlong(2,1) = 'Climatology interpolated zonal velocity' 
     192            fbdata%caddlong(2,2) = 'Climatology interpolated meridional velocity' 
     193            fbdata%caddunit(2,1) = 'm/s' 
     194            fbdata%caddunit(2,2) = 'm/s' 
     195         ENDIF          
    172196         fbdata%cgrid(1)      = 'U'  
    173197         fbdata%cgrid(2)      = 'V' 
    174198         DO ja = 1, iadd 
    175             fbdata%caddname(1+ja) = padd%cdname(ja) 
    176             fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
    177             fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
     199            fbdata%caddname(iadd_exp+ja) = padd%cdname(ja) 
     200            fbdata%caddlong(iadd_exp+ja,1) = padd%cdlong(ja,1) 
     201            fbdata%caddunit(iadd_exp+ja,1) = padd%cdunit(ja,1) 
    178202         END DO 
    179203 
     
    246270         & ( TRIM(profdata%cvars(1)) /= 'UVEL' ) ) THEN 
    247271         CALL alloc_obfbdata( fbdata, 1, profdata%nprof, ilevel, & 
    248             &                 1 + iadd, iext, .TRUE. ) 
     272            &                 iadd_expt + iadd, iext, .TRUE. ) 
    249273         fbdata%cname(1)      = profdata%cvars(1) 
    250274         fbdata%coblong(1)    = cllongname 
     
    252276         fbdata%caddlong(1,1) = 'Model interpolated ' // TRIM(cllongname) 
    253277         fbdata%caddunit(1,1) = clunits 
     278         IF ( ln_output_clim ) THEN 
     279            fbdata%caddlong(2,1) = 'Climatological interpolated ' // TRIM(cllongname) 
     280            fbdata%caddunit(2,1) = clunits 
     281         ENDIF          
    254282         fbdata%cgrid(:)      = clgrid 
    255283         DO je = 1, iext 
     
    259287         END DO 
    260288         DO ja = 1, iadd 
    261             fbdata%caddname(1+ja) = padd%cdname(ja) 
    262             fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
    263             fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
     289            fbdata%caddname(iadd_expt+ja) = padd%cdname(ja) 
     290            fbdata%caddlong(iadd_expt+ja,1) = padd%cdlong(ja,1) 
     291            fbdata%caddunit(iadd_expt+ja,1) = padd%cdunit(ja,1) 
    264292         END DO 
    265293      ENDIF 
    266294 
    267295      fbdata%caddname(1)   = 'Hx' 
    268  
     296      IF ( ln_output_clim ) fbdata%caddname(2)   = 'CLM' 
     297       
    269298      WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc 
    270299 
     
    320349               ik = profdata%var(jvar)%nvlidx(jk) 
    321350               fbdata%padd(ik,jo,1,jvar) = profdata%var(jvar)%vmod(jk) 
     351               IF ( ln_output_clim ) THEN            
     352                  fbdata%padd(ik,jo,2,jvar) = profdata%var(jvar)%vclm(jk)      
     353               ENDIF               
    322354               fbdata%pob(ik,jo,jvar)    = profdata%var(jvar)%vobs(jk) 
    323355               fbdata%pdep(ik,jo)        = profdata%var(jvar)%vdep(jk) 
     
    334366               fbdata%iobsk(ik,jo,jvar)  = profdata%var(jvar)%mvk(jk) 
    335367               DO ja = 1, iadd 
    336                   fbdata%padd(ik,jo,1+ja,jvar) = & 
     368                  fbdata%padd(ik,jo,iadd_exp+ja,jvar) = & 
    337369                     & profdata%var(jvar)%vext(jk,padd%ipoint(ja)) 
    338370               END DO 
     
    417449      INTEGER :: je 
    418450      INTEGER :: iadd 
     451      INTEGER :: iadd_exp 
    419452      INTEGER :: iext 
    420453      INTEGER :: indx_std 
    421454      INTEGER :: iadd_std 
    422  
     455      INTEGER :: iadd_clm       
     456 
     457 
     458      ! Set up number of additional variables to be ouput: 
     459      ! Hx, CLIM, ... 
     460      iadd_exp = 1   ! Hx 
     461      IF ( ln_output_clim ) iadd_exp = iadd_exp + 1 
     462  
    423463      IF ( PRESENT( padd ) ) THEN 
    424464         iadd = padd%inum 
     
    444484      ENDIF 
    445485       
     486      iadd_clm = 0 
     487       
    446488      CALL init_obfbdata( fbdata ) 
    447489 
    448490      SELECT CASE ( TRIM(surfdata%cvars(1)) ) 
    449491      CASE('SLA') 
    450  
     492          
    451493         ! SLA needs special treatment because of MDT, so is all done here 
    452494         ! Other variables are done more generically 
     495         ! No climatology for SLA, MDT is our best estimate of that and is already output. 
    453496 
    454497         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
     
    485528         clunits    = 'Degree centigrade' 
    486529         clgrid     = 'T' 
    487  
     530         IF ( ln_output_clim ) iadd_clm = 1 
     531          
    488532      CASE('ICECONC') 
    489533 
     
    499543         clunits    = 'psu' 
    500544         clgrid     = 'T' 
    501  
     545         IF ( ln_output_clim ) iadd_clm = 1 
     546          
    502547      CASE('SLCHLTOT','LOGCHL','LogChl','logchl') 
    503548 
     
    610655       
    611656         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
    612             &                 1 + iadd + iadd_std, iext, .TRUE. ) 
     657            &                 1 + iadd + iadd_std + iadd_clm, iext, .TRUE. ) 
    613658 
    614659         fbdata%cname(1)      = surfdata%cvars(1) 
     
    628673         fbdata%cgrid(1)      = clgrid 
    629674         DO ja = 1, iadd 
    630             fbdata%caddname(1+iadd_std+ja) = padd%cdname(ja) 
    631             fbdata%caddlong(1+iadd_std+ja,1) = padd%cdlong(ja,1) 
    632             fbdata%caddunit(1+iadd_std+ja,1) = padd%cdunit(ja,1) 
     675            fbdata%caddname(1+iadd_std+iadd_clm+ja) = padd%cdname(ja) 
     676            fbdata%caddlong(1+iadd_std+iadd_clm+ja,1) = padd%cdlong(ja,1) 
     677            fbdata%caddunit(1+iadd_std+iadd_clm+ja,1) = padd%cdunit(ja,1) 
    633678         END DO 
    634679 
     
    642687         fbdata%caddunit(1+iadd_std,1) = fbdata%cobunit(1) 
    643688      ENDIF 
     689       
     690      IF ( ln_output_clim .AND. ( iadd_clm > 0 ) ) THEN 
     691         fbdata%caddname(1+iadd_std+iadd_clm)   = 'CLM' 
     692         fbdata%caddlong(1+iadd_std+iadd_clm,1) = 'Climatology' 
     693         fbdata%caddunit(1+iadd_std+iadd_clm,1) = fbdata%cobunit(1) 
     694      ENDIF 
     695       
    644696       
    645697      WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc 
     
    690742            &           krefdate = 19500101 ) 
    691743         fbdata%padd(1,jo,1,1) = surfdata%rmod(jo,1) 
    692          IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) fbdata%padd(1,jo,2,1) = surfdata%rext(jo,1) 
     744         IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) THEN 
     745            fbdata%padd(1,jo,2,1) = surfdata%rext(jo,1) 
     746         ENDIF     
     747         IF ( ln_output_clim .AND. ( iadd_clm > 0 ) ) THEN 
     748            fbdata%padd(1,jo,2,1) = surfdata%rclm(jo,1) 
     749         ENDIF 
     750                      
    693751         fbdata%pob(1,jo,1)    = surfdata%robs(jo,1)  
    694752         fbdata%pdep(1,jo)     = 0.0 
Note: See TracChangeset for help on using the changeset viewer.