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.
trcsms_cfc.F90 in branches/UKMO/dev_r5518_MEDUSA_optim_RH/NEMOGCM/NEMO/TOP_SRC/CFC – NEMO

source: branches/UKMO/dev_r5518_MEDUSA_optim_RH/NEMOGCM/NEMO/TOP_SRC/CFC/trcsms_cfc.F90 @ 7693

Last change on this file since 7693 was 7693, checked in by frrh, 7 years ago

Merge in branches/NERC/dev_r5518_NOC_MEDUSA_Stable from revision range
r 5711:7611 (i.e. the entire branch)
File size: 17.6 KB
Line 
1MODULE trcsms_cfc
2   !!======================================================================
3   !!                      ***  MODULE trcsms_cfc  ***
4   !! TOP : CFC main model
5   !!======================================================================
6   !! History :  OPA  !  1999-10  (JC. Dutay)  original code
7   !!  NEMO      1.0  !  2004-03  (C. Ethe) free form + modularity
8   !!            2.0  !  2007-12  (C. Ethe, G. Madec)  reorganisation
9   !!----------------------------------------------------------------------
10#if defined key_cfc
11   !!----------------------------------------------------------------------
12   !!   'key_cfc'                                               CFC tracers
13   !!----------------------------------------------------------------------
14   !!   trc_sms_cfc  :  compute and add CFC suface forcing to CFC trends
15   !!   cfc_init     :  sets constants for CFC surface forcing computation
16   !!----------------------------------------------------------------------
17   USE dom_oce       ! ocean space and time domain
18   USE oce_trc       ! Ocean variables
19   USE par_trc       ! TOP parameters
20   USE trc           ! TOP variables
21   USE trd_oce
22   USE trdtrc
23   USE iom           ! I/O library
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   trc_sms_cfc         ! called in ???   
29   PUBLIC   trc_sms_cfc_alloc   ! called in trcini_cfc.F90
30
31   INTEGER , PUBLIC, PARAMETER ::   jphem  =   2   ! parameter for the 2 hemispheres
32   INTEGER , PUBLIC            ::   jpyear         ! Number of years read in CFC1112 file
33   INTEGER , PUBLIC            ::   ndate_beg      ! initial calendar date (aammjj) for CFC
34   INTEGER , PUBLIC            ::   simu_type      ! Kind of simulation: 1- Spin-up
35                                                   !                     2- Hindcast/projection
36   INTEGER , PUBLIC            ::   nyear_res      ! restoring time constant (year)
37   INTEGER , PUBLIC            ::   nyear_beg      ! initial year (aa)
38   
39   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   p_cfc    ! partial hemispheric pressure for CFC
40   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   xphem    ! spatial interpolation factor for patm
41   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qtr_cfc  ! flux at surface
42   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qint_cfc ! cumulative flux
43   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   patm     ! atmospheric function
44
45   REAL(wp), DIMENSION(4,2) ::   soa   ! coefficient for solubility of CFC [mol/l/atm]
46   REAL(wp), DIMENSION(3,2) ::   sob   !    "               "
47   REAL(wp), DIMENSION(4,2) ::   sca   ! coefficients for schmidt number in degre Celcius
48     
49   !                          ! coefficients for conversion
50   REAL(wp) ::   xconv1 = 1.0          ! conversion from to
51   REAL(wp) ::   xconv2 = 0.01/3600.   ! conversion from cm/h to m/s:
52   REAL(wp) ::   xconv3 = 1.0e+3       ! conversion from mol/l/atm to mol/m3/atm
53   REAL(wp) ::   xconv4 = 1.0e-12      ! conversion from mol/m3/atm to mol/m3/pptv
54
55   !! * Substitutions
56#  include "top_substitute.h90"
57   !!----------------------------------------------------------------------
58   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
59   !! $Id$
60   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
61   !!----------------------------------------------------------------------
62CONTAINS
63
64   SUBROUTINE trc_sms_cfc( kt )
65      !!----------------------------------------------------------------------
66      !!                     ***  ROUTINE trc_sms_cfc  ***
67      !!
68      !! ** Purpose :   Compute the surface boundary contition on CFC 11
69      !!             passive tracer associated with air-mer fluxes and add it
70      !!             to the general trend of tracers equations.
71      !!
72      !! ** Method  : - get the atmospheric partial pressure - given in pico -
73      !!              - computation of solubility ( in 1.e-12 mol/l then in 1.e-9 mol/m3)
74      !!              - computation of transfert speed ( given in cm/hour ----> cm/s )
75      !!              - the input function is given by :
76      !!                speed * ( concentration at equilibrium - concentration at surface )
77      !!              - the input function is in pico-mol/m3/s and the
78      !!                CFC concentration in pico-mol/m3
79      !!----------------------------------------------------------------------
80      !
81      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
82      !
83      INTEGER  ::   ji, jj, jn, jl, jm, js
84      INTEGER  ::   iyear_beg, iyear_end, iyear_tmp
85      INTEGER  ::   im1, im2, ierr
86      REAL(wp) ::   ztap, zdtap       
87      REAL(wp) ::   zt1, zt2, zt3, zv2
88      REAL(wp) ::   zsol      ! solubility
89      REAL(wp) ::   zsch      ! schmidt number
90      REAL(wp) ::   zpp_cfc   ! atmospheric partial pressure of CFC
91      REAL(wp) ::   zca_cfc   ! concentration at equilibrium
92      REAL(wp) ::   zak_cfc   ! transfert coefficients
93      REAL(wp), ALLOCATABLE, DIMENSION(:,:)  ::   zpatm     ! atmospheric function
94      !!----------------------------------------------------------------------
95      !
96      !
97      IF( nn_timing == 1 )  CALL timing_start('trc_sms_cfc')
98      !
99      ALLOCATE( zpatm(jphem,jp_cfc), STAT=ierr )
100      IF( ierr > 0 ) THEN
101         CALL ctl_stop( 'trc_sms_cfc: unable to allocate zpatm array' )   ;   RETURN
102      ENDIF
103
104      IF( kt == nittrc000 )   CALL cfc_init
105
106      ! Temporal interpolation
107      ! ----------------------
108      !! JPALM -- 15-06-2016 -- define 2 kind of CFC run.
109      !!                     we want to make cycle experiments,
110      !!                     to periodically compare the ocean dynamic within
111      !!                     1- the SPIN-UP and 2- Hincast/Projections
112      !!                     -- main difference is the way to define the year of
113      !!                     simulation, that determine the atm pCFC.
114      !!                     1-- Spin-up: our atm forcing is of 30y we cycle on.
115      !!                     So we do 90y CFC cycles to be in good
116      !!                     correspondance with the atmosphere
117      !!                     2-- Hindcast/proj, instead of nyear-1900 we keep
118      !!                     the 2 last digit, and enable 3 cycle from 1800 to 2100. 
119      !!----------------------------------------------------------------------
120      !! 1 -- SPIN-UP
121      IF (simu_type==1) THEN
122         iyear_tmp = nyear - nyear_res  !! JPALM -- in our spin-up, nyear_res is 1000
123         iyear_beg = MOD( iyear_tmp , 90 )
124         !! JPALM -- the pCFC file only got 78 years.
125         !!       So if iyear_beg > 78 then we set pCFC to 0
126         !!             iyear_beg = 0 as well -- must try to avoid obvious problems
127         !!             as Pcfc is set to 0.00 up to year 32, let set iyear_beg to year 10
128         !!          else, must add 30 to iyear_beg to match with P_cfc indices
129         !!---------------------------------------
130         IF ((iyear_beg > 77) .OR. (iyear_beg==0)) THEN
131            iyear_beg = 10
132         ELSE
133            iyear_beg = iyear_beg + 30
134         ENDIF
135      !!
136      !! 2 -- Hindcast/proj
137      ELSEIF (simu_type==2) THEN
138         iyear_beg = MOD(nyear, 100)
139         IF (iyear_beg < 9)  iyear_beg = iyear_beg + 100
140         !! JPALM -- Same than previously, if iyear_beg is out of P_cfc range,
141         !!       we want to set p_CFC to 0.00 --> set iyear_beg = 10
142         IF ((iyear_beg < 30) .OR. (iyear_beg > 107)) iyear_beg = 10             
143      ENDIF
144      !!
145      IF ( nmonth <= 6 ) THEN
146         iyear_beg = iyear_beg - 1
147         im1       =  6 - nmonth + 1
148         im2       =  6 + nmonth - 1
149      ELSE
150         im1       = 12 - nmonth + 7
151         im2       =      nmonth - 7
152      ENDIF
153      iyear_end = iyear_beg + 1
154
155      !                                                  !------------!
156      DO jl = 1, jp_cfc                                  !  CFC loop  !
157         !                                               !------------!
158         jn = jp_cfc0 + jl - 1
159         ! time interpolation at time kt
160         DO jm = 1, jphem
161            zpatm(jm,jl) = (  p_cfc(iyear_beg, jm, jl) * FLOAT (im1)  &
162               &           +  p_cfc(iyear_end, jm, jl) * FLOAT (im2) ) / 12.
163         END DO
164         
165         !                                                         !------------!
166         DO jj = 1, jpj                                            !  i-j loop  !
167            DO ji = 1, jpi                                         !------------!
168 
169               ! space interpolation
170               zpp_cfc  =       xphem(ji,jj)   * zpatm(1,jl)   &
171                  &     + ( 1.- xphem(ji,jj) ) * zpatm(2,jl)
172
173               ! Computation of concentration at equilibrium : in picomol/l
174               ! coefficient for solubility for CFC-11/12 in  mol/l/atm
175               IF( tmask(ji,jj,1) .GE. 0.5 ) THEN
176                  ztap  = ( tsn(ji,jj,1,jp_tem) + 273.16 ) * 0.01
177                  zdtap = sob(1,jl) + ztap * ( sob(2,jl) + ztap * sob(3,jl) ) 
178                  zsol  =  EXP( soa(1,jl) + soa(2,jl) / ztap + soa(3,jl) * LOG( ztap )   &
179                     &                    + soa(4,jl) * ztap * ztap + tsn(ji,jj,1,jp_sal) * zdtap ) 
180               ELSE
181                  zsol  = 0.e0
182               ENDIF
183               ! conversion from mol/l/atm to mol/m3/atm and from mol/m3/atm to mol/m3/pptv   
184               zsol = xconv4 * xconv3 * zsol * tmask(ji,jj,1) 
185               ! concentration at equilibrium
186               zca_cfc = xconv1 * zpp_cfc * zsol * tmask(ji,jj,1)             
187 
188               ! Computation of speed transfert
189               !    Schmidt number
190               zt1  = tsn(ji,jj,1,jp_tem)
191               zt2  = zt1 * zt1 
192               zt3  = zt1 * zt2
193               zsch = sca(1,jl) + sca(2,jl) * zt1 + sca(3,jl) * zt2 + sca(4,jl) * zt3
194
195               !    speed transfert : formulae of wanninkhof 1992
196               zv2     = wndm(ji,jj) * wndm(ji,jj)
197               zsch    = zsch / 660.
198               zak_cfc = ( 0.39 * xconv2 * zv2 / SQRT(zsch) ) * tmask(ji,jj,1)
199
200               ! Input function  : speed *( conc. at equil - concen at surface )
201               ! trn in pico-mol/l idem qtr; ak in en m/a
202               qtr_cfc(ji,jj,jl) = -zak_cfc * ( trb(ji,jj,1,jn) - zca_cfc )   &
203#if defined key_degrad
204                  &                         * facvol(ji,jj,1)                           &
205#endif
206                  &                         * tmask(ji,jj,1) * ( 1. - fr_i(ji,jj) )
207               ! Add the surface flux to the trend
208               tra(ji,jj,1,jn) = tra(ji,jj,1,jn) + qtr_cfc(ji,jj,jl) / fse3t(ji,jj,1) 
209
210               ! cumulation of surface flux at each time step
211               qint_cfc(ji,jj,jl) = qint_cfc(ji,jj,jl) + qtr_cfc(ji,jj,jl) * rdt
212               !                                               !----------------!
213            END DO                                             !  end i-j loop  !
214         END DO                                                !----------------!
215         !                                                  !----------------!
216      END DO                                                !  end CFC loop  !
217         !
218      IF( kt == nittrc000 ) THEN
219         DO jl = 1, jp_cfc   
220             WRITE(NUMOUT,*) ' '
221             WRITE(NUMOUT,*) 'CFC interpolation verification '  !! Jpalm 
222             WRITE(NUMOUT,*) '################################## '
223             WRITE(NUMOUT,*) ' '
224               if (jl.EQ.1) then
225                   WRITE(NUMOUT,*) 'Traceur = CFC11: '
226               elseif (jl.EQ.2) then
227                   WRITE(NUMOUT,*) 'Traceur = CFC12: '
228               endif
229             WRITE(NUMOUT,*) 'nyear    = ', nyear
230             WRITE(NUMOUT,*) 'nmonth   = ', nmonth
231             WRITE(NUMOUT,*) 'iyear_beg= ', iyear_beg
232             WRITE(NUMOUT,*) 'iyear_end= ', iyear_end
233             WRITE(NUMOUT,*) 'p_cfc(iyear_beg)= ',p_cfc(iyear_beg, 1, jl)
234             WRITE(NUMOUT,*) 'p_cfc(iyear_end)= ',p_cfc(iyear_end, 1, jl)
235             WRITE(NUMOUT,*) 'Im1= ',im1
236             WRITE(NUMOUT,*) 'Im2= ',im2
237             WRITE(NUMOUT,*) 'zpp_cfc = ',zpp_cfc
238             WRITE(NUMOUT,*) ' '
239         END DO 
240# if defined key_debug_medusa
241         CALL flush(numout)
242# endif
243      ENDIF
244        !
245      !IF( lrst_trc ) THEN
246      !   IF(lwp) WRITE(numout,*)
247      !   IF(lwp) WRITE(numout,*) 'trc_sms_cfc : cumulated input function fields written in ocean restart file ',   &
248      !      &                    'at it= ', kt,' date= ', ndastp
249      !   IF(lwp) WRITE(numout,*) '~~~~'
250      !   DO jn = jp_cfc0, jp_cfc1
251      !      CALL iom_rstput( kt, nitrst, numrtw, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jn) )
252      !   END DO
253      !ENDIF                                           
254      !
255      IF( lk_iomput ) THEN
256         CALL iom_put( "qtrCFC11"  , qtr_cfc (:,:,1) )
257         CALL iom_put( "qintCFC11" , qint_cfc(:,:,1) )
258      ELSE
259         IF( ln_diatrc ) THEN
260            trc2d(:,:,jp_cfc0_2d    ) = qtr_cfc (:,:,1)
261            trc2d(:,:,jp_cfc0_2d + 1) = qint_cfc(:,:,1)
262         END IF
263      END IF
264      !
265      IF( l_trdtrc ) THEN
266          DO jn = jp_cfc0, jp_cfc1
267            CALL trd_trc( tra(:,:,:,jn), jn, jptra_sms, kt )   ! save trends
268          END DO
269      END IF
270      !
271# if defined key_debug_medusa
272      IF(lwp) WRITE(numout,*) '   CFC - Check: nn_timing = ', nn_timing
273      CALL flush(numout)
274# endif
275      IF( nn_timing == 1 )  CALL timing_stop('trc_sms_cfc')
276      !
277   END SUBROUTINE trc_sms_cfc
278
279
280   SUBROUTINE cfc_init
281      !!---------------------------------------------------------------------
282      !!                     ***  cfc_init  *** 
283      !!
284      !! ** Purpose : sets constants for CFC model
285      !!---------------------------------------------------------------------
286      INTEGER :: jl, jn, iyear_beg, iyear_tmp
287
288      ! coefficient for CFC11
289      !----------------------
290
291      ! Solubility
292      soa(1,1) = -229.9261 
293      soa(2,1) =  319.6552
294      soa(3,1) =  119.4471
295      soa(4,1) =  -1.39165
296
297      sob(1,1) =  -0.142382
298      sob(2,1) =   0.091459
299      sob(3,1) =  -0.0157274
300
301      ! Schmidt number
302      sca(1,1) = 3501.8
303      sca(2,1) = -210.31
304      sca(3,1) =  6.1851
305      sca(4,1) = -0.07513
306
307      ! coefficient for CFC12
308      !----------------------
309
310      ! Solubility
311      soa(1,2) = -218.0971
312      soa(2,2) =  298.9702
313      soa(3,2) =  113.8049
314      soa(4,2) =  -1.39165
315
316      sob(1,2) =  -0.143566
317      sob(2,2) =   0.091015
318      sob(3,2) =  -0.0153924
319
320      ! schmidt number
321      sca(1,2) =  3845.4 
322      sca(2,2) =  -228.95
323      sca(3,2) =  6.1908 
324      sca(4,2) =  -0.067430
325
326      !!---------------------------------------------
327      !! JPALM -- re-initialize CFC fields and diags if restart a CFC cycle,
328      !!       Or if out of P_cfc range
329      IF (simu_type==1) THEN
330         iyear_tmp = nyear - nyear_res  !! JPALM -- in our spin-up, nyear_res is 1000
331         iyear_beg = MOD( iyear_tmp , 90 )
332         !!---------------------------------------
333         IF ((iyear_beg > 77) .OR. (iyear_beg==0)) THEN
334            qtr_cfc(:,:,:) = 0._wp
335            IF(lwp) THEN
336               WRITE(numout,*) 
337               WRITE(numout,*) 'restart a CFC cycle or out of P_cfc year bounds zero --'
338               WRITE(numout,*) '                          --    set qtr_CFC = 0.00   --'
339               WRITE(numout,*) '                          --   set qint_CFC = 0.00   --'
340               WRITE(numout,*) '                          --   set trn(CFC) = 0.00   --'
341            ENDIF
342            qtr_cfc(:,:,:) = 0._wp
343            qint_cfc(:,:,:) = 0._wp
344            DO jl = 1, jp_cfc
345              jn = jp_cfc0 + jl - 1
346              trn(:,:,:,jn) = 0._wp
347              trb(:,:,:,jn) = 0._wp
348            END DO
349         ENDIF
350      !!
351      !! 2 -- Hindcast/proj
352      ELSEIF (simu_type==2) THEN
353         iyear_beg = MOD(nyear, 100)
354         IF (iyear_beg < 9)  iyear_beg = iyear_beg + 100
355         IF ((iyear_beg < 30) .OR. (iyear_beg > 107)) THEN
356            qtr_cfc(:,:,:) = 0._wp
357            IF(lwp) THEN
358               WRITE(numout,*)
359               WRITE(numout,*) 'restart a CFC cycle or out of P_cfc year bounds zero --'
360               WRITE(numout,*) '                          --    set qtr_CFC = 0.00   --'
361               WRITE(numout,*) '                          --   set qint_CFC = 0.00   --'
362               WRITE(numout,*) '                          --   set trn(CFC) = 0.00   --'
363            ENDIF
364            qtr_cfc(:,:,:) = 0._wp
365            qint_cfc(:,:,:) = 0._wp
366            DO jl = 1, jp_cfc
367              jn = jp_cfc0 + jl - 1
368              trn(:,:,:,jn) = 0._wp
369              trb(:,:,:,jn) = 0._wp
370            END DO
371         ENDIF
372      ENDIF
373
374      IF(lwp) WRITE(numout,*)
375      !
376   END SUBROUTINE cfc_init
377
378
379   INTEGER FUNCTION trc_sms_cfc_alloc()
380      !!----------------------------------------------------------------------
381      !!                     ***  ROUTINE trc_sms_cfc_alloc  ***
382      !!----------------------------------------------------------------------
383      ALLOCATE( xphem   (jpi,jpj)        ,     &
384         &      qtr_cfc (jpi,jpj,jp_cfc) ,     &
385         &      qint_cfc(jpi,jpj,jp_cfc) , STAT=trc_sms_cfc_alloc )
386         !
387      IF( trc_sms_cfc_alloc /= 0 ) CALL ctl_warn('trc_sms_cfc_alloc : failed to allocate arrays.')
388      !
389   END FUNCTION trc_sms_cfc_alloc
390
391#else
392   !!----------------------------------------------------------------------
393   !!   Dummy module                                         No CFC tracers
394   !!----------------------------------------------------------------------
395CONTAINS
396   SUBROUTINE trc_sms_cfc( kt )       ! Empty routine
397      WRITE(*,*) 'trc_sms_cfc: You should not have seen this print! error?', kt
398   END SUBROUTINE trc_sms_cfc
399#endif
400
401   !!======================================================================
402END MODULE trcsms_cfc
Note: See TracBrowser for help on using the repository browser.