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 branches/UKMO/dev_r5518_25hr_mean_assim_bkg/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

source: branches/UKMO/dev_r5518_25hr_mean_assim_bkg/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90 @ 6762

Last change on this file since 6762 was 6762, checked in by kingr, 8 years ago

Added code to asmbkg.F90 and asminc.F90 to allow writing of average fields (over namelist-specified window) into assimilation background file.

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   !!   'key_asminc' : Switch on the assimilation increment interface
21   !!----------------------------------------------------------------------
22   !!   asm_bkg_wri  : Write out the background state
23   !!   asm_trj_wri  : Write out the model state trajectory (used with 4D-Var)
24   !!----------------------------------------------------------------------
25   USE oce                ! Dynamics and active tracers defined in memory
26   USE sbc_oce            ! Ocean surface boundary conditions
27   USE zdf_oce            ! Vertical mixing variables
28   USE zdfddm             ! Double diffusion mixing parameterization
29   USE ldftra_oce         ! Lateral tracer mixing coefficient defined in memory
30   USE ldfslp             ! Slopes of neutral surfaces
31   USE tradmp             ! Tracer damping
32#if defined key_zdftke
33   USE zdftke             ! TKE vertical physics
34#endif
35   USE eosbn2             ! Equation of state (eos_bn2 routine)
36   USE zdfmxl             ! Mixed layer depth
37   USE dom_oce, ONLY :   ndastp
38   USE sol_oce, ONLY :   gcx   ! Solver variables defined in memory
39   USE in_out_manager     ! I/O manager
40   USE iom                ! I/O module
41   USE asmpar             ! Parameters for the assmilation interface
42   USE zdfmxl             ! mixed layer depth
43#if defined key_traldf_c2d
44   USE ldfeiv             ! eddy induced velocity coef.      (ldf_eiv routine)
45#endif
46#if defined key_lim2
47   USE ice_2
48#endif
49#if defined key_lim3
50   USE ice
51#endif
52   USE asminc, ONLY: ln_avgbkg
53   IMPLICIT NONE
54   PRIVATE
55   
56   PUBLIC   asm_bkg_wri   !: Write out the background state
57
58  !! * variables for calculating time means
59   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   tn_tavg  , sn_tavg 
60   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   un_tavg  , vn_tavg
61   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   avt_tavg
62#if defined key_zdfgls || key_zdftke
63   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   en_tavg
64#endif
65   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:)   ::   sshn_tavg
66   REAL(wp),SAVE :: gcx_tavg
67   REAL(wp),SAVE :: numtimes_tavg     ! No of times to average over
68
69   !!----------------------------------------------------------------------
70   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
71   !! $Id$
72   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
73   !!----------------------------------------------------------------------
74CONTAINS
75
76   SUBROUTINE asm_bkg_wri( kt )
77      !!-----------------------------------------------------------------------
78      !!                  ***  ROUTINE asm_bkg_wri ***
79      !!
80      !! ** Purpose : Write to file the background state for later use in the
81      !!              inner loop of data assimilation or for direct initialization
82      !!              in the outer loop.
83      !!
84      !! ** Method  : Write out the background state for use in the Jb term
85      !!              in the cost function and for use with direct initialization
86      !!              at analysis time.
87      !!-----------------------------------------------------------------------
88      INTEGER, INTENT( IN ) :: kt               ! Current time-step
89      !
90      CHARACTER (LEN=50) :: cl_asmbkg
91      CHARACTER (LEN=50) :: cl_asmdin
92      LOGICAL :: llok          ! Check if file exists
93      INTEGER :: inum          ! File unit number
94      REAL(wp) :: zdate        ! Date
95      !!-----------------------------------------------------------------------
96
97
98      ! If creating an averaged assim bkg, initialise on first timestep
99      IF ( ln_avgbkg .AND. kt == ( nn_it000 - 1) ) THEN
100         ! Allocate memory
101         ALLOCATE( tn_tavg(jpi,jpj,jpk), STAT=ierror )
102         IF( ierror > 0 ) THEN
103            CALL ctl_stop( 'asm_wri_bkg: unable to allocate tn_tavg' )   ;   RETURN
104         ENDIF
105         tn_tavg=0
106         ALLOCATE( sn_tavg(jpi,jpj,jpk), STAT=ierror )
107         IF( ierror > 0 ) THEN
108            CALL ctl_stop( 'asm_wri_bkg: unable to allocate sn_tavg' )   ;   RETURN
109         ENDIF
110         sn_tavg=0
111         ALLOCATE( un_tavg(jpi,jpj,jpk), STAT=ierror )
112         IF( ierror > 0 ) THEN
113            CALL ctl_stop( 'asm_wri_bkg: unable to allocate un_tavg' )   ;   RETURN
114         ENDIF
115         un_tavg=0
116         ALLOCATE( vn_tavg(jpi,jpj,jpk), STAT=ierror )
117         IF( ierror > 0 ) THEN
118            CALL ctl_stop( 'asm_wri_bkg: unable to allocate vn_tavg' )   ;   RETURN
119         ENDIF
120         vn_tavg=0
121         ALLOCATE( ssh_tavg(jpi,jpj), STAT=ierror )
122         IF( ierror > 0 ) THEN
123            CALL ctl_stop( 'asm_wri_bkg: unable to allocate ssh_tavg' )   ;   RETURN
124         ENDIF
125         ssh_tavg=0
126         ALLOCATE( en_tavg(jpi,jpj,jpk), STAT=ierror )
127         IF( ierror > 0 ) THEN
128            CALL ctl_stop( 'asm_wri_bkg: unable to allocate en_tavg' )   ;   RETURN
129         ENDIF
130         en_tavg=0
131         ALLOCATE( avt_tavg(jpi,jpj,jpk), STAT=ierror )
132         IF( ierror > 0 ) THEN
133            CALL ctl_stop( 'asm_wri_bkg: unable to allocate avt_tavg' )   ;   RETURN
134         ENDIF
135         avt_tavg=0
136         gcx_tavg=0   
137         
138         numtimes_tavg = REAL ( nitavgbkg_r -  nn_it000 + 1 )
139      ENDIF   
140
141      ! If creating an averaged assim bkg, sum the contribution every timestep
142      IF ( ln_avgbkg ) THEN
143         IF (lwp) THEN
144              WRITE(numout,*) 'asm_wri_bkg : Summing assim bkg fields at timestep ',kt
145              WRITE(numout,*) '~~~~~~~~~~~~ '
146         ENDIF
147
148         tn_tavg(:,:,:)        = tn_tavg(:,:,:) + tsn(:,:,:,jp_tem) / numtimes_tavg
149         sn_tavg(:,:,:)        = sn_tavg(:,:,:) + tsn(:,:,:,jp_sal) / numtimes_tavg
150         sshn_tavg(:,:)        = sshn_tavg(:,:) + sshn (:,:) / numtimes_tavg
151         un_tavg(:,:,:)        = un_tavg(:,:,:) + un(:,:,:) / numtimes_tavg
152         vn_tavg(:,:,:)        = vn_tavg(:,:,:) + vn(:,:,:) / numtimes_tavg
153         gcx_tavg              = gcx_tavg       + gcx / numtimes_tavg
154         avt_tavg(:,:,:)        = avt_tavg(:,:,:) + avt(:,:,:) / numtimes_tavg
155#if defined key_zdftke
156         en_tavg(:,:,:)       = en_tavg(:,:,:) + en(:,:,:) / numtimes_tavg
157#endif
158      ENDIF
159     
160
161      ! Write out background at time step nitbkg_r or nitavgbkg_r
162      IF( .NOT. ln_avgbkg .AND. (kt == nitbkg_r) ) .OR. &
163      &       ( ln_avgbkg .AND. (kt == nitavgbkg_r) ) THEN
164         !
165         WRITE(cl_asmbkg, FMT='(A,".nc")' ) TRIM( c_asmbkg )
166         cl_asmbkg = TRIM( cl_asmbkg )
167         INQUIRE( FILE = cl_asmbkg, EXIST = llok )
168         !
169         IF( .NOT. llok ) THEN
170            IF(lwp) WRITE(numout,*) ' Setting up assimilation background file '// TRIM( c_asmbkg )
171            !
172            !                                      ! Define the output file       
173            CALL iom_open( c_asmbkg, inum, ldwrt = .TRUE., kiolib = jprstlib)
174            !
175            IF( nitbkg_r == nit000 - 1 ) THEN      ! Treat special case when nitbkg = 0
176               zdate = REAL( ndastp )
177#if defined key_zdftke
178               ! lk_zdftke=T :   Read turbulent kinetic energy ( en )
179               IF(lwp) WRITE(numout,*) ' Reading TKE (en) from restart...'
180               CALL tke_rst( nit000, 'READ' )               ! lk_zdftke=T :   Read turbulent kinetic energy ( en )
181
182#endif
183            ELSE
184               zdate = REAL( ndastp )
185            ENDIF
186            !
187            ! Write the information
188            CALL iom_rstput( kt, nitbkg_r, inum, 'rdastp' , zdate   )
189           
190            IF ( ln_avgbkg ) THEN
191               CALL iom_rstput( kt, nitbkg_r, inum, 'un'     , un_tavg )
192               CALL iom_rstput( kt, nitbkg_r, inum, 'vn'     , vn_tavg )
193               CALL iom_rstput( kt, nitbkg_r, inum, 'tn'     , tn_tavg )
194               CALL iom_rstput( kt, nitbkg_r, inum, 'sn'     , sn_tavg )
195               CALL iom_rstput( kt, nitbkg_r, inum, 'sshn'   , sshn    )
196#if defined key_zdftke
197               CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en_tavg )
198#endif
199               CALL iom_rstput( kt, nitbkg_r, inum, 'gcx'    , gcx_tavg)
200               CALL iom_rstput( kt, nitbkg_r, inum, 'avt'    , avt_tavg)
201               !
202            ELSE
203               CALL iom_rstput( kt, nitbkg_r, inum, 'un'     , un                )
204               CALL iom_rstput( kt, nitbkg_r, inum, 'vn'     , vn                )
205               CALL iom_rstput( kt, nitbkg_r, inum, 'tn'     , tsn(:,:,:,jp_tem) )
206               CALL iom_rstput( kt, nitbkg_r, inum, 'sn'     , tsn(:,:,:,jp_sal) )
207               CALL iom_rstput( kt, nitbkg_r, inum, 'sshn'   , sshn              )
208#if defined key_zdftke
209               CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en                )
210#endif
211               CALL iom_rstput( kt, nitbkg_r, inum, 'gcx'    , gcx               )
212               CALL iom_rstput( kt, nitbkg_r, inum, 'avt'    , avt               )
213               !
214            ENDIF
215           
216            CALL iom_close( inum )
217
218         ENDIF
219         !
220      ENDIF
221
222      !                                !-------------------------------------------
223      IF( kt == nitdin_r ) THEN        ! Write out background at time step nitdin_r
224         !                             !-----------------------------------========
225         !
226         WRITE(cl_asmdin, FMT='(A,".nc")' ) TRIM( c_asmdin )
227         cl_asmdin = TRIM( cl_asmdin )
228         INQUIRE( FILE = cl_asmdin, EXIST = llok )
229         !
230         IF( .NOT. llok ) THEN
231            IF(lwp) WRITE(numout,*) ' Setting up assimilation background file '// TRIM( c_asmdin )
232            !
233            !                                      ! Define the output file       
234            CALL iom_open( c_asmdin, inum, ldwrt = .TRUE., kiolib = jprstlib)
235            !
236            IF( nitdin_r == nit000 - 1 ) THEN      ! Treat special case when nitbkg = 0
237
238               zdate = REAL( ndastp )
239            ELSE
240               zdate = REAL( ndastp )
241            ENDIF
242            !
243            !                                      ! Write the information
244            CALL iom_rstput( kt, nitdin_r, inum, 'rdastp' , zdate             )
245            CALL iom_rstput( kt, nitdin_r, inum, 'un'     , un                )
246            CALL iom_rstput( kt, nitdin_r, inum, 'vn'     , vn                )
247            CALL iom_rstput( kt, nitdin_r, inum, 'tn'     , tsn(:,:,:,jp_tem) )
248            CALL iom_rstput( kt, nitdin_r, inum, 'sn'     , tsn(:,:,:,jp_sal) )
249            CALL iom_rstput( kt, nitdin_r, inum, 'sshn'   , sshn              )
250#if defined key_lim2 || defined key_lim3
251            IF(( nn_ice == 2 ) .OR. ( nn_ice == 3 )) THEN
252          IF(ALLOCATED(frld)) THEN
253                  CALL iom_rstput( kt, nitdin_r, inum, 'iceconc', 1.0 - frld(:,:)   )
254               ELSE
255        CALL ctl_warn('Ice concentration not written to background as ice variable frld not allocated on this timestep')
256          ENDIF
257            ENDIF
258#endif
259            !
260            CALL iom_close( inum )
261         ENDIF
262         !
263      ENDIF
264      !                   
265   END SUBROUTINE asm_bkg_wri
266
267   !!======================================================================
268END MODULE asmbkg
Note: See TracBrowser for help on using the repository browser.