source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90 @ 10574

Last change on this file since 10574 was 10574, checked in by dford, 19 months ago

Merge in functionality added to GO6 at r10149.

File size: 12.8 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#if defined key_top
54   USE asmbgc, ONLY: asm_bgc_bkg_wri
55#endif
56   IMPLICIT NONE
57   PRIVATE
58   
59   PUBLIC   asm_bkg_wri   !: Write out the background state
60
61  !! * variables for calculating time means
62   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   tn_tavg  , sn_tavg 
63   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   un_tavg  , vn_tavg
64   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   avt_tavg
65#if defined key_zdfgls || key_zdftke
66   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   en_tavg
67#endif
68   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:)   ::   sshn_tavg
69   REAL(wp),SAVE :: numtimes_tavg     ! No of times to average over
70
71   !!----------------------------------------------------------------------
72   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
73   !! $Id$
74   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
75   !!----------------------------------------------------------------------
76CONTAINS
77
78   SUBROUTINE asm_bkg_wri( kt )
79      !!-----------------------------------------------------------------------
80      !!                  ***  ROUTINE asm_bkg_wri ***
81      !!
82      !! ** Purpose : Write to file the background state for later use in the
83      !!              inner loop of data assimilation or for direct initialization
84      !!              in the outer loop.
85      !!
86      !! ** Method  : Write out the background state for use in the Jb term
87      !!              in the cost function and for use with direct initialization
88      !!              at analysis time.
89      !!-----------------------------------------------------------------------
90      INTEGER, INTENT( IN ) :: kt               ! Current time-step
91      !
92      CHARACTER (LEN=50) :: cl_asmbkg
93      CHARACTER (LEN=50) :: cl_asmdin
94      LOGICAL :: llok          ! Check if file exists
95      INTEGER :: inum          ! File unit number
96      REAL(wp) :: zdate        ! Date
97      INTEGER :: ierror
98      !!-----------------------------------------------------------------------
99
100      ! If creating an averaged assim bkg, initialise on first timestep
101      IF ( ln_avgbkg .AND. kt == ( nn_it000 - 1) ) THEN
102         ! Allocate memory
103         ALLOCATE( tn_tavg(jpi,jpj,jpk), STAT=ierror )
104         IF( ierror > 0 ) THEN
105            CALL ctl_stop( 'asm_wri_bkg: unable to allocate tn_tavg' )   ;   RETURN
106         ENDIF
107         tn_tavg(:,:,:)=0
108         ALLOCATE( sn_tavg(jpi,jpj,jpk), STAT=ierror )
109         IF( ierror > 0 ) THEN
110            CALL ctl_stop( 'asm_wri_bkg: unable to allocate sn_tavg' )   ;   RETURN
111         ENDIF
112         sn_tavg(:,:,:)=0
113         ALLOCATE( un_tavg(jpi,jpj,jpk), STAT=ierror )
114         IF( ierror > 0 ) THEN
115            CALL ctl_stop( 'asm_wri_bkg: unable to allocate un_tavg' )   ;   RETURN
116         ENDIF
117         un_tavg(:,:,:)=0
118         ALLOCATE( vn_tavg(jpi,jpj,jpk), STAT=ierror )
119         IF( ierror > 0 ) THEN
120            CALL ctl_stop( 'asm_wri_bkg: unable to allocate vn_tavg' )   ;   RETURN
121         ENDIF
122         vn_tavg(:,:,:)=0
123         ALLOCATE( sshn_tavg(jpi,jpj), STAT=ierror )
124         IF( ierror > 0 ) THEN
125            CALL ctl_stop( 'asm_wri_bkg: unable to allocate sshn_tavg' )   ;   RETURN
126         ENDIF
127         sshn_tavg(:,:)=0
128#if defined key_zdftke
129         ALLOCATE( en_tavg(jpi,jpj,jpk), STAT=ierror )
130         IF( ierror > 0 ) THEN
131            CALL ctl_stop( 'asm_wri_bkg: unable to allocate en_tavg' )   ;   RETURN
132         ENDIF
133         en_tavg(:,:,:)=0
134#endif
135         ALLOCATE( avt_tavg(jpi,jpj,jpk), STAT=ierror )
136         IF( ierror > 0 ) THEN
137            CALL ctl_stop( 'asm_wri_bkg: unable to allocate avt_tavg' )   ;   RETURN
138         ENDIF
139         avt_tavg(:,:,:)=0
140         
141         numtimes_tavg = REAL ( nitavgbkg_r -  nn_it000 + 1 )
142      ENDIF   
143
144      ! If creating an averaged assim bkg, sum the contribution every timestep
145      IF ( ln_avgbkg ) THEN
146         IF (lwp) THEN
147              WRITE(numout,*) 'asm_wri_bkg : Summing assim bkg fields at timestep ',kt
148              WRITE(numout,*) '~~~~~~~~~~~~ '
149         ENDIF
150
151         tn_tavg(:,:,:)        = tn_tavg(:,:,:) + tsn(:,:,:,jp_tem) / numtimes_tavg
152         sn_tavg(:,:,:)        = sn_tavg(:,:,:) + tsn(:,:,:,jp_sal) / numtimes_tavg
153         sshn_tavg(:,:)        = sshn_tavg(:,:) + sshn (:,:) / numtimes_tavg
154         un_tavg(:,:,:)        = un_tavg(:,:,:) + un(:,:,:) / numtimes_tavg
155         vn_tavg(:,:,:)        = vn_tavg(:,:,:) + vn(:,:,:) / numtimes_tavg
156         avt_tavg(:,:,:)        = avt_tavg(:,:,:) + avt(:,:,:) / numtimes_tavg
157#if defined key_zdftke
158         en_tavg(:,:,:)       = en_tavg(:,:,:) + en(:,:,:) / numtimes_tavg
159#endif
160      ENDIF
161     
162
163      ! Write out background at time step nitbkg_r or nitavgbkg_r
164      IF ( ( .NOT. ln_avgbkg .AND. (kt == nitbkg_r) ) .OR. &
165      &          ( ln_avgbkg .AND. (kt == nitavgbkg_r) ) ) THEN
166         !
167         WRITE(cl_asmbkg, FMT='(A,".nc")' ) TRIM( c_asmbkg )
168         cl_asmbkg = TRIM( cl_asmbkg )
169         INQUIRE( FILE = cl_asmbkg, EXIST = llok )
170         !
171         IF( .NOT. llok ) THEN
172            IF(lwp) WRITE(numout,*) ' Setting up assimilation background file '// TRIM( c_asmbkg )
173            !
174            !                                      ! Define the output file       
175            CALL iom_open( c_asmbkg, inum, ldwrt = .TRUE., kiolib = jprstlib)
176            !
177            !
178            ! Write the information
179            IF ( ln_avgbkg ) THEN
180               IF( nitavgbkg_r == nit000 - 1 ) THEN      ! Treat special case when nitavgbkg = 0
181                  zdate = REAL( ndastp )
182#if defined key_zdftke
183                  ! lk_zdftke=T :   Read turbulent kinetic energy ( en )
184                  IF(lwp) WRITE(numout,*) ' Reading TKE (en) from restart...'
185                  CALL tke_rst( nit000, 'READ' )               ! lk_zdftke=T :   Read turbulent kinetic energy ( en )
186
187#endif
188               ELSE
189                  zdate = REAL( ndastp )
190               ENDIF
191               CALL iom_rstput( kt, nitavgbkg_r, inum, 'rdastp' , zdate   )
192               CALL iom_rstput( kt, nitavgbkg_r, inum, 'un'     , un_tavg )
193               CALL iom_rstput( kt, nitavgbkg_r, inum, 'vn'     , vn_tavg )
194               CALL iom_rstput( kt, nitavgbkg_r, inum, 'tn'     , tn_tavg )
195               CALL iom_rstput( kt, nitavgbkg_r, inum, 'sn'     , sn_tavg )
196               CALL iom_rstput( kt, nitavgbkg_r, inum, 'sshn'   , sshn_tavg)
197#if defined key_zdftke
198               CALL iom_rstput( kt, nitavgbkg_r, inum, 'en'     , en_tavg )
199#endif
200               CALL iom_rstput( kt, nitavgbkg_r, inum, 'avt'    , avt_tavg)
201               !
202            ELSE
203               IF( nitbkg_r == nit000 - 1 ) THEN      ! Treat special case when nitbkg = 0
204                  zdate = REAL( ndastp )
205#if defined key_zdftke
206                  ! lk_zdftke=T :   Read turbulent kinetic energy ( en )
207                  IF(lwp) WRITE(numout,*) ' Reading TKE (en) from restart...'
208                  CALL tke_rst( nit000, 'READ' )               ! lk_zdftke=T :   Read turbulent kinetic energy ( en )
209
210#endif
211               ELSE
212                  zdate = REAL( ndastp )
213               ENDIF
214               CALL iom_rstput( kt, nitbkg_r, inum, 'rdastp' , zdate   )
215               CALL iom_rstput( kt, nitbkg_r, inum, 'un'     , un                )
216               CALL iom_rstput( kt, nitbkg_r, inum, 'vn'     , vn                )
217               CALL iom_rstput( kt, nitbkg_r, inum, 'tn'     , tsn(:,:,:,jp_tem) )
218               CALL iom_rstput( kt, nitbkg_r, inum, 'sn'     , tsn(:,:,:,jp_sal) )
219               CALL iom_rstput( kt, nitbkg_r, inum, 'sshn'   , sshn              )
220#if defined key_zdftke
221               CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en                )
222#endif
223               CALL iom_rstput( kt, nitbkg_r, inum, 'avt'    , avt               )
224               !
225            ENDIF
226           
227#if defined key_top
228            CALL asm_bgc_bkg_wri( kt, inum )
229#endif
230            CALL iom_close( inum )
231
232         ENDIF
233         !
234      ENDIF
235
236      !                                !-------------------------------------------
237      IF( kt == nitdin_r ) THEN        ! Write out background at time step nitdin_r
238         !                             !-----------------------------------========
239         !
240         WRITE(cl_asmdin, FMT='(A,".nc")' ) TRIM( c_asmdin )
241         cl_asmdin = TRIM( cl_asmdin )
242         INQUIRE( FILE = cl_asmdin, EXIST = llok )
243         !
244         IF( .NOT. llok ) THEN
245            IF(lwp) WRITE(numout,*) ' Setting up assimilation background file '// TRIM( c_asmdin )
246            !
247            !                                      ! Define the output file       
248            CALL iom_open( c_asmdin, inum, ldwrt = .TRUE., kiolib = jprstlib)
249            !
250            IF( nitdin_r == nit000 - 1 ) THEN      ! Treat special case when nitbkg = 0
251
252               zdate = REAL( ndastp )
253            ELSE
254               zdate = REAL( ndastp )
255            ENDIF
256            !
257            !                                      ! Write the information
258            CALL iom_rstput( kt, nitdin_r, inum, 'rdastp' , zdate             )
259            CALL iom_rstput( kt, nitdin_r, inum, 'un'     , un                )
260            CALL iom_rstput( kt, nitdin_r, inum, 'vn'     , vn                )
261            CALL iom_rstput( kt, nitdin_r, inum, 'tn'     , tsn(:,:,:,jp_tem) )
262            CALL iom_rstput( kt, nitdin_r, inum, 'sn'     , tsn(:,:,:,jp_sal) )
263            CALL iom_rstput( kt, nitdin_r, inum, 'sshn'   , sshn              )
264#if defined key_lim2 || defined key_lim3
265            IF(( nn_ice == 2 ) .OR. ( nn_ice == 3 )) THEN
266          IF(ALLOCATED(frld)) THEN
267                  CALL iom_rstput( kt, nitdin_r, inum, 'iceconc', 1.0 - frld(:,:)   )
268               ELSE
269        CALL ctl_warn('Ice concentration not written to background as ice variable frld not allocated on this timestep')
270          ENDIF
271            ENDIF
272#endif
273            !
274            CALL iom_close( inum )
275         ENDIF
276         !
277      ENDIF
278      !                   
279   END SUBROUTINE asm_bkg_wri
280
281   !!======================================================================
282END MODULE asmbkg
Note: See TracBrowser for help on using the repository browser.