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_under_ice_relax_dr_hook/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcini_medusa.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 4 years ago

The Dr Hook changes from my perl code.

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