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.
asmbkg.F90 in NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/ASM – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_CO9_OBS_ASM/src/OCE/ASM/asmbkg.F90 @ 15183

Last change on this file since 15183 was 15183, checked in by kingr, 3 years ago

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.
File size: 12.2 KB
Line 
1MODULE asmbkg
2   !!======================================================================
3   !!                       ***  MODULE asmtrj -> asmbkg  ***
4   !! Assimilation trajectory interface: Write to file the background state and the model state trajectory
5   !!======================================================================
6   !! History :       ! 2007-03 (M. Martin)  Met. Office version
7   !!                 ! 2007-04 (A. Weaver)  asm_trj_wri, original code
8   !!                 ! 2007-03 (K. Mogensen)  Adapt to NEMOVAR and use IOM instead of IOIPSL
9   !!                 ! 2007-04 (A. Weaver)  Name change (formally asmbkg.F90). Distinguish
10   !!                                        background states in Jb term and at analysis time.
11   !!                                        Include state trajectory routine (currently empty)
12   !!                 ! 2007-07 (A. Weaver)  Add tke_rst and flt_rst for case nitbkg=0
13   !!                 ! 2009-03 (F. Vigilant)  Add hmlp (zdfmxl) for no tracer nmldp=2
14   !!                 ! 2009-06 (F. Vigilant) asm_trj_wri: special case when kt=nit000-1
15   !!                 ! 2009-07 (F. Vigilant) asm_trj_wri: add computation of eiv at restart
16   !!                 ! 2010-01 (A. Vidard) split asm_trj_wri into tam_trj_wri and asm_bkg_wri
17   !!----------------------------------------------------------------------
18
19   !!----------------------------------------------------------------------
20   !!   asm_bkg_wri  : Write out the background state
21   !!   asm_trj_wri  : Write out the model state trajectory (used with 4D-Var)
22   !!----------------------------------------------------------------------
23   USE oce                ! Dynamics and active tracers defined in memory
24   USE sbc_oce            ! Ocean surface boundary conditions
25   USE zdf_oce            ! Vertical mixing variables
26   USE zdfddm             ! Double diffusion mixing parameterization
27   USE ldftra             ! Lateral diffusion: eddy diffusivity coefficients
28   USE ldfslp             ! Lateral diffusion: slopes of neutral surfaces
29   USE tradmp             ! Tracer damping
30   USE zdftke             ! TKE vertical physics
31   USE eosbn2             ! Equation of state (eos_bn2 routine)
32   USE zdfmxl             ! Mixed layer depth
33   USE dom_oce     , ONLY :   ndastp
34   USE in_out_manager     ! I/O manager
35   USE iom                ! I/O module
36   USE asmpar             ! Parameters for the assmilation interface
37   USE zdfmxl             ! mixed layer depth
38#if defined key_si3
39   USE ice
40#endif
41   USE asminc, ONLY: ln_avgbkg
42
43   IMPLICIT NONE
44   PRIVATE
45   
46   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
55
56   !!----------------------------------------------------------------------
57   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
58   !! $Id$
59   !! Software governed by the CeCILL license (see ./LICENSE)
60   !!----------------------------------------------------------------------
61CONTAINS
62
63   SUBROUTINE asm_bkg_wri( kt )
64      !!-----------------------------------------------------------------------
65      !!                  ***  ROUTINE asm_bkg_wri ***
66      !!
67      !! ** Purpose : Write to file the background state for later use in the
68      !!              inner loop of data assimilation or for direct initialization
69      !!              in the outer loop.
70      !!
71      !! ** Method  : Write out the background state for use in the Jb term
72      !!              in the cost function and for use with direct initialization
73      !!              at analysis time.
74      !!-----------------------------------------------------------------------
75      INTEGER, INTENT( IN ) :: kt               ! Current time-step
76      !
77      CHARACTER (LEN=50) :: cl_asmbkg
78      CHARACTER (LEN=50) :: cl_asmdin
79      LOGICAL :: llok          ! Check if file exists
80      INTEGER :: inum          ! File unit number
81      REAL(wp) :: zdate        ! Date
82      INTEGER :: ierror
83      !!-----------------------------------------------------------------------
84
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
151         !
152         WRITE(cl_asmbkg, FMT='(A,".nc")' ) TRIM( c_asmbkg )
153         cl_asmbkg = TRIM( cl_asmbkg )
154         INQUIRE( FILE = cl_asmbkg, EXIST = llok )
155         !
156         IF( .NOT. llok ) THEN
157            IF(lwp) WRITE(numout,*) ' Setting up assimilation background file '// TRIM( c_asmbkg )
158            !
159            !                                      ! Define the output file       
160            CALL iom_open( c_asmbkg, inum, ldwrt = .TRUE. )
161            !
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 )
172               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               !
182            ELSE
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
203         ENDIF
204         !
205      ENDIF
206
207      !                                !-------------------------------------------
208      IF( kt == nitdin_r ) THEN        ! Write out background at time step nitdin_r
209         !                             !-----------------------------------========
210         !
211         WRITE(cl_asmdin, FMT='(A,".nc")' ) TRIM( c_asmdin )
212         cl_asmdin = TRIM( cl_asmdin )
213         INQUIRE( FILE = cl_asmdin, EXIST = llok )
214         !
215         IF( .NOT. llok ) THEN
216            IF(lwp) WRITE(numout,*) ' Setting up assimilation background file '// TRIM( c_asmdin )
217            !
218            !                                      ! Define the output file       
219            CALL iom_open( c_asmdin, inum, ldwrt = .TRUE. )
220            !
221            IF( nitdin_r == nit000 - 1 ) THEN      ! Treat special case when nitbkg = 0
222
223               zdate = REAL( ndastp )
224            ELSE
225               zdate = REAL( ndastp )
226            ENDIF
227            !
228            !                                      ! Write the information
229            CALL iom_rstput( kt, nitdin_r, inum, 'rdastp' , zdate             )
230            CALL iom_rstput( kt, nitdin_r, inum, 'un'     , un                )
231            CALL iom_rstput( kt, nitdin_r, inum, 'vn'     , vn                )
232            CALL iom_rstput( kt, nitdin_r, inum, 'tn'     , tsn(:,:,:,jp_tem) )
233            CALL iom_rstput( kt, nitdin_r, inum, 'sn'     , tsn(:,:,:,jp_sal) )
234            CALL iom_rstput( kt, nitdin_r, inum, 'sshn'   , sshn              )
235#if defined key_si3
236            IF( nn_ice == 2 ) THEN
237               IF( ALLOCATED(at_i) ) THEN
238                  CALL iom_rstput( kt, nitdin_r, inum, 'iceconc', at_i(:,:)   )
239               ELSE
240                  CALL ctl_warn('asm_bkg_wri: Ice concentration not written to background ',   &
241                     &          'as ice variable at_i not allocated on this timestep')
242               ENDIF
243            ENDIF
244#endif
245            !
246            CALL iom_close( inum )
247         ENDIF
248         !
249      ENDIF
250      !                   
251   END SUBROUTINE asm_bkg_wri
252
253   !!======================================================================
254END MODULE asmbkg
Note: See TracBrowser for help on using the repository browser.