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 5063 for branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/SAO_SRC/sao_read.F90 – NEMO

Ignore:
Timestamp:
2015-02-05T17:29:55+01:00 (9 years ago)
Author:
andrewryan
Message:

gross simplification of stand alone observation operator

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/SAO_SRC/sao_read.F90

    r4849 r5063  
    1  
    21MODULE sao_read 
    32   !!================================================================== 
    … …  
    54   !! Read routines : I/O for Stand Alone Observation operator 
    65   !!================================================================== 
    7  
    86   USE mppini 
    97   USE lib_mpp 
    … …  
    2826      !! 
    2927      !! Purpose : To choose appropriate read method 
    30       !! Method  :  
     28      !! Method  : 
    3129      !! 
    3230      !! Author  : A. Ryan Oct 2013 
    … …  
    3634              & kfile         !: File number 
    3735      CHARACTER(len=lc) :: & 
    38               & cdfilename, & !: File name 
    39               & cmatchname    !: Match name 
     36              & cdfilename    !: File name 
    4037      INTEGER :: & 
    41               & kindex       !: File index to read 
    42   
    43       !! Filename, index and match-up kind 
     38              & kindex        !: File index to read 
     39 
    4440      cdfilename = TRIM(sao_files(kfile)) 
    45       cmatchname = TRIM(cl4_vars(kfile)) 
    4641      kindex = nn_sao_idx(kfile) 
    47  
    48       !! Update model fields 
    49       !! Class 4 variables: forecast, persistence, 
    50       !!                    nrt_analysis, best_estimate 
    51       !! Feedback variables: empty string 
    52       IF ( (TRIM(cmatchname) == 'forecast') .OR. & 
    53          & (TRIM(cmatchname) == 'persistence') .OR. & 
    54          & (TRIM(cmatchname) == 'nrt_analysis') .OR. & 
    55          & (TRIM(cmatchname) == 'best_estimate').OR. & 
    56          & (TRIM(cmatchname) == '') ) THEN 
    57          CALL sao_read_file(TRIM(cdfilename), kindex) 
    58          CALL sao_read_juld(TRIM(cdfilename), kindex, cl4_modjuld) 
    59       ELSE IF (TRIM(cmatchname) == 'climatology') THEN 
    60          WRITE(numout,*) 'Interpolating climatologies' 
    61       ELSE IF (TRIM(cmatchname) == 'altimeter') THEN 
    62          CALL sao_read_altbias(TRIM(cdfilename)) 
    63          CALL sao_read_juld(TRIM(cdfilename), kindex, cl4_modjuld) 
    64       END IF 
     42      CALL sao_read_file(TRIM(cdfilename), kindex) 
    6543 
    6644   END SUBROUTINE sao_rea_dri 
    67  
    68    SUBROUTINE sao_read_altbias(filename) 
    69       !!------------------------------------------------------------------------ 
    70       !!                      *** sao_read_altbias *** 
    71       !! 
    72       !! Purpose : To read altimeter bias and set tn,sn to missing values 
    73       !! Method  : Use subdomain indices to create start and count matrices 
    74       !!           for netcdf read. 
    75       !! 
    76       !! Author  : A. Ryan Sept 2012 
    77       !!------------------------------------------------------------------------ 
    78       CHARACTER(len=*), INTENT(IN) :: filename 
    79       INTEGER                      :: ncid, & 
    80                                     & varid,& 
    81                                     & istat,& 
    82                                     & ntimes,& 
    83                                     & tdim, & 
    84                                     & xdim, & 
    85                                     & ydim, & 
    86                                     & zdim 
    87       INTEGER                      :: ii, ij, ik 
    88       INTEGER, DIMENSION(3)        :: start_s, & 
    89                                     & count_s 
    90       REAL(fbdp), DIMENSION(:,:),  ALLOCATABLE :: temp_sshn 
    91       REAL(fbdp)                     :: fill_val 
    92  
    93       IF (TRIM(filename) == 'nofile') THEN 
    94          tsn(:,:,:,:) = fbrmdi 
    95          sshn(:,:) = fbrmdi  
    96       ELSE 
    97          ! Open Netcdf file to find dimension id 
    98          istat = nf90_open(trim(filename),nf90_nowrite,ncid) 
    99          istat = nf90_inq_dimid(ncid,'x',xdim) 
    100          istat = nf90_inq_dimid(ncid,'y',ydim) 
    101          istat = nf90_inq_dimid(ncid,'deptht',zdim) 
    102          istat = nf90_inq_dimid(ncid,'time',tdim) 
    103          istat = nf90_inquire_dimension(ncid, tdim, len=ntimes) 
    104          ! Allocate temporary temperature array 
    105          ALLOCATE(temp_sshn(nlci,nlcj)) 
    106          ! Create start and count arrays 
    107          start_s = (/ nimpp, njmpp, 1 /) 
    108          count_s = (/ nlci,  nlcj,  1 /) 
    109            
    110          ! Altimeter bias 
    111          istat = nf90_inq_varid(ncid,'altbias',varid) 
    112          istat = nf90_get_att(ncid, varid, '_FillValue', fill_val) 
    113          istat = nf90_get_var(ncid, varid, temp_sshn, start_s, count_s) 
    114          WHERE(temp_sshn(:,:) == fill_val) temp_sshn(:,:) = fbrmdi 
    115     
    116          ! Initialise tsn, sshn to fbrmdi 
    117          tsn(:,:,:,:) = fbrmdi 
    118          sshn(:,:) = fbrmdi  
    119  
    120          ! Fill sshn with altimeter bias  
    121          sshn(1:nlci,1:nlcj) = temp_sshn(:,:) * tmask(1:nlci,1:nlcj,1) 
    122  
    123          ! Remove halo from tmask, sshn to prevent double obs counting 
    124          IF (jpi > nlci) THEN 
    125              tmask(nlci+1:,:,:) = 0 
    126              sshn(nlci+1:,:) = 0 
    127          END IF 
    128          IF (jpj > nlcj) THEN 
    129              tmask(:,nlcj+1:,:) = 0 
    130              sshn(:,nlcj+1:) = 0 
    131          END IF 
    132  
    133          ! Deallocate arrays 
    134          DEALLOCATE(temp_sshn) 
    135  
    136          ! Close netcdf file 
    137          istat = nf90_close(ncid) 
    138       END IF 
    139     
    140    END SUBROUTINE sao_read_altbias 
    14145 
    14246   SUBROUTINE sao_read_file(filename, ifcst) 
    … …  
    17680      IF (TRIM(filename) == 'nofile') THEN 
    17781         tsn(:,:,:,:) = fbrmdi 
    178          sshn(:,:) = fbrmdi  
     82         sshn(:,:) = fbrmdi 
    17983      ELSE 
    18084         WRITE(numout,*) "Opening :", TRIM(filename) 
    … …  
    19094         istat = nf90_inq_dimid(ncid,'time_counter',tdim) 
    19195         istat = nf90_inquire_dimension(ncid, tdim, len=ntimes) 
    192          IF (ifcst .LE. ntimes) THEN    
     96         IF (ifcst .LE. ntimes) THEN 
    19397            ! Allocate temporary temperature array 
    19498            ALLOCATE(temp_tn(nlci,nlcj,jpk)) 
    19599            ALLOCATE(temp_sn(nlci,nlcj,jpk)) 
    196100            ALLOCATE(temp_sshn(nlci,nlcj)) 
    197        
     101 
    198102            ! Set temp_tn, temp_sn to 0. 
    199103            temp_tn(:,:,:) = fbrmdi 
    200104            temp_sn(:,:,:) = fbrmdi 
    201105            temp_sshn(:,:) = fbrmdi 
    202        
     106 
    203107            ! Create start and count arrays 
    204108            start_n = (/ nimpp, njmpp, 1,   ifcst /) 
    … …  
    206110            start_s = (/ nimpp, njmpp, ifcst /) 
    207111            count_s = (/ nlci,  nlcj,  1     /) 
    208               
     112 
    209113            ! Read information into temporary arrays 
    210114            ! retrieve varid and read in temperature 
    … …  
    213117            istat = nf90_get_var(ncid, varid, temp_tn, start_n, count_n) 
    214118            WHERE(temp_tn(:,:,:) == fill_val) temp_tn(:,:,:) = fbrmdi 
    215        
     119 
    216120            ! retrieve varid and read in salinity 
    217121            istat = nf90_inq_varid(ncid,'vosaline',varid) 
    … …  
    219123            istat = nf90_get_var(ncid, varid, temp_sn, start_n, count_n) 
    220124            WHERE(temp_sn(:,:,:) == fill_val) temp_sn(:,:,:) = fbrmdi 
    221        
     125 
    222126            ! retrieve varid and read in SSH 
    223127            istat = nf90_inq_varid(ncid,'sossheig',varid) 
    … …  
    226130               istat = nf90_inq_varid(ncid,'altbias',varid) 
    227131            END IF 
    228        
     132 
    229133            istat = nf90_get_att(ncid, varid, '_FillValue', fill_val) 
    230134            istat = nf90_get_var(ncid, varid, temp_sshn, start_s, count_s) 
    231135            WHERE(temp_sshn(:,:) == fill_val) temp_sshn(:,:) = fbrmdi 
    232     
     136 
    233137            ! Initialise tsn, sshn to fbrmdi 
    234138            tsn(:,:,:,:) = fbrmdi 
    235             sshn(:,:) = fbrmdi  
     139            sshn(:,:) = fbrmdi 
    236140 
    237141            ! Mask out missing data index 
    … …  
    256160            ! Deallocate arrays 
    257161            DEALLOCATE(temp_tn, temp_sn, temp_sshn) 
    258          ELSE    
     162         ELSE 
    259163            ! Mark all as missing data 
    260164            tsn(:,:,:,:) = fbrmdi 
    261             sshn(:,:) = fbrmdi  
     165            sshn(:,:) = fbrmdi 
    262166         ENDIF 
    263167         ! Close netcdf file 
    … …  
    266170      END IF 
    267171   END SUBROUTINE sao_read_file 
    268  
    269    SUBROUTINE sao_read_juld(filename, ifcst, julian) 
    270       USE calendar 
    271       !!-------------------------------------------------------------------- 
    272       !!                 *** sao_read_juld *** 
    273       !! 
    274       !!   Purpose : To read model julian day information from file 
    275       !!   Author  : A. Ryan Nov 2010 
    276       !!-------------------------------------------------------------------- 
    277  
    278       !! Routine arguments 
    279       CHARACTER(len=*), INTENT(IN)  :: filename 
    280       INTEGER,          INTENT(IN)  :: ifcst 
    281       REAL,             INTENT(OUT) :: julian    !: Julian day 
    282  
    283       !! Local variables 
    284       INTEGER :: year,  &   !: Date information 
    285                & month, & 
    286                & day,   & 
    287                & hour,  & 
    288                & minute,& 
    289                & second 
    290       INTEGER :: istat, &   !: Netcdf variables 
    291                & ncid,  & 
    292                & dimid, & 
    293                & varid, & 
    294                & ntime       
    295       REAL,DIMENSION(:),ALLOCATABLE :: r_sec     !: Remainder seconds 
    296       CHARACTER(len=120) :: time_str  !: time string 
    297  
    298       time_str='' 
    299       !! Read in time_counter and remainder seconds 
    300       istat = nf90_open(trim(filename),nf90_nowrite,ncid) 
    301       istat = nf90_inq_dimid(ncid,'time_counter',dimid) 
    302       IF (istat /= nf90_noerr) THEN 
    303           istat = nf90_inq_dimid(ncid,'time',dimid) 
    304       ENDIF 
    305       istat = nf90_inquire_dimension(ncid,dimid,len=ntime) 
    306       istat = nf90_inq_varid(ncid,'time_counter',varid) 
    307       IF (istat /= nf90_noerr) THEN 
    308           istat = nf90_inq_dimid(ncid,'time',dimid) 
    309       ENDIF 
    310       istat = nf90_get_att(ncid,varid,'units',time_str) 
    311       ALLOCATE(r_sec(ntime)) 
    312       istat = nf90_get_var(ncid,varid, r_sec) 
    313       istat = nf90_close(ncid) 
    314  
    315       !! Fill yyyy-mm-dd hh:mm:ss 
    316       !! format(('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2)) 
    317       100 format((14x, I4.4,1x,I2.2,1x,I2.2,1x,I2.2,1x,I2.2,1x,I2.2)) 
    318       READ( time_str, 100 ) year, month, day, hour, minute, second 
    319  
    320       CALL ymds2ju(year, month, day, r_sec(ifcst), julian) 
    321  
    322       !! To take a comment from the ymds2ju subroutine 
    323  
    324    !- In the case of the Gregorian calendar we have chosen to use 
    325    !- the Lilian day numbers. This is the day counter which starts 
    326    !- on the 15th October 1582. 
    327    !- This is the day at which Pope Gregory XIII introduced the 
    328    !- Gregorian calendar. 
    329    !- Compared to the true Julian calendar, which starts some 
    330    !- 7980 years ago, the Lilian days are smaler and are dealt with 
    331    !- easily on 32 bit machines. With the true Julian days you can only 
    332    !- the fraction of the day in the real part to a precision of 
    333    !- a 1/4 of a day with 32 bits. 
    334        
    335       !! The obs operator routines prefer to calculate Julian days since  
    336       !! 01/01/1950 00:00:00 
    337       !! In order to convert to the 1950 version we must adjust by the number 
    338       !! of days between 15th October 1582 and 1st Jan 1950 
    339  
    340       julian = julian - 134123. 
    341        
    342       DEALLOCATE(r_sec) 
    343        
    344    END SUBROUTINE sao_read_juld 
    345  
    346 END MODULE sao_read  
    347  
     172END MODULE sao_read 
Note: See TracChangeset for help on using the changeset viewer.