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.
trcini_medusa.F90 in branches/UKMO/dev_r5518_GO6_package_FOAMv14_Hs_CO2/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_GO6_package_FOAMv14_Hs_CO2/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcini_medusa.F90 @ 15463

Last change on this file since 15463 was 15463, checked in by dford, 8 months ago

CO2 flux based on input significant wave height.

File size: 24.1 KB
Line 
1MODULE trcini_medusa
2   !!======================================================================
3   !!                         ***  MODULE trcini_medusa  ***
4   !! TOP :   initialisation of the MEDUSA tracers
5   !!======================================================================
6   !! History :   2.0  !  2007-12  (C. Ethe, G. Madec) Original code
7   !!              -   !  2008-08  (K. Popova) adaptation for MEDUSA
8   !!              -   !  2008-11  (A. Yool) continuing adaptation for MEDUSA
9   !!              -   !  2010-03  (A. Yool) updated for branch inclusion
10   !!              -   !  2011-04  (A. Yool) updated for ROAM project
11   !!----------------------------------------------------------------------
12#if defined key_medusa
13   !!----------------------------------------------------------------------
14   !!   'key_medusa'                                         MEDUSA tracers
15   !!----------------------------------------------------------------------
16   !! trc_ini_medusa   : MEDUSA model initialisation
17   !!----------------------------------------------------------------------
18   USE par_trc         ! TOP parameters
19   USE oce_trc
20   USE trc
21   USE in_out_manager
22   !! AXY (04/11/13): add this in for initialisation stuff
23   USE iom
24   USE par_medusa
25   !! AXY (13/01/12): add this in for sediment variables
26   USE sms_medusa
27   !! AXY (04/11/13): add this in for initialisation stuff
28   USE trcsed_medusa
29   USE sbc_oce, ONLY: lk_oasis
30   USE oce,     ONLY: CO2Flux_out_cpl, DMS_out_cpl, chloro_out_cpl  !! Coupling variable
31
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   trc_ini_medusa   ! called by trcini.F90 module
37
38   !! AXY (25/02/10)
39   LOGICAL, PUBLIC ::                  &
40      bocalccd = .TRUE.
41   !! JPALM (14/09/15)
42   LOGICAL, PUBLIC ::                  &
43      ln_ccd = .TRUE.
44
45   INTEGER ::                          &
46      numccd, numswh
47
48   !! AXY (25/02/10)
49   INTEGER ::                          &
50      numriv
51
52   !!----------------------------------------------------------------------
53   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
54   !! $Id$
55   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
56   !!----------------------------------------------------------------------
57
58CONTAINS
59
60   SUBROUTINE trc_ini_medusa
61      !!----------------------------------------------------------------------
62      !!                     ***  trc_ini_medusa  *** 
63      !!
64      !! ** Purpose :   initialization for MEDUSA model
65      !!
66      !! ** Method  : - Read the namcfc namelist and check the parameter values
67      !!----------------------------------------------------------------------
68      !!----------------------------------------------------------------------
69
70      !! vertical array index
71      INTEGER  ::    jk, ierr
72      !! AXY (19/07/12): added jk2 to set up friver_dep array
73      INTEGER            :: jk2
74      !! AXY (19/07/12): added tfthk to set up friver_dep array
75      REAL(wp)           :: fthk, tfthk
76      !! AXY (04/11/13): add in temporary variables for checks
77      REAL(wp)           :: fq0, fq1, fq2
78
79      IF(lwp) WRITE(numout,*)
80      IF(lwp) WRITE(numout,*) ' trc_ini_medusa: initialisation of MEDUSA model'
81      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
82# if defined key_debug_medusa
83            CALL flush(numout)
84# endif
85
86                                           ! Allocate MEDUSA arrays
87      ierr =         sms_medusa_alloc()
88# if defined key_debug_medusa
89            IF (lwp) write (numout,*) '------------------------------'
90            IF (lwp) write (numout,*) 'Jpalm - debug'
91            IF (lwp) write (numout,*) 'in trc_ini_medusa, just after array allocate'
92            IF (lwp) write (numout,*) ' '
93            CALL flush(numout)
94# endif
95
96!!
97!! AXY (19/07/12): setup array to control distribution of river nutrients
98      friver_dep(:,:) = 0.
99      DO jk = 1,jpk
100         tfthk = 0.
101         DO jk2 = 1,jriver_dep
102            fthk  = e3t_1d(jk2)
103            if (jk2 .le. jk) then
104               tfthk = tfthk + fthk
105               friver_dep(jk2,jk) = fthk
106            endif
107         ENDDO
108         DO jk2 = 1,jriver_dep
109            friver_dep(jk2,jk) = friver_dep(jk2,jk) / tfthk
110         ENDDO
111      ENDDO
112!!
113!! Have a look at the result of this for a single depth (jriver_dep + 1)
114      IF(lwp) THEN
115          WRITE(numout,*) '=== River nutrient fraction by depth (for a water column of jpk depth)'
116          DO jk = 1,jpk
117             WRITE(numout,*)     &
118             &   ' cell = ', jk, ', friver_dep value = ', friver_dep(jk,jpk)
119          ENDDO
120          IF(lwp) CALL flush(numout)
121       ENDIF
122
123#if defined key_roam
124!! ROAM 3D and 2D carbonate system fields (calculated on first time
125!! step, then monthly)
126      f3_pH(:,:,:)    = 0.
127      f3_h2co3(:,:,:) = 0.
128      f3_hco3(:,:,:)  = 0.
129      f3_co3(:,:,:)   = 0.
130      f3_omcal(:,:,:) = 0.
131      f3_omarg(:,:,:) = 0.
132!!
133      f2_ccd_cal(:,:) = 0.
134      f2_ccd_arg(:,:) = 0.
135      IF(lwp) WRITE(numout,*) ' trc_ini_medusa: carbonate fields initialised to zero'
136#endif
137      IF(lwp) CALL flush(numout)
138
139      !!----------------------------------------------------------------------
140      !! State variable initial conditions (all mmol / m3)
141      !!----------------------------------------------------------------------
142      !!     
143      !! biological and detrital components are initialised to nominal
144      !! values above 100 m depth and zero below; the latter condition
145      !! is applied since non-linear loss processes allow significant
146      !! concentrations of these components to persist at depth
147      !!
148      trn(:,:,:,jpchn) = 0.
149      trn(:,:,:,jpchd) = 0.
150      trn(:,:,:,jpphn) = 0.
151      trn(:,:,:,jpphd) = 0.
152      trn(:,:,:,jppds) = 0.
153      trn(:,:,:,jpzmi) = 0.
154      trn(:,:,:,jpzme) = 0.
155      trn(:,:,:,jpdet) = 0.
156      !!
157      DO jk = 1,13
158         !! non-diatom chlorophyll         (nominal)
159         trn(:,:,jk,jpchn) = 0.01
160         !!
161         !! diatom chlorophyll             (nominal)
162         trn(:,:,jk,jpchd) = 0.01
163         !!
164         !! non-diatom                     (nominal)
165         trn(:,:,jk,jpphn) = 0.01
166         !!
167         !! diatom                         (nominal)
168         trn(:,:,jk,jpphd) = 0.01
169         !!
170         !! diatom silicon                 (nominal)
171         trn(:,:,jk,jppds) = 0.01
172         !!
173         !! microzooplankton               (nominal)
174         trn(:,:,jk,jpzmi) = 0.01
175         !!
176         !! mesozooplankton                (nominal)
177         trn(:,:,jk,jpzme) = 0.01
178         !!
179         !! detrital nitrogen              (nominal)
180         trn(:,:,jk,jpdet) = 0.01
181      ENDDO
182      !!
183      !! dissolved inorganic nitrogen     (nominal average value; typically initialised from climatology)
184      trn(:,:,:,jpdin) = 30.
185      !!
186      !! dissolved silicic acid           (nominal average value; typically initialised from climatology)
187      trn(:,:,:,jpsil) = 90.
188      !!
189      !! dissolved "total" iron           (nominal; typically initialised from model-derived climatology)
190      trn(:,:,:,jpfer) = 1.0e-4           !! = 0.1 umol Fe / m3
191      !!
192      IF(lwp) WRITE(numout,*) ' trc_ini_medusa: MEDUSA-1 fields initialised to defaults'
193# if defined key_roam
194      !!
195      !! detrital carbon                  (nominal)
196      trn(:,:,:,jpdtc) = 0.
197      DO jk = 1,13
198         trn(:,:,jk,jpdtc) = 0.06625
199      ENDDO
200      !!
201      !! dissolved inorganic carbon (DIC) (nominal average value; typically initialised from climatology)
202      trn(:,:,:,jpdic) = 2330.
203      !!
204      !! total alkalinity                 (nominal average value; typically initialised from climatology)
205      trn(:,:,:,jpalk) = 2450.
206      !!
207      !! dissolved oxygen                 (nominal average value; typically initialised from climatology)
208      trn(:,:,:,jpoxy) = 175.
209      !!
210      IF(lwp) WRITE(numout,*) ' trc_ini_medusa: MEDUSA-2 fields initialised to defaults'
211# endif
212      IF(lwp) CALL flush(numout)
213
214      !!----------------------------------------------------------------------
215      !! Sediment pools initial conditions (all mmol / m2)
216      !!----------------------------------------------------------------------
217      !!     
218      !! these pools store biogenic material that has sunk to the seabed,
219      !! and act as a temporary reservoir
220      zb_sed_n(:,:)  = 0.0  !! organic N
221      zn_sed_n(:,:)  = 0.0
222      za_sed_n(:,:)  = 0.0
223      zb_sed_fe(:,:) = 0.0  !! organic Fe
224      zn_sed_fe(:,:) = 0.0
225      za_sed_fe(:,:) = 0.0
226      zb_sed_si(:,:) = 0.0  !! inorganic Si
227      zn_sed_si(:,:) = 0.0
228      za_sed_si(:,:) = 0.0
229      zb_sed_c(:,:)  = 0.0  !! organic C
230      zn_sed_c(:,:)  = 0.0
231      za_sed_c(:,:)  = 0.0
232      zb_sed_ca(:,:) = 0.0  !! inorganic C
233      zn_sed_ca(:,:) = 0.0
234      za_sed_ca(:,:) = 0.0
235      !!
236      IF(lwp) WRITE(numout,*) ' trc_ini_medusa: benthic fields initialised to zero'
237      IF(lwp) CALL flush(numout)
238     
239      !!----------------------------------------------------------------------
240      !! Averaged properties for DMS calculations (various units)
241      !!----------------------------------------------------------------------
242      !!     
243      !! these store temporally averaged properties for DMS calculations (AXY, 07/07/15)
244      zb_dms_chn(:,:)  = 0.0  !! CHN
245      zn_dms_chn(:,:)  = 0.0
246      za_dms_chn(:,:)  = 0.0
247      zb_dms_chd(:,:)  = 0.0  !! CHD
248      zn_dms_chd(:,:)  = 0.0
249      za_dms_chd(:,:)  = 0.0
250      zb_dms_mld(:,:)  = 0.0  !! MLD
251      zn_dms_mld(:,:)  = 0.0
252      za_dms_mld(:,:)  = 0.0
253      zb_dms_qsr(:,:)  = 0.0  !! QSR
254      zn_dms_qsr(:,:)  = 0.0
255      za_dms_qsr(:,:)  = 0.0
256      zb_dms_din(:,:)  = 0.0  !! DIN
257      zn_dms_din(:,:)  = 0.0
258      za_dms_din(:,:)  = 0.0
259      !!
260      IF(lwp) WRITE(numout,*) ' trc_ini_medusa: average fields for DMS initialised to zero'
261      IF(lwp) CALL flush(numout)
262      !!
263      !!---------------------------------------------------------------------
264      !!JPALM (14-06-2016): init dms and co2 flux for coupling with atm (UKESM)
265      !!---------------------------------------------------------------------
266      !!
267      zb_co2_flx(:,:)  = 0.0  !! CO2 flx
268      zn_co2_flx(:,:)  = 0.0
269      za_co2_flx(:,:)  = 0.0
270      zb_dms_srf(:,:)  = 0.0  !! DMS srf
271      zn_dms_srf(:,:)  = 0.0
272      za_dms_srf(:,:)  = 0.0
273      zn_chl_srf(:,:)  = 2.0E-8 !! Chl cpl - set first as surf
274      !!
275      IF(lwp) WRITE(numout,*) ' trc_ini_medusa: DMS and CO2 flux (UKESM) initialised to zero'
276      IF(lwp) CALL flush(numout)
277      IF (lk_oasis) THEN
278         CO2Flux_out_cpl(:,:) =  zn_co2_flx(:,:)   !! Coupling variable
279         DMS_out_cpl(:,:)     =  zn_dms_srf(:,:)   !! Coupling variable
280         chloro_out_cpl(:,:)  =  zn_chl_srf(:,:) * scl_chl   !! Coupling variable
281      END IF
282      !!
283      !!----------------------------------------------------------------------
284      !! AXY (04/11/13): initialise fields previously done by trc_sed_medusa
285      !!----------------------------------------------------------------------
286      !!     
287      IF(lwp) WRITE(numout,*) ' trc_ini_medusa: initialising dust deposition fields'
288      CALL trc_sed_medusa_sbc( nit000 )
289      !!
290      IF(lwp) WRITE(numout,*) ' trc_ini_medusa: initialising ocean CCD array'
291      CALL trc_ini_medusa_ccd( nit000 )
292      fq0 = MINVAL(ocal_ccd(:,:))
293      fq1 = MAXVAL(ocal_ccd(:,:))
294      if (lwp) write (numout,'(a,f10.3,a,f10.3)') & 
295         & 'CCD: min ', fq0, ' max ', fq1
296      !!
297      IF(lwp) WRITE(numout,*) ' trc_ini_medusa: initialising input_swh array'
298      CALL trc_ini_medusa_swh( nit000 )
299      !!
300      IF(lwp) WRITE(numout,*) ' trc_ini_medusa: initialising riverine nutrient arrays'
301      riv_n(:,:)   = 0.0
302      riv_si(:,:)  = 0.0
303      riv_c(:,:)   = 0.0
304      riv_alk(:,:) = 0.0 
305      !!
306      CALL trc_ini_medusa_river( nit000 )
307      fq0 = MINVAL(riv_n(:,:))
308      fq1 = MAXVAL(riv_n(:,:))
309      if (lwp) write (numout,'(a,f10.3,a,f10.3)') & 
310         & 'RIV_N:   min ', fq0, ' max ', fq1
311      fq0 = MINVAL(riv_si(:,:))
312      fq1 = MAXVAL(riv_si(:,:))
313      if (lwp) write (numout,'(a,f10.3,a,f10.3)') & 
314         & 'RIV_SI:  min ', fq0, ' max ', fq1
315      fq0 = MINVAL(riv_c(:,:))
316      fq1 = MAXVAL(riv_c(:,:))
317      if (lwp) write (numout,'(a,f10.3,a,f10.3)') & 
318         & 'RIV_C:   min ', fq0, ' max ', fq1
319      fq0 = MINVAL(riv_alk(:,:))
320      fq1 = MAXVAL(riv_alk(:,:))
321      if (lwp) write (numout,'(a,f10.3,a,f10.3)') & 
322         & 'RIV_ALK: min ', fq0, ' max ', fq1
323      IF(lwp) CALL flush(numout)
324
325      IF(lwp) WRITE(numout,*)
326      IF(lwp) WRITE(numout,*) ' trc_ini_medusa: MEDUSA initialised'
327      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
328      IF(lwp) CALL flush(numout)
329      !!
330      !!----------------------------------------------------------------------
331      !! JPALM (23-01-2017): new way to initialize CO2-atm for cmip6
332      !!                     initially done in trcsms_medusa
333      !!----------------------------------------------------------------------
334      !!
335      IF( ( .NOT.lk_oasis ) .AND. ( .NOT.lk_pi_co2 ) .AND. ( xobs_xco2a <= 0.0 ) ) THEN
336         IF(lwp) WRITE(numout,*) ' trc_ini_medusa: initialisating atm CO2 record'
337         CALL trc_ini_medusa_co2atm
338      ENDIF
339
340   END SUBROUTINE trc_ini_medusa
341
342   !! ======================================================================
343   !! ======================================================================
344   !! ======================================================================
345
346   !! AXY (25/02/10)
347   SUBROUTINE trc_ini_medusa_ccd(kt)
348
349      !!----------------------------------------------------------------------
350      !!                  ***  ROUTINE trc_ini_medusa_ccd  ***
351      !!
352      !! ** Purpose :   Read CCD field
353      !!
354      !! ** Method  :   Read the file
355      !!
356      !! ** input   :   external netcdf files
357      !!
358      !!----------------------------------------------------------------------
359      !! * arguments
360      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
361
362      !!---------------------------------------------------------------------
363
364      !! Open the file
365      !! -------------
366      !!
367      !!!! JPALM -- 14-09-2015 --
368      !!!!       -- to test on ORCA2 with Christian, no file available, so initiate to 0
369      IF (ln_ccd) THEN
370         IF(lwp) WRITE(numout,*) ' '
371         IF(lwp) WRITE(numout,*) ' **** Routine trc_ini_medusa_ccd'
372         CALL iom_open ( 'ccd_ocal_nemo.nc', numccd )
373         IF(lwp) WRITE(numout,*) ' **** trc_ini_medusa_ccd: ccd_ocal_nemo.nc opened'
374
375      !! Read the data
376      !! -------------
377      !!
378         CALL iom_get ( numccd, jpdom_data, 'OCAL_CCD', ocal_ccd )
379         IF(lwp) WRITE(numout,*) ' **** trc_ini_medusa_ccd: data read'
380
381      !! Close the file
382      !! --------------
383      !!
384         CALL iom_close ( numccd )
385         IF(lwp) WRITE(numout,*) ' **** trc_ini_medusa_ccd: ccd_ocal_nemo.nc closed'
386         IF(lwp) CALL flush(numout)
387      ELSE
388         IF(lwp) WRITE(numout,*) ' '
389         IF(lwp) WRITE(numout,*) ' **** Routine trc_ini_medusa_ccd'
390         IF(lwp) WRITE(numout,*) ' **** trc_ini_medusa_ccd: do not read ccd_ocal_nemo.nc'
391         IF(lwp) WRITE(numout,*) ' **** ln_ccd = FALSE and ocal_ccd = 0.0 ---'
392         ocal_ccd(:,:) = 0.0 
393      ENDIF
394 
395   END SUBROUTINE trc_ini_medusa_ccd
396
397   !! ======================================================================
398   !! ======================================================================
399   !! ======================================================================
400
401   !! AXY (25/02/10)
402   SUBROUTINE trc_ini_medusa_swh(kt)
403
404      !!----------------------------------------------------------------------
405      !!                  ***  ROUTINE trc_ini_medusa_swh  ***
406      !!
407      !! ** Purpose :   Read Hs field
408      !!
409      !! ** Method  :   Read the file
410      !!
411      !! ** input   :   external netcdf files
412      !!
413      !!----------------------------------------------------------------------
414      !! * arguments
415      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
416
417      !!---------------------------------------------------------------------
418
419      !! Open the file
420      !! -------------
421      IF(lwp) WRITE(numout,*) ' '
422      IF(lwp) WRITE(numout,*) ' **** Routine trc_ini_medusa_swh'
423      CALL iom_open ( 'gridded_swh.nc', numswh )
424      IF(lwp) WRITE(numout,*) ' **** trc_ini_medusa_swh: gridded_swh.nc opened'
425
426      !! Read the data
427      !! -------------
428      !!
429      CALL iom_get ( numswh, jpdom_data, 'swh', input_swh )
430      IF(lwp) WRITE(numout,*) ' **** trc_ini_medusa_swh: data read'
431
432      !! Close the file
433      !! --------------
434      !!
435      CALL iom_close ( numswh )
436      IF(lwp) WRITE(numout,*) ' **** trc_ini_medusa_swh: gridded_swh.nc closed'
437      IF(lwp) CALL flush(numout)
438 
439   END SUBROUTINE trc_ini_medusa_swh
440
441   !! ======================================================================
442   !! ======================================================================
443   !! ======================================================================
444
445   !! AXY (26/01/12)
446   SUBROUTINE trc_ini_medusa_river(kt)
447
448      !!----------------------------------------------------------------------
449      !!                  ***  ROUTINE trc_ini_medusa_river  ***
450      !!
451      !! ** Purpose :   Read riverine nutrient fields
452      !!
453      !! ** Method  :   Read the file
454      !!
455      !! ** input   :   external netcdf files
456      !!
457      !!----------------------------------------------------------------------
458      !! * arguments
459      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
460
461      !!---------------------------------------------------------------------
462
463      IF(lwp) THEN
464         WRITE(numout,*) ' '
465         WRITE(numout,*) ' **** Routine trc_ini_medusa_river'
466         WRITE(numout,*) ' '
467      ENDIF
468
469      !! Open and read the files
470      !! -----------------------
471      !!
472      if (jriver_n.gt.0) then
473         if (jriver_n.eq.1) CALL iom_open ( 'river_N_conc_orca100.nc', numriv )
474         if (jriver_n.eq.2) CALL iom_open ( 'river_N_flux_orca100.nc', numriv )
475         CALL iom_get  ( numriv, jpdom_data, 'RIV_N', riv_n )
476         IF(lwp) THEN
477            if (jriver_n.eq.1) WRITE(numout,*) ' **** trc_ini_medusa_river: N CONC data read'
478            if (jriver_n.eq.2) WRITE(numout,*) ' **** trc_ini_medusa_river: N FLUX data read'
479         ENDIF
480         CALL iom_close ( numriv )
481         IF(lwp) WRITE(numout,*) ' **** trc_ini_medusa_river: river N file closed'
482      else
483         IF(lwp) THEN
484            WRITE(numout,*) ' **** trc_ini_medusa_river: N data NOT read'
485         ENDIF
486      endif
487      !!
488      if (jriver_si.gt.0) then
489         if (jriver_si.eq.1) CALL iom_open ( 'river_Si_conc_orca100.nc', numriv )
490         if (jriver_si.eq.2) CALL iom_open ( 'river_Si_flux_orca100.nc', numriv )
491         CALL iom_get  ( numriv, jpdom_data, 'RIV_SI', riv_si )
492         IF(lwp) THEN
493            if (jriver_si.eq.1) WRITE(numout,*) ' **** trc_ini_medusa_river: Si CONC data read'
494            if (jriver_si.eq.2) WRITE(numout,*) ' **** trc_ini_medusa_river: Si FLUX data read'
495         ENDIF
496         CALL iom_close ( numriv )
497         IF(lwp) WRITE(numout,*) ' **** trc_ini_medusa_river: river Si file closed'
498      else
499         IF(lwp) THEN
500            WRITE(numout,*) ' **** trc_ini_medusa_river: Si data NOT read'
501         ENDIF
502      endif
503      !!
504      if (jriver_c.gt.0) then
505         if (jriver_c.eq.1) CALL iom_open ( 'river_C_conc_orca100.nc', numriv )
506         if (jriver_c.eq.2) CALL iom_open ( 'river_C_flux_orca100.nc', numriv )
507         CALL iom_get  ( numriv, jpdom_data, 'RIV_C', riv_c )
508         IF(lwp) THEN
509            if (jriver_c.eq.1) WRITE(numout,*) ' **** trc_ini_medusa_river: C CONC data read'
510            if (jriver_c.eq.2) WRITE(numout,*) ' **** trc_ini_medusa_river: C FLUX data read'
511         ENDIF
512         CALL iom_close ( numriv )
513         IF(lwp) WRITE(numout,*) ' **** trc_ini_medusa_river: river C file closed'
514      else
515         IF(lwp) THEN
516            WRITE(numout,*) ' **** trc_ini_medusa_river: C data NOT read'
517         ENDIF
518      endif
519      !!
520      if (jriver_alk.gt.0) then
521         if (jriver_alk.eq.1) CALL iom_open ( 'river_alk_conc_orca100.nc', numriv )
522         if (jriver_alk.eq.2) CALL iom_open ( 'river_alk_flux_orca100.nc', numriv )
523         CALL iom_get  ( numriv, jpdom_data, 'RIV_ALK', riv_alk )
524         IF(lwp) THEN
525            if (jriver_alk.eq.1) WRITE(numout,*) ' **** trc_ini_medusa_river: alkalinity CONC data read'
526            if (jriver_alk.eq.2) WRITE(numout,*) ' **** trc_ini_medusa_river: alkalinity FLUX data read'
527         ENDIF
528         CALL iom_close ( numriv )
529         IF(lwp) WRITE(numout,*) ' **** trc_ini_medusa_river: river alkalinity file closed'
530      else
531         IF(lwp) THEN
532            WRITE(numout,*) ' **** trc_ini_medusa_river: alkalinity data NOT read'
533         ENDIF
534      endif
535      IF(lwp) CALL flush(numout)
536
537   END SUBROUTINE trc_ini_medusa_river
538   
539   SUBROUTINE trc_ini_medusa_co2atm
540      !!----------------------------------------------------------------------
541      !!                     ***  trc_ini_medusa_co2atm  *** 
542      !!
543      !! ** Purpose :   initialization atmospheric co2 record
544      !!
545      !! ** Method  : - Read the xco2 file
546      !!----------------------------------------------------------------------
547      INTEGER                       ::  jn, jm, io, ierr, inum, iostatus
548      INTEGER, PARAMETER            ::  iskip = 4   ! number of 1st descriptor lines
549      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   zyy !: xCO2 record years
550      CHARACTER (len=10), PARAMETER ::  clname = 'xco2.atm'  !! atm CO2 record file
551      !!----------------------------------------------------------------------
552
553      IF(lwp) WRITE(numout,*)
554      IF(lwp) WRITE(numout,*) ' trc_ini_medusa_co2atm: initialisation of atm CO2 historical record'
555      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
556
557
558      IF(lwp) WRITE(numout,*) 'read of formatted file xco2.atm'
559
560      CALL ctl_opn( inum, clname, 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
561      !!!
562      ! -Compute the number of year in the file
563      ! -File starts in co2_yinit, jn represents the record number in the file.
564      ! -Remove the file head (iskip lines) to jn
565      ! -The year is jn + yinit - 1
566      !! Determine the number of lines in xCO2 input file
567      iostatus = 0
568      jn = 1
569      DO WHILE ( iostatus == 0 )
570        READ(inum,'(1x)', IOSTAT=iostatus, END=100)
571        jn = jn + 1
572      ENDDO
573      IF( iostatus .NE. 0 ) THEN
574        !! Error while reading xCO2 input file
575        CALL ctl_stop('trc_ini_medusa_co2atm: &
576                      & Error on the 1st reading of xco2.atm')
577        RETURN
578      ENDIF
579 100  co2_rec = jn - 1 - iskip
580      IF ( lwp) WRITE(numout,*) '    ', co2_rec ,' years read in the file'
581      !                                ! Allocate CO2 hist arrays
582      ierr = 0 
583      ALLOCATE( hist_pco2(co2_rec),zyy(co2_rec), STAT=ierr )
584      IF( ierr > 0 ) THEN
585         CALL ctl_stop( 'trc_ini_medusa_co2atm: unable to allocate  array' ) 
586         RETURN
587      ENDIF
588
589      REWIND(inum)
590
591      DO jm = 1, iskip        ! Skip over 1st six descriptor lines
592         READ(inum,'(1x)')
593      END DO
594      ! file starts in 1931 do jn represent the year in the century.jhh
595      ! Read file till the end
596      ! allocate start and end year of the file
597      DO jn = 1, co2_rec
598        READ(inum,'(F6.1,F12.7)', IOSTAT=io) zyy(jn), hist_pco2(jn)
599        IF( io .NE. 0 ) THEN
600          !! Error while reading xCO2 input file
601          CALL ctl_stop('trc_ini_medusa_co2atm: &
602                        & Error on the 2nd reading of xco2.atm')
603          RETURN
604        ENDIF
605
606        IF(jn==1) co2_yinit = zyy(jn)
607      END DO
608      co2_yend = co2_yinit + real(co2_rec) - 1.
609
610      IF(lwp) THEN        ! Control print
611         WRITE(numout,*)
612         WRITE(numout,*) 'CO2 hist start year: ', co2_yinit
613         WRITE(numout,*) 'CO2 hist end   year: ', co2_yend
614         WRITE(numout,*) ' Year   xCO2 atm '
615         DO jn = 1, co2_rec
616            WRITE(numout, '(F6.1,F12.7)') zyy(jn), hist_pco2(jn)
617         END DO
618      ENDIF
619
620   END SUBROUTINE trc_ini_medusa_co2atm
621
622
623#else
624   !!----------------------------------------------------------------------
625   !!   Dummy module                                        No MEDUSA model
626   !!----------------------------------------------------------------------
627CONTAINS
628   SUBROUTINE trc_ini_medusa             ! Empty routine
629   END SUBROUTINE trc_ini_medusa
630#endif
631
632   !!======================================================================
633END MODULE trcini_medusa
Note: See TracBrowser for help on using the repository browser.