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.
p4zflx.F90 in branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z – NEMO

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90 @ 4148

Last change on this file since 4148 was 4148, checked in by cetlod, 10 years ago

merge in trunk changes between r3853 and r3940 and commit the changes, see ticket #1169

File size: 17.5 KB
Line 
1MODULE p4zflx
2   !!======================================================================
3   !!                         ***  MODULE p4zflx  ***
4   !! TOP :   PISCES CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
5   !!======================================================================
6   !! History :    -   !  1988-07  (E. MAIER-REIMER) Original code
7   !!              -   !  1998     (O. Aumont) additions
8   !!              -   !  1999     (C. Le Quere) modifications
9   !!             1.0  !  2004     (O. Aumont) modifications
10   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
11   !!                  !  2011-02  (J. Simeon, J. Orr) Include total atm P correction
12   !!----------------------------------------------------------------------
13#if defined key_pisces
14   !!----------------------------------------------------------------------
15   !!   'key_pisces'                                       PISCES bio-model
16   !!----------------------------------------------------------------------
17   !!   p4z_flx       :   CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
18   !!   p4z_flx_init  :   Read the namelist
19   !!   p4z_patm      :   Read sfc atm pressure [atm] for each grid cell
20   !!----------------------------------------------------------------------
21   USE oce_trc                      !  shared variables between ocean and passive tracers
22   USE trc                          !  passive tracers common variables
23   USE sms_pisces                   !  PISCES Source Minus Sink variables
24   USE p4zche                       !  Chemical model
25   USE prtctl_trc                   !  print control for debugging
26   USE iom                          !  I/O manager
27   USE fldread                      !  read input fields
28#if defined key_cpl_carbon_cycle
29   USE sbc_oce, ONLY :  atm_co2     !  atmospheric pCO2               
30#endif
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   p4z_flx 
36   PUBLIC   p4z_flx_init 
37   PUBLIC   p4z_flx_alloc 
38
39   !                               !!** Namelist  nampisext  **
40   REAL(wp)          ::  atcco2     !: pre-industrial atmospheric [co2] (ppm)   
41   LOGICAL           ::  ln_co2int  !: flag to read in a file and interpolate atmospheric pco2 or not
42   CHARACTER(len=34) ::  clname     !: filename of pco2 values
43   INTEGER           ::  nn_offset  !: Offset model-data start year (default = 0)
44
45   !!  Variables related to reading atmospheric CO2 time history   
46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: atcco2h, years
47   INTEGER  :: nmaxrec, numco2
48
49   !                               !!* nampisatm namelist (Atmospheric PRessure) *
50   LOGICAL, PUBLIC ::   ln_presatm  !: ref. pressure: global mean Patm (F) or a constant (F)
51
52   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:)  ::  patm      ! atmospheric pressure at kt                 [N/m2]
53   TYPE(FLD), ALLOCATABLE,       DIMENSION(:)    ::  sf_patm   ! structure of input fields (file informations, fields read)
54
55
56   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: oce_co2   !: ocean carbon flux
57   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: satmco2   !: atmospheric pco2
58
59   REAL(wp) ::  xconv  = 0.01_wp / 3600._wp !: coefficients for conversion
60
61   !!* Substitution
62#  include "top_substitute.h90"
63   !!----------------------------------------------------------------------
64   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
65   !! $Id: p4zflx.F90 3294 2012-01-28 16:44:18Z rblod $
66   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
67   !!----------------------------------------------------------------------
68CONTAINS
69
70   SUBROUTINE p4z_flx ( kt )
71      !!---------------------------------------------------------------------
72      !!                     ***  ROUTINE p4z_flx  ***
73      !!
74      !! ** Purpose :   CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
75      !!
76      !! ** Method  :
77      !!              - Include total atm P correction via Esbensen & Kushnir (1981)
78      !!              - Pressure correction NOT done for key_cpl_carbon_cycle
79      !!              - Remove Wanninkhof chemical enhancement;
80      !!              - Add option for time-interpolation of atcco2.txt 
81      !!---------------------------------------------------------------------
82      !
83      INTEGER, INTENT(in) ::   kt   !
84      !
85      INTEGER  ::   ji, jj, jm, iind, iindm1
86      REAL(wp) ::   ztc, ztc2, ztc3, zws, zkgwan
87      REAL(wp) ::   zfld, zflu, zfld16, zflu16, zfact
88      REAL(wp) ::   zph, zah2, zbot, zdic, zalk, zsch_o2, zalka, zsch_co2
89      REAL(wp) ::   zyr_dec, zdco2dt
90      CHARACTER (len=25) :: charout
91      REAL(wp), POINTER, DIMENSION(:,:) :: zkgco2, zkgo2, zh2co3, zoflx 
92      !!---------------------------------------------------------------------
93      !
94      IF( nn_timing == 1 )  CALL timing_start('p4z_flx')
95      !
96      CALL wrk_alloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx )
97      !
98
99      ! SURFACE CHEMISTRY (PCO2 AND [H+] IN
100      !     SURFACE LAYER); THE RESULT OF THIS CALCULATION
101      !     IS USED TO COMPUTE AIR-SEA FLUX OF CO2
102
103      IF( kt /= nit000 ) CALL p4z_patm( kt )    ! Get sea-level pressure (E&K [1981] climatology) for use in flux calcs
104
105      IF( ln_co2int ) THEN 
106         ! Linear temporal interpolation  of atmospheric pco2.  atcco2.txt has annual values.
107         ! Caveats: First column of .txt must be in years, decimal  years preferably.
108         ! For nn_offset, if your model year is iyy, nn_offset=(years(1)-iyy)
109         ! then the first atmospheric CO2 record read is at years(1)
110         zyr_dec = REAL( nyear + nn_offset, wp ) + REAL( nday_year, wp ) / REAL( nyear_len(1), wp )
111         jm = 1
112         DO WHILE( jm <= nmaxrec .AND. years(jm) < zyr_dec ) ;  jm = jm + 1 ;  END DO
113         iind = jm  ;   iindm1 = jm - 1
114         zdco2dt = ( atcco2h(iind) - atcco2h(iindm1) ) / ( years(iind) - years(iindm1) + rtrn )
115         atcco2  = zdco2dt * ( zyr_dec - years(iindm1) ) + atcco2h(iindm1)
116         satmco2(:,:) = atcco2 
117      ENDIF
118
119#if defined key_cpl_carbon_cycle
120      satmco2(:,:) = atm_co2(:,:)
121#endif
122
123      DO jm = 1, 10
124!CDIR NOVERRCHK
125         DO jj = 1, jpj
126!CDIR NOVERRCHK
127            DO ji = 1, jpi
128
129               ! DUMMY VARIABLES FOR DIC, H+, AND BORATE
130               zbot  = borat(ji,jj,1)
131               zfact = rhop(ji,jj,1) / 1000. + rtrn
132               zdic  = trn(ji,jj,1,jpdic) / zfact
133               zph   = MAX( hi(ji,jj,1), 1.e-10 ) / zfact
134               zalka = trn(ji,jj,1,jptal) / zfact
135
136               ! CALCULATE [ALK]([CO3--], [HCO3-])
137               zalk  = zalka - (  akw3(ji,jj,1) / zph - zph + zbot / ( 1.+ zph / akb3(ji,jj,1) )  )
138
139               ! CALCULATE [H+] AND [H2CO3]
140               zah2   = SQRT(  (zdic-zalk)**2 + 4.* ( zalk * ak23(ji,jj,1)   &
141                  &                                        / ak13(ji,jj,1) ) * ( 2.* zdic - zalk )  )
142               zah2   = 0.5 * ak13(ji,jj,1) / zalk * ( ( zdic - zalk ) + zah2 )
143               zh2co3(ji,jj) = ( 2.* zdic - zalk ) / ( 2.+ ak13(ji,jj,1) / zah2 ) * zfact
144               hi(ji,jj,1)   = zah2 * zfact
145            END DO
146         END DO
147      END DO
148
149
150      ! --------------
151      ! COMPUTE FLUXES
152      ! --------------
153
154      ! FIRST COMPUTE GAS EXCHANGE COEFFICIENTS
155      ! -------------------------------------------
156
157!CDIR NOVERRCHK
158      DO jj = 1, jpj
159!CDIR NOVERRCHK
160         DO ji = 1, jpi
161            ztc  = MIN( 35., tsn(ji,jj,1,jp_tem) )
162            ztc2 = ztc * ztc
163            ztc3 = ztc * ztc2 
164            ! Compute the schmidt Number both O2 and CO2
165            zsch_co2 = 2073.1 - 125.62 * ztc + 3.6276 * ztc2 - 0.043126 * ztc3
166            zsch_o2  = 1953.4 - 128.0  * ztc + 3.9918 * ztc2 - 0.050091 * ztc3
167            !  wind speed
168            zws  = wndm(ji,jj) * wndm(ji,jj)
169            ! Compute the piston velocity for O2 and CO2
170            zkgwan = 0.3 * zws  + 2.5 * ( 0.5246 + 0.016256 * ztc + 0.00049946  * ztc2 )
171            zkgwan = zkgwan * xconv * ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1)
172# if defined key_degrad
173            zkgwan = zkgwan * facvol(ji,jj,1)
174#endif 
175            ! compute gas exchange for CO2 and O2
176            zkgco2(ji,jj) = zkgwan * SQRT( 660./ zsch_co2 )
177            zkgo2 (ji,jj) = zkgwan * SQRT( 660./ zsch_o2 )
178         END DO
179      END DO
180
181      DO jj = 1, jpj
182         DO ji = 1, jpi
183            ! Compute CO2 flux for the sea and air
184            zfld = satmco2(ji,jj) * patm(ji,jj) * tmask(ji,jj,1) * chemc(ji,jj,1) * zkgco2(ji,jj)   ! (mol/L) * (m/s)
185            zflu = zh2co3(ji,jj) * tmask(ji,jj,1) * zkgco2(ji,jj)                                   ! (mol/L) (m/s) ?
186            oce_co2(ji,jj) = ( zfld - zflu ) * rfact * e1e2t(ji,jj) * tmask(ji,jj,1) * 1000.
187            ! compute the trend
188            tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) / fse3t(ji,jj,1)
189
190            ! Compute O2 flux
191            zfld16 = atcox * patm(ji,jj) * chemc(ji,jj,2) * tmask(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s)
192            zflu16 = trn(ji,jj,1,jpoxy) * tmask(ji,jj,1) * zkgo2(ji,jj)
193            zoflx(ji,jj) = zfld16 - zflu16
194            tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy) + zoflx(ji,jj) / fse3t(ji,jj,1)
195         END DO
196      END DO
197
198      t_oce_co2_flx = t_oce_co2_flx + glob_sum( oce_co2(:,:) )      ! Cumulative Total Flux of Carbon
199      t_atm_co2_flx = glob_sum( satmco2(:,:) * e1e2t(:,:) )         ! Total atmospheric pCO2
200
201      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
202         WRITE(charout, FMT="('flx ')")
203         CALL prt_ctl_trc_info(charout)
204         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
205      ENDIF
206
207      IF( ln_diatrc ) THEN
208         IF( lk_iomput ) THEN
209            CALL iom_put( "Cflx" , oce_co2(:,:) / e1e2t(:,:) / rfact ) 
210            CALL iom_put( "Oflx" , zoflx(:,:) * 1000 * tmask(:,:,1)  )
211            CALL iom_put( "Kg"   , zkgco2(:,:) * tmask(:,:,1) )
212            CALL iom_put( "Dpco2", ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) )
213            CALL iom_put( "Dpo2" , ( atcox * patm(:,:) - trn(:,:,1,jpoxy) / ( chemc(:,:,2) + rtrn ) ) * tmask(:,:,1) )
214         ELSE
215            trc2d(:,:,jp_pcs0_2d    ) = oce_co2(:,:) / e1e2t(:,:) / rfact 
216            trc2d(:,:,jp_pcs0_2d + 1) = zoflx(:,:) * 1000 * tmask(:,:,1) 
217            trc2d(:,:,jp_pcs0_2d + 2) = zkgco2(:,:) * tmask(:,:,1) 
218            trc2d(:,:,jp_pcs0_2d + 3) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) 
219         ENDIF
220      ENDIF
221      !
222      CALL wrk_dealloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx )
223      !
224      IF( nn_timing == 1 )  CALL timing_stop('p4z_flx')
225      !
226   END SUBROUTINE p4z_flx
227
228
229   SUBROUTINE p4z_flx_init
230      !!----------------------------------------------------------------------
231      !!                  ***  ROUTINE p4z_flx_init  ***
232      !!
233      !! ** Purpose :   Initialization of atmospheric conditions
234      !!
235      !! ** Method  :   Read the nampisext namelist and check the parameters
236      !!      called at the first timestep (nittrc000)
237      !! ** input   :   Namelist nampisext
238      !!----------------------------------------------------------------------
239      NAMELIST/nampisext/ln_co2int, atcco2, clname, nn_offset
240      INTEGER :: jm
241      INTEGER :: ios                 ! Local integer output status for namelist read
242      !!----------------------------------------------------------------------
243      !
244
245      REWIND( numnatp_ref )              ! Namelist nampisext in reference namelist : Pisces atm. conditions
246      READ  ( numnatp_ref, nampisext, IOSTAT = ios, ERR = 901)
247901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisext in reference namelist', lwp )
248
249      REWIND( numnatp_cfg )              ! Namelist nampisext in configuration namelist : Pisces atm. conditions
250      READ  ( numnatp_cfg, nampisext, IOSTAT = ios, ERR = 902 )
251902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisext in configuration namelist', lwp )
252      WRITE ( numonp, nampisext )
253      !
254      IF(lwp) THEN                         ! control print
255         WRITE(numout,*) ' '
256         WRITE(numout,*) ' Namelist parameters for air-sea exchange, nampisext'
257         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
258         WRITE(numout,*) '    Choice for reading in the atm pCO2 file or constant value, ln_co2int =', ln_co2int
259         WRITE(numout,*) ' '
260      ENDIF
261      IF( .NOT.ln_co2int ) THEN
262         IF(lwp) THEN                         ! control print
263            WRITE(numout,*) '    Constant Atmospheric pCO2 value  atcco2    =', atcco2
264            WRITE(numout,*) ' '
265         ENDIF
266         satmco2(:,:)  = atcco2      ! Initialisation of atmospheric pco2
267      ELSE
268         IF(lwp)  THEN
269            WRITE(numout,*) '    Atmospheric pCO2 value  from file clname      =', TRIM( clname )
270            WRITE(numout,*) '    Offset model-data start year      nn_offset   =', nn_offset
271            WRITE(numout,*) ' '
272         ENDIF
273         CALL ctl_opn( numco2, TRIM( clname) , 'OLD', 'FORMATTED', 'SEQUENTIAL', -1 , numout, lwp )
274         jm = 0                      ! Count the number of record in co2 file
275         DO
276           READ(numco2,*,END=100) 
277           jm = jm + 1
278         END DO
279 100     nmaxrec = jm - 1 
280         ALLOCATE( years  (nmaxrec) )     ;      years  (:) = 0._wp
281         ALLOCATE( atcco2h(nmaxrec) )     ;      atcco2h(:) = 0._wp
282
283         REWIND(numco2)
284         DO jm = 1, nmaxrec          ! get  xCO2 data
285            READ(numco2, *)  years(jm), atcco2h(jm)
286            IF(lwp) WRITE(numout, '(f6.0,f7.2)')  years(jm), atcco2h(jm)
287         END DO
288         CLOSE(numco2)
289      ENDIF
290      !
291      oce_co2(:,:)  = 0._wp                ! Initialization of Flux of Carbon
292      t_atm_co2_flx = 0._wp
293      t_oce_co2_flx = 0._wp
294      !
295      CALL p4z_patm( nit000 )
296      !
297   END SUBROUTINE p4z_flx_init
298
299   SUBROUTINE p4z_patm( kt )
300
301      !!----------------------------------------------------------------------
302      !!                  ***  ROUTINE p4z_atm  ***
303      !!
304      !! ** Purpose :   Read and interpolate the external atmospheric sea-levl pressure
305      !! ** Method  :   Read the files and interpolate the appropriate variables
306      !!
307      !!----------------------------------------------------------------------
308      !! * arguments
309      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
310      !
311      INTEGER            ::  ierr
312      INTEGER            ::  ios      ! Local integer output status for namelist read
313      CHARACTER(len=100) ::  cn_dir   ! Root directory for location of ssr files
314      TYPE(FLD_N)        ::  sn_patm  ! informations about the fields to be read
315      !!
316      NAMELIST/nampisatm/ ln_presatm, sn_patm, cn_dir
317
318      !                                         ! ----------------------- !
319      IF( kt == nit000 ) THEN                   ! First call kt=nittrc000 !
320         !                                      ! ----------------------- !
321
322         REWIND( numnatp_ref )              ! Namelist nampisatm in reference namelist : Pisces atm. sea level pressure file
323         READ  ( numnatp_ref, nampisatm, IOSTAT = ios, ERR = 901)
324901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisatm in reference namelist', lwp )
325
326         REWIND( numnatp_cfg )              ! Namelist nampisatm in configuration namelist : Pisces atm. sea level pressure file
327         READ  ( numnatp_cfg, nampisatm, IOSTAT = ios, ERR = 902 )
328902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisatm in configuration namelist', lwp )
329         WRITE ( numonp, nampisatm )
330         !
331         !
332         IF(lwp) THEN                                 !* control print
333            WRITE(numout,*)
334            WRITE(numout,*) '   Namelist nampisatm : Atmospheric Pressure as external forcing'
335            WRITE(numout,*) '      constant atmopsheric pressure (F) or from a file (T)  ln_presatm = ', ln_presatm
336            WRITE(numout,*)
337         ENDIF
338         !
339         IF( ln_presatm ) THEN
340            ALLOCATE( sf_patm(1), STAT=ierr )           !* allocate and fill sf_patm (forcing structure) with sn_patm
341            IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_flx: unable to allocate sf_patm structure' )
342            !
343            CALL fld_fill( sf_patm, (/ sn_patm /), cn_dir, 'p4z_flx', 'Atmospheric pressure ', 'nampisatm' )
344                                   ALLOCATE( sf_patm(1)%fnow(jpi,jpj,1)   )
345            IF( sn_patm%ln_tint )  ALLOCATE( sf_patm(1)%fdta(jpi,jpj,1,2) )
346         ENDIF
347         !                                         
348         IF( .NOT.ln_presatm )   patm(:,:) = 1.e0    ! Initialize patm if no reading from a file
349         !
350      ENDIF
351      !
352      IF( ln_presatm ) THEN
353         CALL fld_read( kt, 1, sf_patm )               !* input Patm provided at kt + 1/2
354         patm(:,:) = sf_patm(1)%fnow(:,:,1)                        ! atmospheric pressure
355      ENDIF
356      !
357   END SUBROUTINE p4z_patm
358
359   INTEGER FUNCTION p4z_flx_alloc()
360      !!----------------------------------------------------------------------
361      !!                     ***  ROUTINE p4z_flx_alloc  ***
362      !!----------------------------------------------------------------------
363      ALLOCATE( oce_co2(jpi,jpj), satmco2(jpi,jpj), patm(jpi,jpj), STAT=p4z_flx_alloc )
364      !
365      IF( p4z_flx_alloc /= 0 )   CALL ctl_warn('p4z_flx_alloc : failed to allocate arrays')
366      !
367   END FUNCTION p4z_flx_alloc
368
369#else
370   !!======================================================================
371   !!  Dummy module :                                   No PISCES bio-model
372   !!======================================================================
373CONTAINS
374   SUBROUTINE p4z_flx( kt )                   ! Empty routine
375      INTEGER, INTENT( in ) ::   kt
376      WRITE(*,*) 'p4z_flx: You should not have seen this print! error?', kt
377   END SUBROUTINE p4z_flx
378#endif 
379
380   !!======================================================================
381END MODULE  p4zflx
Note: See TracBrowser for help on using the repository browser.