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 15183 – NEMO

Changeset 15183


Ignore:
Timestamp:
2021-08-12T15:18:46+02:00 (3 years ago)
Author:
kingr
Message:

Adding changes required to write out time-averaged assimilation background. Required moving call to asm_bkg_wri from asminc.F90 to nemogcm.F90 to avoid circular USE statements.

Location:
NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/ASM/asmbkg.F90

    r15181 r15183  
    3939   USE ice 
    4040#endif 
     41   USE asminc, ONLY: ln_avgbkg 
    4142 
    4243   IMPLICIT NONE 
     
    4445    
    4546   PUBLIC   asm_bkg_wri   !: Write out the background state 
     47 
     48  !! * variables for calculating time means 
     49   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   tn_tavg  , sn_tavg   
     50   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   un_tavg  , vn_tavg 
     51   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   avt_tavg 
     52   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   en_tavg 
     53   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:)   ::   sshn_tavg 
     54   REAL(wp),SAVE :: numtimes_tavg     ! No of times to average over 
    4655 
    4756   !!---------------------------------------------------------------------- 
     
    7180      INTEGER :: inum          ! File unit number 
    7281      REAL(wp) :: zdate        ! Date 
     82      INTEGER :: ierror 
    7383      !!----------------------------------------------------------------------- 
    7484 
    75       !                                !------------------------------------------- 
    76       IF( kt == nitbkg_r ) THEN        ! Write out background at time step nitbkg_r 
    77          !                             !-----------------------------------======== 
     85      ! If creating an averaged assim bkg, initialise on first timestep 
     86      IF ( ln_avgbkg .AND. kt == ( nn_it000 - 1) ) THEN 
     87         ! Allocate memory  
     88         ALLOCATE( tn_tavg(jpi,jpj,jpk), STAT=ierror ) 
     89         IF( ierror > 0 ) THEN 
     90            CALL ctl_stop( 'asm_bkg_wri: unable to allocate tn_tavg' )   ;   RETURN 
     91         ENDIF 
     92         tn_tavg(:,:,:)=0 
     93         ALLOCATE( sn_tavg(jpi,jpj,jpk), STAT=ierror ) 
     94         IF( ierror > 0 ) THEN 
     95            CALL ctl_stop( 'asm_bkg_wri: unable to allocate sn_tavg' )   ;   RETURN 
     96         ENDIF 
     97         sn_tavg(:,:,:)=0 
     98         ALLOCATE( un_tavg(jpi,jpj,jpk), STAT=ierror ) 
     99         IF( ierror > 0 ) THEN 
     100            CALL ctl_stop( 'asm_bkg_wri: unable to allocate un_tavg' )   ;   RETURN 
     101         ENDIF 
     102         un_tavg(:,:,:)=0 
     103         ALLOCATE( vn_tavg(jpi,jpj,jpk), STAT=ierror ) 
     104         IF( ierror > 0 ) THEN 
     105            CALL ctl_stop( 'asm_bkg_wri: unable to allocate vn_tavg' )   ;   RETURN 
     106         ENDIF 
     107         vn_tavg(:,:,:)=0 
     108         ALLOCATE( sshn_tavg(jpi,jpj), STAT=ierror ) 
     109         IF( ierror > 0 ) THEN 
     110            CALL ctl_stop( 'asm_bkg_wri: unable to allocate sshn_tavg' )   ;   RETURN 
     111         ENDIF 
     112         sshn_tavg(:,:)=0 
     113         IF( ln_zdftke ) THEN 
     114            ALLOCATE( en_tavg(jpi,jpj,jpk), STAT=ierror ) 
     115            IF( ierror > 0 ) THEN 
     116               CALL ctl_stop( 'asm_bkg_wri: unable to allocate en_tavg' )   ;   RETURN 
     117            ENDIF 
     118            en_tavg(:,:,:)=0 
     119         ENDIF 
     120         ALLOCATE( avt_tavg(jpi,jpj,jpk), STAT=ierror ) 
     121         IF( ierror > 0 ) THEN 
     122            CALL ctl_stop( 'asm_bkg_wri: unable to allocate avt_tavg' )   ;   RETURN 
     123         ENDIF 
     124         avt_tavg(:,:,:)=0 
     125          
     126         numtimes_tavg = REAL ( nitavgbkg_r -  nn_it000 + 1 ) 
     127      ENDIF 
     128       
     129      ! If creating an averaged assim bkg, sum the contribution every timestep 
     130      IF ( ln_avgbkg ) THEN  
     131         IF (lwp) THEN 
     132              WRITE(numout,*) 'asm_bkg_wri : Summing assim bkg fields at timestep ',kt 
     133              WRITE(numout,*) '~~~~~~~~~~~~ ' 
     134         ENDIF 
     135 
     136         tn_tavg(:,:,:)  = tn_tavg(:,:,:)  + tsn(:,:,:,jp_tem) / numtimes_tavg 
     137         sn_tavg(:,:,:)  = sn_tavg(:,:,:)  + tsn(:,:,:,jp_sal) / numtimes_tavg 
     138         sshn_tavg(:,:)  = sshn_tavg(:,:)  + sshn (:,:)        / numtimes_tavg 
     139         un_tavg(:,:,:)  = un_tavg(:,:,:)  + un(:,:,:)         / numtimes_tavg 
     140         vn_tavg(:,:,:)  = vn_tavg(:,:,:)  + vn(:,:,:)         / numtimes_tavg 
     141         avt_tavg(:,:,:) = avt_tavg(:,:,:) + avt(:,:,:)        / numtimes_tavg 
     142         IF( ln_zdftke ) THEN 
     143            en_tavg(:,:,:)  = en_tavg(:,:,:)  + en(:,:,:)         / numtimes_tavg 
     144         ENDIF 
     145      ENDIF 
     146      
     147 
     148      ! Write out background at time step nitbkg_r or nitavgbkg_r 
     149      IF ( ( .NOT. ln_avgbkg .AND. (kt == nitbkg_r) ) .OR. & 
     150      &          ( ln_avgbkg .AND. (kt == nitavgbkg_r) ) ) THEN 
    78151         ! 
    79152         WRITE(cl_asmbkg, FMT='(A,".nc")' ) TRIM( c_asmbkg ) 
     
    87160            CALL iom_open( c_asmbkg, inum, ldwrt = .TRUE. ) 
    88161            ! 
    89             IF( nitbkg_r == nit000 - 1 ) THEN      ! Treat special case when nitbkg = 0 
    90                zdate = REAL( ndastp ) 
    91                IF( ln_zdftke ) THEN                   ! read turbulent kinetic energy ( en ) 
    92                   IF(lwp) WRITE(numout,*) ' Reading TKE (en) from restart...' 
    93                   CALL tke_rst( nit000, 'READ' ) 
     162            !                                      ! Write the information 
     163            IF ( ln_avgbkg ) THEN 
     164               IF( nitavgbkg_r == nit000 - 1 ) THEN      ! Treat special case when nitavgbkg = 0 
     165                  zdate = REAL( ndastp ) 
     166                  IF( ln_zdftke ) THEN 
     167                     IF(lwp) WRITE(numout,*) ' Reading TKE (en) from restart...' 
     168                     CALL tke_rst( nit000, 'READ' )               ! lk_zdftke=T :   Read turbulent kinetic energy ( en ) 
     169                  ENDIF 
     170               ELSE 
     171                  zdate = REAL( ndastp ) 
    94172               ENDIF 
     173               CALL iom_rstput( kt, nitavgbkg_r, inum, 'rdastp' , zdate   ) 
     174               CALL iom_rstput( kt, nitavgbkg_r, inum, 'un'     , un_tavg ) 
     175               CALL iom_rstput( kt, nitavgbkg_r, inum, 'vn'     , vn_tavg ) 
     176               CALL iom_rstput( kt, nitavgbkg_r, inum, 'tn'     , tn_tavg ) 
     177               CALL iom_rstput( kt, nitavgbkg_r, inum, 'sn'     , sn_tavg ) 
     178               CALL iom_rstput( kt, nitavgbkg_r, inum, 'sshn'   , sshn_tavg) 
     179               IF( ln_zdftke ) CALL iom_rstput( kt, nitavgbkg_r, inum, 'en'     , en_tavg ) 
     180               CALL iom_rstput( kt, nitavgbkg_r, inum, 'avt'    , avt_tavg) 
     181               ! 
    95182            ELSE 
    96                zdate = REAL( ndastp ) 
    97             ENDIF 
    98             ! 
    99             !                                      ! Write the information 
    100             CALL iom_rstput( kt, nitbkg_r, inum, 'rdastp' , zdate             ) 
    101             CALL iom_rstput( kt, nitbkg_r, inum, 'un'     , un                ) 
    102             CALL iom_rstput( kt, nitbkg_r, inum, 'vn'     , vn                ) 
    103             CALL iom_rstput( kt, nitbkg_r, inum, 'tn'     , tsn(:,:,:,jp_tem) ) 
    104             CALL iom_rstput( kt, nitbkg_r, inum, 'sn'     , tsn(:,:,:,jp_sal) ) 
    105             CALL iom_rstput( kt, nitbkg_r, inum, 'sshn'   , sshn              ) 
    106             IF( ln_zdftke )   CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en                ) 
    107             CALL iom_rstput( kt, nitbkg_r, inum, 'avt'    , avt               ) 
    108             ! 
    109             CALL iom_close( inum ) 
     183               IF( nitbkg_r == nit000 - 1 ) THEN      ! Treat special case when nitbkg = 0 
     184                  zdate = REAL( ndastp ) 
     185                  IF( ln_zdftke ) THEN 
     186                     IF(lwp) WRITE(numout,*) ' Reading TKE (en) from restart...' 
     187                     CALL tke_rst( nit000, 'READ' )               ! lk_zdftke=T :   Read turbulent kinetic energy ( en ) 
     188                  ENDIF 
     189               ELSE 
     190                  zdate = REAL( ndastp ) 
     191               ENDIF 
     192               CALL iom_rstput( kt, nitbkg_r, inum, 'rdastp' , zdate             ) 
     193               CALL iom_rstput( kt, nitbkg_r, inum, 'un'     , un                ) 
     194               CALL iom_rstput( kt, nitbkg_r, inum, 'vn'     , vn                ) 
     195               CALL iom_rstput( kt, nitbkg_r, inum, 'tn'     , tsn(:,:,:,jp_tem) ) 
     196               CALL iom_rstput( kt, nitbkg_r, inum, 'sn'     , tsn(:,:,:,jp_sal) ) 
     197               CALL iom_rstput( kt, nitbkg_r, inum, 'sshn'   , sshn              ) 
     198               IF( ln_zdftke )  CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en                ) 
     199               CALL iom_rstput( kt, nitbkg_r, inum, 'avt'    , avt               ) 
     200               ! 
     201               CALL iom_close( inum ) 
     202            ENDIF 
    110203         ENDIF 
    111204         ! 
  • NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/ASM/asminc.F90

    r14075 r15183  
    3131   USE zpshde          ! Partial step : Horizontal Derivative 
    3232   USE asmpar          ! Parameters for the assmilation interface 
    33    USE asmbkg          !  
    3433   USE c1d             ! 1D initialization 
    3534   USE sbc_oce         ! Surface boundary condition variables. 
     
    5958#endif 
    6059   LOGICAL, PUBLIC :: ln_bkgwri     !: No output of the background state fields 
     60   LOGICAL, PUBLIC :: ln_avgbkg     !: No output of the mean background state fields 
    6161   LOGICAL, PUBLIC :: ln_asmiau     !: No applying forcing with an assimilation increment 
    6262   LOGICAL, PUBLIC :: ln_asmdin     !: No direct initialization 
     
    8282   INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval  
    8383   INTEGER , PUBLIC ::   nitiaufin   !: Time step of the end of the IAU interval 
     84   INTEGER , PUBLIC ::   nitavgbkg   !: Number of timesteps to average assim bkg [0,nitavgbkg] 
    8485   !  
    8586   INTEGER , PUBLIC ::   niaufn      !: Type of IAU weighing function: = 0   Constant weighting 
     
    122123      REAL(KIND=dp) :: ditiaustr_date  ! Date YYYYMMDD.HHMMSS of IAU interval start time step 
    123124      REAL(KIND=dp) :: ditiaufin_date  ! Date YYYYMMDD.HHMMSS of IAU interval final time step 
     125      REAL(KIND=dp) :: ditavgbkg_date  ! Date YYYYMMDD.HHMMSS of end of assim bkg averaging period 
    124126 
    125127      REAL(wp) :: znorm        ! Normalization factor for IAU weights 
     
    132134      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zhdiv   ! 2D workspace 
    133135      !! 
    134       NAMELIST/nam_asminc/ ln_bkgwri,                                      & 
     136      NAMELIST/nam_asminc/ ln_bkgwri, ln_avgbkg,                           & 
    135137         &                 ln_trainc, ln_dyninc, ln_sshinc,                & 
    136138         &                 ln_asmdin, ln_asmiau,                           & 
    137139         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   & 
    138          &                 ln_salfix, salfixmin, nn_divdmp 
     140         &                 ln_salfix, salfixmin, nn_divdmp, nitavgbkg 
    139141      !!---------------------------------------------------------------------- 
    140142 
     
    142144      ! Read Namelist nam_asminc : assimilation increment interface 
    143145      !----------------------------------------------------------------------- 
    144       ln_seaiceinc   = .FALSE. 
     146      ! Set default values 
     147      ln_bkgwri = .FALSE. 
     148      ln_avgbkg = .FALSE. 
     149      ln_trainc = .FALSE. 
     150      ln_dyninc = .FALSE. 
     151      ln_sshinc = .FALSE. 
     152      ln_asmdin = .FALSE. 
     153      ln_asmiau = .TRUE. 
     154      ln_salfix = .FALSE. 
     155      ln_seaiceinc = .FALSE. 
    145156      ln_temnofreeze = .FALSE. 
     157      nitbkg    = 0 
     158      nitdin    = 0       
     159      nitiaustr = 1 
     160      nitiaufin = 150 
     161      niaufn    = 0 
     162      salfixmin = -9999 
     163      nitavgbkg = 1 
    146164 
    147165      REWIND( numnam_ref )              ! Namelist nam_asminc in reference namelist : Assimilation increment 
     
    160178         WRITE(numout,*) '   Namelist namasm : set assimilation increment parameters' 
    161179         WRITE(numout,*) '      Logical switch for writing out background state          ln_bkgwri = ', ln_bkgwri 
     180         WRITE(numout,*) '      Logical switch for writing mean background state         ln_avgbkg = ', ln_avgbkg 
    162181         WRITE(numout,*) '      Logical switch for applying tracer increments            ln_trainc = ', ln_trainc 
    163182         WRITE(numout,*) '      Logical switch for applying velocity increments          ln_dyninc = ', ln_dyninc 
     
    170189         WRITE(numout,*) '      Timestep of start of IAU interval in [0,nitend-nit000-1] nitiaustr = ', nitiaustr 
    171190         WRITE(numout,*) '      Timestep of end of IAU interval in [0,nitend-nit000-1]   nitiaufin = ', nitiaufin 
     191         WRITE(numout,*) '      Number of timesteps to average assim bkg [0,nitavgbkg]   nitavgbkg = ', nitavgbkg 
    172192         WRITE(numout,*) '      Type of IAU weighting function                           niaufn    = ', niaufn 
    173193         WRITE(numout,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', ln_salfix 
     
    179199      nitiaustr_r = nitiaustr + nit000 - 1            ! Start of IAU interval referenced to nit000 
    180200      nitiaufin_r = nitiaufin + nit000 - 1            ! End of IAU interval referenced to nit000 
     201      nitavgbkg_r = nitavgbkg + nit000 - 1            ! Averaging period referenced to nit000 
    181202 
    182203      iiauper     = nitiaufin_r - nitiaustr_r + 1     ! IAU interval length 
     
    188209      CALL calc_date( nitiaustr_r, ditiaustr_date )   ! IAU start time referenced to ndate0 
    189210      CALL calc_date( nitiaufin_r, ditiaufin_date )   ! IAU end time referenced to ndate0 
     211      CALL calc_date( nitavgbkg_r, ditavgbkg_date )   ! End of assim bkg averaging period referenced to ndate0 
    190212 
    191213      IF(lwp) THEN 
     
    199221         WRITE(numout,*) '       nitiaustr_r = ', nitiaustr_r 
    200222         WRITE(numout,*) '       nitiaufin_r = ', nitiaufin_r 
     223         WRITE(numout,*) '       nitavgbkg_r = ', nitavgbkg_r 
    201224         WRITE(numout,*) 
    202225         WRITE(numout,*) '   Dates referenced to current cycle:' 
     
    209232         WRITE(numout,*) '       ditiaustr_date = ', ditiaustr_date 
    210233         WRITE(numout,*) '       ditiaufin_date = ', ditiaufin_date 
     234         WRITE(numout,*) '       ditavgbkg_date = ', ditavgbkg_date 
    211235      ENDIF 
    212236 
     
    249273         &                ' the cycle interval') 
    250274 
     275      IF ( nitavgbkg_r > nitend ) & 
     276         & CALL ctl_stop( ' nitavgbkg_r :',  & 
     277         &                ' Assim bkg averaging period is outside', & 
     278         &                ' the cycle interval') 
     279       
    251280      IF ( nstop > 0 ) RETURN       ! if there are any errors then go no further 
    252281 
     
    494523      ! 
    495524      IF( lk_asminc ) THEN                            !==  data assimilation  ==! 
    496          IF( ln_bkgwri )   CALL asm_bkg_wri( nit000 - 1 )      ! Output background fields 
     525!         IF( ln_bkgwri )   CALL asm_bkg_wri( nit000 - 1 )      ! Output background fields 
    497526         IF( ln_asmdin ) THEN                                  ! Direct initialization 
    498527            IF( ln_trainc )   CALL tra_asm_inc( nit000 - 1 )      ! Tracers 
  • NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/ASM/asmpar.F90

    r14075 r15183  
    2121   INTEGER, PUBLIC ::   nitiaufin_r   !: IAU final time step referenced to nit000 
    2222   INTEGER, PUBLIC ::   nittrjfrq     !: Frequency of trajectory output for 4D-VAR 
     23   INTEGER, PUBLIC ::   nitavgbkg_r   !: Averaging period for assim bkg referenced to nit000 
    2324 
    2425   !!---------------------------------------------------------------------- 
  • NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/nemogcm.F90

    r14075 r15183  
    166166      END DO 
    167167#else 
     168      IF( lk_asminc .AND. ln_bkgwri ) CALL asm_bkg_wri( nit000 - 1 )    ! Output background fields 
    168169      ! 
    169170# if defined key_agrif 
  • NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/step.F90

    r14075 r15183  
    296296#endif 
    297297 
     298      IF( ln_bkgwri )    CALL asm_bkg_wri  ( kstp )      ! output background fields 
    298299      IF( ln_diaobs  )   CALL dia_obs      ( kstp )      ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 
    299300 
Note: See TracChangeset for help on using the changeset viewer.