source: branches/UKMO/dev_r5518_v3.6_asm_nemovar_community/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90 @ 7731

Last change on this file since 7731 was 7731, checked in by dford, 3 years ago

Merge in revisions 6625:7726 of dev_r5518_v3.4_asm_nemovar_community, so this branch will be identical to revison 7726 of dev_r5518_v3.6_asm_nemovar_community.

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