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.
icbtrj.F90 in NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ICB – NEMO

source: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ICB/icbtrj.F90 @ 14038

Last change on this file since 14038 was 14038, checked in by laurent, 4 years ago

Catch up with trunk at rev r14037
  • Property svn:keywords set to Id
File size: 17.0 KB
Line 
1MODULE icbtrj
2   !!======================================================================
3   !!                       ***  MODULE  icbtrj  ***
4   !! Ocean physics:  trajectory I/O routines
5   !!======================================================================
6   !! History :  3.3  !  2010-01  (Martin&Adcroft) Original code
7   !!             -   !  2011-03  (Madec)          Part conversion to NEMO form
8   !!             -   !                            Removal of mapping from another grid
9   !!             -   !  2011-05  (Alderson)       New module to handle trajectory output
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   icb_trj_init  : initialise iceberg trajectory output files
14   !!   icb_trj_write :
15   !!   icb_trj_sync  :
16   !!   icb_trj_end   :
17   !!----------------------------------------------------------------------
18   USE par_oce        ! NEMO parameters
19   USE dom_oce        ! NEMO ocean domain
20   USE phycst         ! NEMO physical constants
21   USE icb_oce        ! define iceberg arrays
22   USE icbutl         ! iceberg utility routines
23   !
24   USE lib_mpp        ! NEMO MPI library, lk_mpp in particular
25   USE in_out_manager ! NEMO IO, numout in particular
26   USE ioipsl  , ONLY : ju2ymds    ! for calendar
27   USE netcdf
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   icb_trj_init    ! routine called in icbini.F90 module
33   PUBLIC   icb_trj_write   ! routine called in icbstp.F90 module
34   PUBLIC   icb_trj_sync    ! routine called in icbstp.F90 module
35   PUBLIC   icb_trj_end     ! routine called in icbstp.F90 module
36
37   INTEGER ::   num_traj = 0
38   INTEGER ::   n_dim, m_dim
39   INTEGER ::   ntrajid
40   INTEGER ::   numberid, nstepid, nscaling_id
41   INTEGER ::   nlonid, nlatid, nxid, nyid, nuvelid, nvvelid, nmassid
42   INTEGER ::   nssuid, nssvid, nuaid, nvaid, nuiid, nviid
43   INTEGER ::   nsshxid, nsshyid, nsstid, ncntid, nthkid
44   INTEGER ::   nthicknessid, nwidthid, nlengthid
45   INTEGER ::   nyearid, ndayid
46   INTEGER ::   nmass_of_bits_id, nheat_density_id
47
48   !!----------------------------------------------------------------------
49   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
50   !! $Id$
51   !! Software governed by the CeCILL license (see ./LICENSE)
52   !!----------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE icb_trj_init( ktend )
56      !!----------------------------------------------------------------------
57      !!                  ***  ROUTINE icb_trj_init  ***
58      !!
59      !! ** Purpose :   initialise iceberg trajectory output files
60      !!----------------------------------------------------------------------
61      INTEGER, INTENT(in) ::   ktend   ! time step index
62      !
63      INTEGER                ::   iret, iyear, imonth, iday
64      INTEGER                ::   idg  ! number of digits
65      REAL(wp)               ::   zfjulday, zsec
66      CHARACTER(len=80)      ::   cl_filename
67      CHARACTER(LEN=12)      ::   clfmt            ! writing format
68      CHARACTER(LEN=8 )      ::   cldate_ini, cldate_end
69      TYPE(iceberg), POINTER ::   this
70      TYPE(point)  , POINTER ::   pt
71      !!----------------------------------------------------------------------
72
73      ! compute initial time step date
74      CALL ju2ymds( fjulday, iyear, imonth, iday, zsec )
75      WRITE(cldate_ini, '(i4.4,2i2.2)') iyear, imonth, iday
76
77      ! compute end time step date
78      zfjulday = fjulday + rn_Dt / rday * REAL( nitend - nit000 + 1 , wp)
79      IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday )   zfjulday = REAL(NINT(zfjulday),wp)   ! avoid truncation error
80      CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec )
81      WRITE(cldate_end, '(i4.4,2i2.2)') iyear, imonth, iday
82
83      ! define trajectory output name
84      cl_filename = 'trajectory_icebergs_'//cldate_ini//'-'//cldate_end
85      IF ( lk_mpp ) THEN
86         idg = MAX( INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1, 4 )          ! how many digits to we need to write? min=4, max=9
87         WRITE(clfmt, "('(a,a,i', i1, '.', i1, ',a)')") idg, idg          ! '(a,a,ix.x,a)'
88         WRITE(cl_filename,  clfmt) TRIM(cl_filename), '_', narea-1, '.nc'
89      ELSE
90         WRITE(cl_filename,'(a,a)') TRIM(cl_filename),               '.nc'
91      ENDIF
92      IF( lwp .AND. nn_verbose_level >= 0 )   WRITE(numout,'(2a)') 'icebergs, icb_trj_init: creating ',TRIM(cl_filename)
93
94      iret = NF90_CREATE( TRIM(cl_filename), NF90_CLOBBER, ntrajid )
95      IF (iret .NE. NF90_NOERR)   CALL ctl_stop('icebergs, icb_trj_init: nf_create failed')
96
97      ! Dimensions
98      iret = NF90_DEF_DIM( ntrajid, 'n', NF90_UNLIMITED, n_dim )
99      IF ( iret /= NF90_NOERR )   CALL ctl_stop('icebergs, icb_trj_init: nf_def_dim n failed')
100      iret = NF90_DEF_DIM( ntrajid, 'k', nkounts, m_dim )
101      IF ( iret /= NF90_NOERR )   CALL ctl_stop('icebergs, icb_trj_init: nf_def_dim k failed')
102
103      ! Variables
104      iret = NF90_DEF_VAR( ntrajid, 'iceberg_number', NF90_INT   , (/m_dim,n_dim/), numberid         )
105      iret = NF90_DEF_VAR( ntrajid, 'timestep'      , NF90_INT   , n_dim          , nstepid          )
106      iret = NF90_DEF_VAR( ntrajid, 'mass_scaling'  , NF90_DOUBLE, n_dim          , nscaling_id      )
107      iret = NF90_DEF_VAR( ntrajid, 'lon'           , NF90_DOUBLE, n_dim          , nlonid           )
108      iret = NF90_DEF_VAR( ntrajid, 'lat'           , NF90_DOUBLE, n_dim          , nlatid           )
109      iret = NF90_DEF_VAR( ntrajid, 'xi'            , NF90_DOUBLE, n_dim          , nxid             )
110      iret = NF90_DEF_VAR( ntrajid, 'yj'            , NF90_DOUBLE, n_dim          , nyid             )
111      iret = NF90_DEF_VAR( ntrajid, 'uvel'          , NF90_DOUBLE, n_dim          , nuvelid          )
112      iret = NF90_DEF_VAR( ntrajid, 'vvel'          , NF90_DOUBLE, n_dim          , nvvelid          )
113      iret = NF90_DEF_VAR( ntrajid, 'ssu'           , NF90_DOUBLE, n_dim          , nssuid           )
114      iret = NF90_DEF_VAR( ntrajid, 'ssv'           , NF90_DOUBLE, n_dim          , nssvid           )
115      iret = NF90_DEF_VAR( ntrajid, 'uta'           , NF90_DOUBLE, n_dim          , nuaid            )
116      iret = NF90_DEF_VAR( ntrajid, 'vta'           , NF90_DOUBLE, n_dim          , nvaid            )
117      iret = NF90_DEF_VAR( ntrajid, 'uti'           , NF90_DOUBLE, n_dim          , nuiid            )
118      iret = NF90_DEF_VAR( ntrajid, 'vti'           , NF90_DOUBLE, n_dim          , nviid            )
119      iret = NF90_DEF_VAR( ntrajid, 'ssh_x'         , NF90_DOUBLE, n_dim          , nsshxid          )
120      iret = NF90_DEF_VAR( ntrajid, 'ssh_y'         , NF90_DOUBLE, n_dim          , nsshyid          )
121      iret = NF90_DEF_VAR( ntrajid, 'sst'           , NF90_DOUBLE, n_dim          , nsstid           )
122      iret = NF90_DEF_VAR( ntrajid, 'icnt'          , NF90_DOUBLE, n_dim          , ncntid           )
123      iret = NF90_DEF_VAR( ntrajid, 'ithk'          , NF90_DOUBLE, n_dim          , nthkid           )
124      iret = NF90_DEF_VAR( ntrajid, 'mass'          , NF90_DOUBLE, n_dim          , nmassid          )
125      iret = NF90_DEF_VAR( ntrajid, 'thickness'     , NF90_DOUBLE, n_dim          , nthicknessid     )
126      iret = NF90_DEF_VAR( ntrajid, 'width'         , NF90_DOUBLE, n_dim          , nwidthid         )
127      iret = NF90_DEF_VAR( ntrajid, 'length'        , NF90_DOUBLE, n_dim          , nlengthid        )
128      iret = NF90_DEF_VAR( ntrajid, 'year'          , NF90_INT   , n_dim          , nyearid          )
129      iret = NF90_DEF_VAR( ntrajid, 'day'           , NF90_DOUBLE, n_dim          , ndayid           )
130      iret = NF90_DEF_VAR( ntrajid, 'mass_of_bits'  , NF90_DOUBLE, n_dim          , nmass_of_bits_id )
131      iret = NF90_DEF_VAR( ntrajid, 'heat_density'  , NF90_DOUBLE, n_dim          , nheat_density_id )
132
133      ! Attributes
134      iret = NF90_PUT_ATT( ntrajid, numberid        , 'long_name', 'iceberg number on this processor' )
135      iret = NF90_PUT_ATT( ntrajid, numberid        , 'units'    , 'count' )
136      iret = NF90_PUT_ATT( ntrajid, nstepid         , 'long_name', 'timestep number kt' )
137      iret = NF90_PUT_ATT( ntrajid, nstepid         , 'units'    , 'count' )
138      iret = NF90_PUT_ATT( ntrajid, nlonid          , 'long_name', 'longitude' )
139      iret = NF90_PUT_ATT( ntrajid, nlonid          , 'units'    , 'degrees_E')
140      iret = NF90_PUT_ATT( ntrajid, nlatid          , 'long_name', 'latitude' )
141      iret = NF90_PUT_ATT( ntrajid, nlatid          , 'units'    , 'degrees_N' )
142      iret = NF90_PUT_ATT( ntrajid, nxid            , 'long_name', 'x grid box position' )
143      iret = NF90_PUT_ATT( ntrajid, nxid            , 'units'    , 'fractional' )
144      iret = NF90_PUT_ATT( ntrajid, nyid            , 'long_name', 'y grid box position' )
145      iret = NF90_PUT_ATT( ntrajid, nyid            , 'units'    , 'fractional' )
146      iret = NF90_PUT_ATT( ntrajid, nuvelid         , 'long_name', 'zonal velocity' )
147      iret = NF90_PUT_ATT( ntrajid, nuvelid         , 'units'    , 'm/s' )
148      iret = NF90_PUT_ATT( ntrajid, nvvelid         , 'long_name', 'meridional velocity' )
149      iret = NF90_PUT_ATT( ntrajid, nvvelid         , 'units'    , 'm/s' )
150      iret = NF90_PUT_ATT( ntrajid, nssuid          , 'long_name', 'ocean u component' )
151      iret = NF90_PUT_ATT( ntrajid, nssuid          , 'units'    , 'm/s' )
152      iret = NF90_PUT_ATT( ntrajid, nssvid          , 'long_name', 'ocean v component' )
153      iret = NF90_PUT_ATT( ntrajid, nssvid          , 'units'    , 'm/s' )
154      iret = NF90_PUT_ATT( ntrajid, nuaid           , 'long_name', 'atmosphere u component' )
155      iret = NF90_PUT_ATT( ntrajid, nuaid           , 'units'    , 'm/s' )
156      iret = NF90_PUT_ATT( ntrajid, nvaid           , 'long_name', 'atmosphere v component' )
157      iret = NF90_PUT_ATT( ntrajid, nvaid           , 'units'    , 'm/s' )
158      iret = NF90_PUT_ATT( ntrajid, nuiid           , 'long_name', 'sea ice u component' )
159      iret = NF90_PUT_ATT( ntrajid, nuiid           , 'units'    , 'm/s' )
160      iret = NF90_PUT_ATT( ntrajid, nviid           , 'long_name', 'sea ice v component' )
161      iret = NF90_PUT_ATT( ntrajid, nviid           , 'units'    , 'm/s' )
162      iret = NF90_PUT_ATT( ntrajid, nsshxid         , 'long_name', 'sea surface height gradient from x points' )
163      iret = NF90_PUT_ATT( ntrajid, nsshxid         , 'units'    , 'm/m' )
164      iret = NF90_PUT_ATT( ntrajid, nsshyid         , 'long_name', 'sea surface height gradient from y points' )
165      iret = NF90_PUT_ATT( ntrajid, nsshyid         , 'units'    , 'm/m' )
166      iret = NF90_PUT_ATT( ntrajid, nsstid          , 'long_name', 'sea surface temperature' )
167      iret = NF90_PUT_ATT( ntrajid, nsstid          , 'units'    , 'degC')
168      iret = NF90_PUT_ATT( ntrajid, ncntid          , 'long_name', 'sea ice concentration' )
169      iret = NF90_PUT_ATT( ntrajid, ncntid          , 'units'    , 'degC')
170      iret = NF90_PUT_ATT( ntrajid, nthkid          , 'long_name', 'sea ice thickness' )
171      iret = NF90_PUT_ATT( ntrajid, nthkid          , 'units'    , 'm'   )
172      iret = NF90_PUT_ATT( ntrajid, nmassid         , 'long_name', 'mass')
173      iret = NF90_PUT_ATT( ntrajid, nmassid         , 'units'    , 'kg'  )
174      iret = NF90_PUT_ATT( ntrajid, nthicknessid    , 'long_name', 'thickness' )
175      iret = NF90_PUT_ATT( ntrajid, nthicknessid    , 'units'    , 'm'   )
176      iret = NF90_PUT_ATT( ntrajid, nwidthid        , 'long_name', 'width' )
177      iret = NF90_PUT_ATT( ntrajid, nwidthid        , 'units'    , 'm'    )
178      iret = NF90_PUT_ATT( ntrajid, nlengthid       , 'long_name', 'length' )
179      iret = NF90_PUT_ATT( ntrajid, nlengthid       , 'units'    , 'm'    )
180      iret = NF90_PUT_ATT( ntrajid, nyearid         , 'long_name', 'calendar year' )
181      iret = NF90_PUT_ATT( ntrajid, nyearid         , 'units'    , 'years' )
182      iret = NF90_PUT_ATT( ntrajid, ndayid          , 'long_name', 'day of year' )
183      iret = NF90_PUT_ATT( ntrajid, ndayid          , 'units'    , 'days' )
184      iret = NF90_PUT_ATT( ntrajid, nscaling_id     , 'long_name', 'scaling factor for mass of berg' )
185      iret = NF90_PUT_ATT( ntrajid, nscaling_id     , 'units'    , 'none' )
186      iret = NF90_PUT_ATT( ntrajid, nmass_of_bits_id, 'long_name', 'mass of bergy bits' )
187      iret = NF90_PUT_ATT( ntrajid, nmass_of_bits_id, 'units'    , 'kg'   )
188      iret = NF90_PUT_ATT( ntrajid, nheat_density_id, 'long_name', 'heat density' )
189      iret = NF90_PUT_ATT( ntrajid, nheat_density_id, 'units'    , 'J/kg' )
190      !
191      ! End define mode
192      iret = NF90_ENDDEF( ntrajid )
193      !
194   END SUBROUTINE icb_trj_init
195
196
197   SUBROUTINE icb_trj_write( kt )
198      !!----------------------------------------------------------------------
199      !!                  ***  ROUTINE icb_trj_write  ***
200      !!
201      !! ** Purpose :   write out iceberg trajectories
202      !!
203      !! ** Method  : - for the moment write out each snapshot of positions later
204      !!                can rewrite so that it is buffered and written out more efficiently
205      !!----------------------------------------------------------------------
206      INTEGER, INTENT( in ) ::   kt   ! time-step index
207      !
208      INTEGER                ::   iret, jn
209      CHARACTER(len=80)      ::   cl_filename
210      TYPE(iceberg), POINTER ::   this
211      TYPE(point  ), POINTER ::   pt
212      !!----------------------------------------------------------------------
213
214      ! Write variables
215      ! sga - just write out the current point of the trajectory
216
217      this => first_berg
218      jn = num_traj
219      DO WHILE( ASSOCIATED(this) )
220         pt => this%current_point
221         jn = jn + 1
222         !
223         iret = NF90_PUT_VAR( ntrajid, numberid        , this%number      , (/1,jn/) , (/nkounts,1/) )
224         iret = NF90_PUT_VAR( ntrajid, nstepid         , kt               , (/ jn /) )
225         iret = NF90_PUT_VAR( ntrajid, nscaling_id     , this%mass_scaling, (/ jn /) )
226         !
227         iret = NF90_PUT_VAR( ntrajid, nlonid          , pt%lon           , (/ jn /) )
228         iret = NF90_PUT_VAR( ntrajid, nlatid          , pt%lat           , (/ jn /) )
229         iret = NF90_PUT_VAR( ntrajid, nxid            , pt%xi            , (/ jn /) )
230         iret = NF90_PUT_VAR( ntrajid, nyid            , pt%yj            , (/ jn /) )
231         iret = NF90_PUT_VAR( ntrajid, nuvelid         , pt%uvel          , (/ jn /) )
232         iret = NF90_PUT_VAR( ntrajid, nvvelid         , pt%vvel          , (/ jn /) )
233         iret = NF90_PUT_VAR( ntrajid, nssuid          , pt%ssu           , (/ jn /) )
234         iret = NF90_PUT_VAR( ntrajid, nssvid          , pt%ssv           , (/ jn /) )
235         iret = NF90_PUT_VAR( ntrajid, nuaid           , pt%ua            , (/ jn /) )
236         iret = NF90_PUT_VAR( ntrajid, nvaid           , pt%va            , (/ jn /) )
237         iret = NF90_PUT_VAR( ntrajid, nuiid           , pt%ui            , (/ jn /) )
238         iret = NF90_PUT_VAR( ntrajid, nviid           , pt%vi            , (/ jn /) )
239         iret = NF90_PUT_VAR( ntrajid, nsshxid         , pt%ssh_x         , (/ jn /) )
240         iret = NF90_PUT_VAR( ntrajid, nsshyid         , pt%ssh_y         , (/ jn /) )
241         iret = NF90_PUT_VAR( ntrajid, nsstid          , pt%sst           , (/ jn /) )
242         iret = NF90_PUT_VAR( ntrajid, ncntid          , pt%cn            , (/ jn /) )
243         iret = NF90_PUT_VAR( ntrajid, nthkid          , pt%hi            , (/ jn /) )
244         iret = NF90_PUT_VAR( ntrajid, nmassid         , pt%mass          , (/ jn /) )
245         iret = NF90_PUT_VAR( ntrajid, nthicknessid    , pt%thickness     , (/ jn /) )
246         iret = NF90_PUT_VAR( ntrajid, nwidthid        , pt%width         , (/ jn /) )
247         iret = NF90_PUT_VAR( ntrajid, nlengthid       , pt%length        , (/ jn /) )
248         iret = NF90_PUT_VAR( ntrajid, nyearid         , pt%year          , (/ jn /) )
249         iret = NF90_PUT_VAR( ntrajid, ndayid          , pt%day           , (/ jn /) )
250         iret = NF90_PUT_VAR( ntrajid, nmass_of_bits_id, pt%mass_of_bits  , (/ jn /) )
251         iret = NF90_PUT_VAR( ntrajid, nheat_density_id, pt%heat_density  , (/ jn /) )
252         !
253         this => this%next
254      END DO
255      IF( lwp .AND. nn_verbose_level > 0 )   WRITE(numout,*) 'trajectory write to frame ', jn
256      num_traj = jn
257      !
258   END SUBROUTINE icb_trj_write
259
260
261   SUBROUTINE icb_trj_sync()
262      !!----------------------------------------------------------------------
263      !!                  ***  ROUTINE icb_trj_sync  ***
264      !!
265      !! ** Purpose :   
266      !!----------------------------------------------------------------------
267      INTEGER ::   iret
268      !!----------------------------------------------------------------------
269      ! flush to file
270      iret = NF90_SYNC( ntrajid )
271      IF ( iret /= NF90_NOERR )   CALL ctl_stop( 'icebergs, icb_trj_sync: nf_sync failed' )
272      !
273   END SUBROUTINE icb_trj_sync
274
275
276   SUBROUTINE icb_trj_end()
277      !!----------------------------------------------------------------------
278      INTEGER ::   iret
279      !!----------------------------------------------------------------------
280      ! Finish up
281      iret = NF90_CLOSE( ntrajid )
282      IF ( iret /= NF90_NOERR )   CALL ctl_stop( 'icebergs, icb_trj_end: nf_close failed' )
283      !
284   END SUBROUTINE icb_trj_end
285
286   !!======================================================================
287END MODULE icbtrj
Note: See TracBrowser for help on using the repository browser.