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 NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/TOP/PISCES/P4Z – NEMO

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/TOP/PISCES/P4Z/p4zflx.F90 @ 11671

Last change on this file since 11671 was 11671, checked in by acc, 5 years ago

Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Final, non-substantive changes to complete this branch. These changes remove all REWIND statements on the old namelist fortran units (now character variables for internal files). These changes have been left until last since they are easily repeated via a script and it may be preferable to use the previous revision for merge purposes and reapply these last changes separately. This branch has been fully SETTE tested.
  • Property svn:keywords set to Id
File size: 18.1 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   !!   p4z_flx       :   CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
14   !!   p4z_flx_init  :   Read the namelist
15   !!   p4z_patm      :   Read sfc atm pressure [atm] for each grid cell
16   !!----------------------------------------------------------------------
17   USE oce_trc        !  shared variables between ocean and passive tracers
18   USE trc            !  passive tracers common variables
19   USE sms_pisces     !  PISCES Source Minus Sink variables
20   USE p4zche         !  Chemical model
21   USE prtctl_trc     !  print control for debugging
22   USE iom            !  I/O manager
23   USE fldread        !  read input fields
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   p4z_flx 
29   PUBLIC   p4z_flx_init 
30   PUBLIC   p4z_flx_alloc 
31
32   !                                 !!** Namelist  nampisext  **
33   REAL(wp)          ::   atcco2      !: pre-industrial atmospheric [co2] (ppm) 
34   LOGICAL           ::   ln_co2int   !: flag to read in a file and interpolate atmospheric pco2 or not
35   CHARACTER(len=34) ::   clname      !: filename of pco2 values
36   INTEGER           ::   nn_offset   !: Offset model-data start year (default = 0)
37
38   !!  Variables related to reading atmospheric CO2 time history   
39   INTEGER                                   ::   nmaxrec, numco2   !
40   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   atcco2h, years    !
41
42   !                                  !!* nampisatm namelist (Atmospheric PRessure) *
43   LOGICAL, PUBLIC ::   ln_presatm     !: ref. pressure: global mean Patm (F) or a constant (F)
44   LOGICAL, PUBLIC ::   ln_presatmco2  !: accounting for spatial atm CO2 in the compuation of carbon flux (T) or not (F)
45
46   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) ::   patm      ! atmospheric pressure at kt                 [N/m2]
47   TYPE(FLD), ALLOCATABLE,       DIMENSION(:)   ::   sf_patm   ! structure of input fields (file informations, fields read)
48   TYPE(FLD), ALLOCATABLE,       DIMENSION(:)   ::   sf_atmco2 ! structure of input fields (file informations, fields read)
49
50   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  satmco2   !: atmospheric pco2
51
52   REAL(wp) ::   xconv  = 0.01_wp / 3600._wp   !: coefficients for conversion
53
54   !!----------------------------------------------------------------------
55   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
56   !! $Id$
57   !! Software governed by the CeCILL license (see ./LICENSE)
58   !!----------------------------------------------------------------------
59CONTAINS
60
61   SUBROUTINE p4z_flx ( kt, knt )
62      !!---------------------------------------------------------------------
63      !!                     ***  ROUTINE p4z_flx  ***
64      !!
65      !! ** Purpose :   CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
66      !!
67      !! ** Method  :
68      !!              - Include total atm P correction via Esbensen & Kushnir (1981)
69      !!              - Remove Wanninkhof chemical enhancement;
70      !!              - Add option for time-interpolation of atcco2.txt 
71      !!---------------------------------------------------------------------
72      INTEGER, INTENT(in) ::   kt, knt   !
73      !
74      INTEGER  ::   ji, jj, jm, iind, iindm1
75      REAL(wp) ::   ztc, ztc2, ztc3, ztc4, zws, zkgwan
76      REAL(wp) ::   zfld, zflu, zfld16, zflu16, zfact
77      REAL(wp) ::   zvapsw, zsal, zfco2, zxc2, xCO2approx, ztkel, zfugcoeff
78      REAL(wp) ::   zph, zdic, zsch_o2, zsch_co2
79      REAL(wp) ::   zyr_dec, zdco2dt
80      CHARACTER (len=25) ::   charout
81      REAL(wp), DIMENSION(jpi,jpj) ::   zkgco2, zkgo2, zh2co3, zoflx,  zpco2atm 
82      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zw2d
83      !!---------------------------------------------------------------------
84      !
85      IF( ln_timing )   CALL timing_start('p4z_flx')
86      !
87      ! SURFACE CHEMISTRY (PCO2 AND [H+] IN
88      !     SURFACE LAYER); THE RESULT OF THIS CALCULATION
89      !     IS USED TO COMPUTE AIR-SEA FLUX OF CO2
90
91      IF( kt /= nit000 .AND. .NOT.l_co2cpl .AND. knt == 1 )   CALL p4z_patm( kt )   ! Get sea-level pressure (E&K [1981] climatology) for use in flux calcs
92
93      IF( ln_co2int .AND. .NOT.ln_presatmco2 .AND. .NOT.l_co2cpl ) THEN 
94         ! Linear temporal interpolation  of atmospheric pco2.  atcco2.txt has annual values.
95         ! Caveats: First column of .txt must be in years, decimal  years preferably.
96         ! For nn_offset, if your model year is iyy, nn_offset=(years(1)-iyy)
97         ! then the first atmospheric CO2 record read is at years(1)
98         zyr_dec = REAL( nyear + nn_offset, wp ) + REAL( nday_year, wp ) / REAL( nyear_len(1), wp )
99         jm = 1
100         DO WHILE( jm <= nmaxrec .AND. years(jm) < zyr_dec ) ;  jm = jm + 1 ;  END DO
101         iind = jm  ;   iindm1 = jm - 1
102         zdco2dt = ( atcco2h(iind) - atcco2h(iindm1) ) / ( years(iind) - years(iindm1) + rtrn )
103         atcco2  = zdco2dt * ( zyr_dec - years(iindm1) ) + atcco2h(iindm1)
104         satmco2(:,:) = atcco2 
105      ENDIF
106
107      IF( l_co2cpl )   satmco2(:,:) = atm_co2(:,:)
108
109      DO jj = 1, jpj
110         DO ji = 1, jpi
111            ! DUMMY VARIABLES FOR DIC, H+, AND BORATE
112            zfact = rhop(ji,jj,1) / 1000. + rtrn
113            zdic  = trb(ji,jj,1,jpdic)
114            zph   = MAX( hi(ji,jj,1), 1.e-10 ) / zfact
115            ! CALCULATE [H2CO3]
116            zh2co3(ji,jj) = zdic/(1. + ak13(ji,jj,1)/zph + ak13(ji,jj,1)*ak23(ji,jj,1)/zph**2)
117         END DO
118      END DO
119
120      ! --------------
121      ! COMPUTE FLUXES
122      ! --------------
123
124      ! FIRST COMPUTE GAS EXCHANGE COEFFICIENTS
125      ! -------------------------------------------
126
127      DO jj = 1, jpj
128         DO ji = 1, jpi
129            ztc  = MIN( 35., tsn(ji,jj,1,jp_tem) )
130            ztc2 = ztc * ztc
131            ztc3 = ztc * ztc2 
132            ztc4 = ztc2 * ztc2 
133            ! Compute the schmidt Number both O2 and CO2
134            zsch_co2 = 2116.8 - 136.25 * ztc + 4.7353 * ztc2 - 0.092307 * ztc3 + 0.0007555 * ztc4
135            zsch_o2  = 1920.4 - 135.6  * ztc + 5.2122 * ztc2 - 0.109390 * ztc3 + 0.0009377 * ztc4
136            !  wind speed
137            zws  = wndm(ji,jj) * wndm(ji,jj)
138            ! Compute the piston velocity for O2 and CO2
139            zkgwan = 0.251 * zws
140            zkgwan = zkgwan * xconv * ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1)
141            ! compute gas exchange for CO2 and O2
142            zkgco2(ji,jj) = zkgwan * SQRT( 660./ zsch_co2 )
143            zkgo2 (ji,jj) = zkgwan * SQRT( 660./ zsch_o2 )
144         END DO
145      END DO
146
147
148      DO jj = 1, jpj
149         DO ji = 1, jpi
150            ztkel = tempis(ji,jj,1) + 273.15
151            zsal  = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35.
152            zvapsw    = EXP(24.4543 - 67.4509*(100.0/ztkel) - 4.8489*LOG(ztkel/100) - 0.000544*zsal)
153            zpco2atm(ji,jj) = satmco2(ji,jj) * ( patm(ji,jj) - zvapsw )
154            zxc2      = ( 1.0 - zpco2atm(ji,jj) * 1E-6 )**2
155            zfugcoeff = EXP( patm(ji,jj) * (chemc(ji,jj,2) + 2.0 * zxc2 * chemc(ji,jj,3) )   &
156            &           / ( 82.05736 * ztkel ))
157            zfco2 = zpco2atm(ji,jj) * zfugcoeff
158
159            ! Compute CO2 flux for the sea and air
160            zfld = zfco2 * chemc(ji,jj,1) * zkgco2(ji,jj)  ! (mol/L) * (m/s)
161            zflu = zh2co3(ji,jj) * zkgco2(ji,jj)                                   ! (mol/L) (m/s) ?
162            oce_co2(ji,jj) = ( zfld - zflu ) * rfact2 * e1e2t(ji,jj) * tmask(ji,jj,1) * 1000.
163            ! compute the trend
164            tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) * rfact2 / e3t_n(ji,jj,1) * tmask(ji,jj,1)
165
166            ! Compute O2 flux
167            zfld16 = patm(ji,jj) * chemo2(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s)
168            zflu16 = trb(ji,jj,1,jpoxy) * zkgo2(ji,jj)
169            zoflx(ji,jj) = ( zfld16 - zflu16 ) * tmask(ji,jj,1)
170            tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy) + zoflx(ji,jj) * rfact2 / e3t_n(ji,jj,1)
171         END DO
172      END DO
173
174      IF( iom_use("tcflx") .OR. iom_use("tcflxcum") .OR. kt == nitrst   &
175         &                 .OR. (ln_check_mass .AND. kt == nitend) )    &
176         t_oce_co2_flx  = glob_sum( 'p4zflx', oce_co2(:,:) )                    !  Total Flux of Carbon
177      t_oce_co2_flx_cum = t_oce_co2_flx_cum + t_oce_co2_flx       !  Cumulative Total Flux of Carbon
178!      t_atm_co2_flx     = glob_sum( 'p4zflx', satmco2(:,:) * e1e2t(:,:) )       ! Total atmospheric pCO2
179      t_atm_co2_flx     =  atcco2      ! Total atmospheric pCO2
180 
181      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
182         WRITE(charout, FMT="('flx ')")
183         CALL prt_ctl_trc_info(charout)
184         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
185      ENDIF
186
187      IF( lk_iomput .AND. knt == nrdttrc ) THEN
188         ALLOCATE( zw2d(jpi,jpj) ) 
189         IF( iom_use( "Cflx"  ) )  THEN
190            zw2d(:,:) = oce_co2(:,:) / e1e2t(:,:) * rfact2r
191            CALL iom_put( "Cflx"     , zw2d ) 
192         ENDIF
193         IF( iom_use( "Oflx"  ) )  THEN
194            zw2d(:,:) =  zoflx(:,:) * 1000 * tmask(:,:,1)
195            CALL iom_put( "Oflx" , zw2d )
196         ENDIF
197         IF( iom_use( "Kg"    ) )  THEN
198            zw2d(:,:) =  zkgco2(:,:) * tmask(:,:,1)
199            CALL iom_put( "Kg"   , zw2d )
200         ENDIF
201         IF( iom_use( "Dpco2" ) ) THEN
202           zw2d(:,:) = ( zpco2atm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1)
203           CALL iom_put( "Dpco2" ,  zw2d )
204         ENDIF
205         IF( iom_use( "Dpo2" ) )  THEN
206           zw2d(:,:) = ( atcox * patm(:,:) - atcox * trb(:,:,1,jpoxy) / ( chemo2(:,:,1) + rtrn ) ) * tmask(:,:,1)
207           CALL iom_put( "Dpo2"  , zw2d )
208         ENDIF
209         CALL iom_put( "tcflx"    , t_oce_co2_flx * rfact2r )   ! molC/s
210         CALL iom_put( "tcflxcum" , t_oce_co2_flx_cum       )   ! molC
211         !
212         DEALLOCATE( zw2d )
213      ENDIF
214      !
215      IF( ln_timing )   CALL timing_stop('p4z_flx')
216      !
217   END SUBROUTINE p4z_flx
218
219
220   SUBROUTINE p4z_flx_init
221      !!----------------------------------------------------------------------
222      !!                  ***  ROUTINE p4z_flx_init  ***
223      !!
224      !! ** Purpose :   Initialization of atmospheric conditions
225      !!
226      !! ** Method  :   Read the nampisext namelist and check the parameters
227      !!      called at the first timestep (nittrc000)
228      !!
229      !! ** input   :   Namelist nampisext
230      !!----------------------------------------------------------------------
231      INTEGER ::   jm, ios   ! Local integer
232      !!
233      NAMELIST/nampisext/ln_co2int, atcco2, clname, nn_offset
234      !!----------------------------------------------------------------------
235      IF(lwp) THEN
236         WRITE(numout,*)
237         WRITE(numout,*) ' p4z_flx_init : atmospheric conditions for air-sea flux calculation'
238         WRITE(numout,*) ' ~~~~~~~~~~~~'
239      ENDIF
240      !
241      READ  ( numnatp_ref, nampisext, IOSTAT = ios, ERR = 901)
242901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampisext in reference namelist' )
243      READ  ( numnatp_cfg, nampisext, IOSTAT = ios, ERR = 902 )
244902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisext in configuration namelist' )
245      IF(lwm) WRITE ( numonp, nampisext )
246      !
247      IF(lwp) THEN                         ! control print
248         WRITE(numout,*) '   Namelist : nampisext --- parameters for air-sea exchange'
249         WRITE(numout,*) '      reading in the atm pCO2 file or constant value   ln_co2int =', ln_co2int
250      ENDIF
251      !
252      CALL p4z_patm( nit000 )
253      !
254      IF( .NOT.ln_co2int .AND. .NOT.ln_presatmco2 ) THEN
255         IF(lwp) THEN                         ! control print
256            WRITE(numout,*) '         Constant Atmospheric pCO2 value               atcco2    =', atcco2
257         ENDIF
258         satmco2(:,:)  = atcco2      ! Initialisation of atmospheric pco2
259      ELSEIF( ln_co2int .AND. .NOT.ln_presatmco2 ) THEN
260         IF(lwp)  THEN
261            WRITE(numout,*) '         Constant Atmospheric pCO2 value               atcco2    =', atcco2
262            WRITE(numout,*) '         Atmospheric pCO2 value  from file             clname    =', TRIM( clname )
263            WRITE(numout,*) '         Offset model-data start year                  nn_offset =', nn_offset
264         ENDIF
265         CALL ctl_opn( numco2, TRIM( clname) , 'OLD', 'FORMATTED', 'SEQUENTIAL', -1 , numout, lwp )
266         jm = 0                      ! Count the number of record in co2 file
267         DO
268           READ(numco2,*,END=100) 
269           jm = jm + 1
270         END DO
271 100     nmaxrec = jm - 1 
272         ALLOCATE( years  (nmaxrec) )   ;   years  (:) = 0._wp
273         ALLOCATE( atcco2h(nmaxrec) )   ;   atcco2h(:) = 0._wp
274         !
275         REWIND(numco2)
276         DO jm = 1, nmaxrec          ! get  xCO2 data
277            READ(numco2, *)  years(jm), atcco2h(jm)
278            IF(lwp) WRITE(numout, '(f6.0,f7.2)')  years(jm), atcco2h(jm)
279         END DO
280         CLOSE(numco2)
281      ELSEIF( .NOT.ln_co2int .AND. ln_presatmco2 ) THEN
282         IF(lwp) WRITE(numout,*) '    Spatialized Atmospheric pCO2 from an external file'
283      ELSE
284         IF(lwp) WRITE(numout,*) '    Spatialized Atmospheric pCO2 from an external file'
285      ENDIF
286      !
287      oce_co2(:,:)  = 0._wp                ! Initialization of Flux of Carbon
288      t_oce_co2_flx = 0._wp
289      t_atm_co2_flx = 0._wp
290      !
291   END SUBROUTINE p4z_flx_init
292
293
294   SUBROUTINE p4z_patm( kt )
295      !!----------------------------------------------------------------------
296      !!                  ***  ROUTINE p4z_atm  ***
297      !!
298      !! ** Purpose :   Read and interpolate the external atmospheric sea-level pressure
299      !! ** Method  :   Read the files and interpolate the appropriate variables
300      !!
301      !!----------------------------------------------------------------------
302      INTEGER, INTENT(in) ::   kt   ! ocean time step
303      !
304      INTEGER            ::   ierr, ios   ! Local integer
305      CHARACTER(len=100) ::   cn_dir      ! Root directory for location of ssr files
306      TYPE(FLD_N)        ::   sn_patm     ! informations about the fields to be read
307      TYPE(FLD_N)        ::   sn_atmco2   ! informations about the fields to be read
308      !!
309      NAMELIST/nampisatm/ ln_presatm, ln_presatmco2, sn_patm, sn_atmco2, cn_dir
310      !!----------------------------------------------------------------------
311      !
312      IF( kt == nit000 ) THEN    !==  First call kt=nittrc000  ==!
313         !
314         IF(lwp) THEN
315            WRITE(numout,*)
316            WRITE(numout,*) ' p4z_patm : sea-level atmospheric pressure'
317            WRITE(numout,*) ' ~~~~~~~~'
318         ENDIF
319         !
320         READ  ( numnatp_ref, nampisatm, IOSTAT = ios, ERR = 901)
321901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisatm in reference namelist' )
322         READ  ( numnatp_cfg, nampisatm, IOSTAT = ios, ERR = 902 )
323902      IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisatm in configuration namelist' )
324         IF(lwm) WRITE ( numonp, nampisatm )
325         !
326         !
327         IF(lwp) THEN                                 !* control print
328            WRITE(numout,*) '   Namelist : nampisatm --- Atmospheric Pressure as external forcing'
329            WRITE(numout,*) '      constant atmopsheric pressure (F) or from a file (T)  ln_presatm    = ', ln_presatm
330            WRITE(numout,*) '      spatial atmopsheric CO2 for flux calcs                ln_presatmco2 = ', ln_presatmco2
331         ENDIF
332         !
333         IF( ln_presatm ) THEN
334            ALLOCATE( sf_patm(1), STAT=ierr )           !* allocate and fill sf_patm (forcing structure) with sn_patm
335            IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_flx: unable to allocate sf_patm structure' )
336            !
337            CALL fld_fill( sf_patm, (/ sn_patm /), cn_dir, 'p4z_flx', 'Atmospheric pressure ', 'nampisatm' )
338                                   ALLOCATE( sf_patm(1)%fnow(jpi,jpj,1)   )
339            IF( sn_patm%ln_tint )  ALLOCATE( sf_patm(1)%fdta(jpi,jpj,1,2) )
340         ENDIF
341         !                                         
342         IF( ln_presatmco2 ) THEN
343            ALLOCATE( sf_atmco2(1), STAT=ierr )           !* allocate and fill sf_atmco2 (forcing structure) with sn_atmco2
344            IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_flx: unable to allocate sf_atmco2 structure' )
345            !
346            CALL fld_fill( sf_atmco2, (/ sn_atmco2 /), cn_dir, 'p4z_flx', 'Atmospheric co2 partial pressure ', 'nampisatm' )
347                                   ALLOCATE( sf_atmco2(1)%fnow(jpi,jpj,1)   )
348            IF( sn_atmco2%ln_tint )  ALLOCATE( sf_atmco2(1)%fdta(jpi,jpj,1,2) )
349         ENDIF
350         !
351         IF( .NOT.ln_presatm )   patm(:,:) = 1._wp    ! Initialize patm if no reading from a file
352         !
353      ENDIF
354      !
355      IF( ln_presatm ) THEN
356         CALL fld_read( kt, 1, sf_patm )               !* input Patm provided at kt + 1/2
357         patm(:,:) = sf_patm(1)%fnow(:,:,1)                        ! atmospheric pressure
358      ENDIF
359      !
360      IF( ln_presatmco2 ) THEN
361         CALL fld_read( kt, 1, sf_atmco2 )               !* input atmco2 provided at kt + 1/2
362         satmco2(:,:) = sf_atmco2(1)%fnow(:,:,1)                        ! atmospheric pressure
363      ELSE
364         satmco2(:,:) = atcco2    ! Initialize atmco2 if no reading from a file
365      ENDIF
366      !
367   END SUBROUTINE p4z_patm
368
369
370   INTEGER FUNCTION p4z_flx_alloc()
371      !!----------------------------------------------------------------------
372      !!                     ***  ROUTINE p4z_flx_alloc  ***
373      !!----------------------------------------------------------------------
374      ALLOCATE( satmco2(jpi,jpj), patm(jpi,jpj), STAT=p4z_flx_alloc )
375      !
376      IF( p4z_flx_alloc /= 0 )   CALL ctl_stop( 'STOP', 'p4z_flx_alloc : failed to allocate arrays' )
377      !
378   END FUNCTION p4z_flx_alloc
379
380   !!======================================================================
381END MODULE p4zflx
Note: See TracBrowser for help on using the repository browser.