source: branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_ersem_hem08/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90 @ 9331

Last change on this file since 9331 was 9331, checked in by dford, 2 years ago

Add balancing code.

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