source: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/TOP/PISCES/P4Z/p4zflx.F90

Last change on this file was 12377, checked in by acc, 10 months ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

  • Property svn:keywords set to Id
File size: 17.6 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   !! * Substitutions
55#  include "do_loop_substitute.h90"
56   !!----------------------------------------------------------------------
57   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
58   !! $Id$
59   !! Software governed by the CeCILL license (see ./LICENSE)
60   !!----------------------------------------------------------------------
61CONTAINS
62
63   SUBROUTINE p4z_flx ( kt, knt, Kbb, Kmm, Krhs )
64      !!---------------------------------------------------------------------
65      !!                     ***  ROUTINE p4z_flx  ***
66      !!
67      !! ** Purpose :   CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
68      !!
69      !! ** Method  :
70      !!              - Include total atm P correction via Esbensen & Kushnir (1981)
71      !!              - Remove Wanninkhof chemical enhancement;
72      !!              - Add option for time-interpolation of atcco2.txt 
73      !!---------------------------------------------------------------------
74      INTEGER, INTENT(in) ::   kt, knt   !
75      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs      ! time level indices
76      !
77      INTEGER  ::   ji, jj, jm, iind, iindm1
78      REAL(wp) ::   ztc, ztc2, ztc3, ztc4, zws, zkgwan
79      REAL(wp) ::   zfld, zflu, zfld16, zflu16, zfact
80      REAL(wp) ::   zvapsw, zsal, zfco2, zxc2, xCO2approx, ztkel, zfugcoeff
81      REAL(wp) ::   zph, zdic, zsch_o2, zsch_co2
82      REAL(wp) ::   zyr_dec, zdco2dt
83      CHARACTER (len=25) ::   charout
84      REAL(wp), DIMENSION(jpi,jpj) ::   zkgco2, zkgo2, zh2co3, zoflx,  zpco2atm 
85      !!---------------------------------------------------------------------
86      !
87      IF( ln_timing )   CALL timing_start('p4z_flx')
88      !
89      ! SURFACE CHEMISTRY (PCO2 AND [H+] IN
90      !     SURFACE LAYER); THE RESULT OF THIS CALCULATION
91      !     IS USED TO COMPUTE AIR-SEA FLUX OF CO2
92
93      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
94
95      IF( ln_co2int .AND. .NOT.ln_presatmco2 .AND. .NOT.l_co2cpl ) THEN 
96         ! Linear temporal interpolation  of atmospheric pco2.  atcco2.txt has annual values.
97         ! Caveats: First column of .txt must be in years, decimal  years preferably.
98         ! For nn_offset, if your model year is iyy, nn_offset=(years(1)-iyy)
99         ! then the first atmospheric CO2 record read is at years(1)
100         zyr_dec = REAL( nyear + nn_offset, wp ) + REAL( nday_year, wp ) / REAL( nyear_len(1), wp )
101         jm = 1
102         DO WHILE( jm <= nmaxrec .AND. years(jm) < zyr_dec ) ;  jm = jm + 1 ;  END DO
103         iind = jm  ;   iindm1 = jm - 1
104         zdco2dt = ( atcco2h(iind) - atcco2h(iindm1) ) / ( years(iind) - years(iindm1) + rtrn )
105         atcco2  = zdco2dt * ( zyr_dec - years(iindm1) ) + atcco2h(iindm1)
106         satmco2(:,:) = atcco2 
107      ENDIF
108
109      IF( l_co2cpl )   satmco2(:,:) = atm_co2(:,:)
110
111      DO_2D_11_11
112         ! DUMMY VARIABLES FOR DIC, H+, AND BORATE
113         zfact = rhop(ji,jj,1) / 1000. + rtrn
114         zdic  = tr(ji,jj,1,jpdic,Kbb)
115         zph   = MAX( hi(ji,jj,1), 1.e-10 ) / zfact
116         ! CALCULATE [H2CO3]
117         zh2co3(ji,jj) = zdic/(1. + ak13(ji,jj,1)/zph + ak13(ji,jj,1)*ak23(ji,jj,1)/zph**2)
118      END_2D
119
120      ! --------------
121      ! COMPUTE FLUXES
122      ! --------------
123
124      ! FIRST COMPUTE GAS EXCHANGE COEFFICIENTS
125      ! -------------------------------------------
126
127      DO_2D_11_11
128         ztc  = MIN( 35., ts(ji,jj,1,jp_tem,Kmm) )
129         ztc2 = ztc * ztc
130         ztc3 = ztc * ztc2 
131         ztc4 = ztc2 * ztc2 
132         ! Compute the schmidt Number both O2 and CO2
133         zsch_co2 = 2116.8 - 136.25 * ztc + 4.7353 * ztc2 - 0.092307 * ztc3 + 0.0007555 * ztc4
134         zsch_o2  = 1920.4 - 135.6  * ztc + 5.2122 * ztc2 - 0.109390 * ztc3 + 0.0009377 * ztc4
135         !  wind speed
136         zws  = wndm(ji,jj) * wndm(ji,jj)
137         ! Compute the piston velocity for O2 and CO2
138         zkgwan = 0.251 * zws
139         zkgwan = zkgwan * xconv * ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1)
140         ! compute gas exchange for CO2 and O2
141         zkgco2(ji,jj) = zkgwan * SQRT( 660./ zsch_co2 )
142         zkgo2 (ji,jj) = zkgwan * SQRT( 660./ zsch_o2 )
143      END_2D
144
145
146      DO_2D_11_11
147         ztkel = tempis(ji,jj,1) + 273.15
148         zsal  = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35.
149         zvapsw    = EXP(24.4543 - 67.4509*(100.0/ztkel) - 4.8489*LOG(ztkel/100) - 0.000544*zsal)
150         zpco2atm(ji,jj) = satmco2(ji,jj) * ( patm(ji,jj) - zvapsw )
151         zxc2      = ( 1.0 - zpco2atm(ji,jj) * 1E-6 )**2
152         zfugcoeff = EXP( patm(ji,jj) * (chemc(ji,jj,2) + 2.0 * zxc2 * chemc(ji,jj,3) )   &
153         &           / ( 82.05736 * ztkel ))
154         zfco2 = zpco2atm(ji,jj) * zfugcoeff
155
156         ! Compute CO2 flux for the sea and air
157         zfld = zfco2 * chemc(ji,jj,1) * zkgco2(ji,jj)  ! (mol/L) * (m/s)
158         zflu = zh2co3(ji,jj) * zkgco2(ji,jj)                                   ! (mol/L) (m/s) ?
159         oce_co2(ji,jj) = ( zfld - zflu ) * tmask(ji,jj,1) 
160         ! compute the trend
161         tr(ji,jj,1,jpdic,Krhs) = tr(ji,jj,1,jpdic,Krhs) + oce_co2(ji,jj) * rfact2 / e3t(ji,jj,1,Kmm)
162
163         ! Compute O2 flux
164         zfld16 = patm(ji,jj) * chemo2(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s)
165         zflu16 = tr(ji,jj,1,jpoxy,Kbb) * zkgo2(ji,jj)
166         zoflx(ji,jj) = ( zfld16 - zflu16 ) * tmask(ji,jj,1)
167         tr(ji,jj,1,jpoxy,Krhs) = tr(ji,jj,1,jpoxy,Krhs) + zoflx(ji,jj) * rfact2 / e3t(ji,jj,1,Kmm)
168      END_2D
169
170      IF( iom_use("tcflx") .OR. iom_use("tcflxcum") .OR. kt == nitrst   &
171         &                 .OR. (ln_check_mass .AND. kt == nitend) )    &
172         t_oce_co2_flx  = glob_sum( 'p4zflx', oce_co2(:,:) * e1e2t(:,:) * 1000. )                    !  Total Flux of Carbon
173      t_oce_co2_flx_cum = t_oce_co2_flx_cum + t_oce_co2_flx       !  Cumulative Total Flux of Carbon
174!      t_atm_co2_flx     = glob_sum( 'p4zflx', satmco2(:,:) * e1e2t(:,:) )       ! Total atmospheric pCO2
175      t_atm_co2_flx     =  atcco2      ! Total atmospheric pCO2
176 
177      IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
178         WRITE(charout, FMT="('flx ')")
179         CALL prt_ctl_trc_info(charout)
180         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
181      ENDIF
182
183      IF( lk_iomput .AND. knt == nrdttrc ) THEN
184         CALL iom_put( "AtmCo2"  , satmco2(:,:) * tmask(:,:,1) )   ! Atmospheric CO2 concentration
185         CALL iom_put( "Cflx"    , oce_co2(:,:) * 1000. ) 
186         CALL iom_put( "Oflx"    , zoflx(:,:) * 1000.  )
187         CALL iom_put( "Kg"      , zkgco2(:,:) * tmask(:,:,1)  )
188         CALL iom_put( "Dpco2"   , ( zpco2atm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) )
189         CALL iom_put( "pCO2sea" , ( zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) )
190         CALL iom_put( "Dpo2"    , ( atcox * patm(:,:) - atcox * tr(:,:,1,jpoxy,Kbb) / ( chemo2(:,:,1) + rtrn ) ) * tmask(:,:,1) )
191         CALL iom_put( "tcflx"   , t_oce_co2_flx     )   ! molC/s
192         CALL iom_put( "tcflxcum", t_oce_co2_flx_cum )   ! molC
193      ENDIF
194      !
195      IF( ln_timing )   CALL timing_stop('p4z_flx')
196      !
197   END SUBROUTINE p4z_flx
198
199
200   SUBROUTINE p4z_flx_init
201      !!----------------------------------------------------------------------
202      !!                  ***  ROUTINE p4z_flx_init  ***
203      !!
204      !! ** Purpose :   Initialization of atmospheric conditions
205      !!
206      !! ** Method  :   Read the nampisext namelist and check the parameters
207      !!      called at the first timestep (nittrc000)
208      !!
209      !! ** input   :   Namelist nampisext
210      !!----------------------------------------------------------------------
211      INTEGER ::   jm, ios   ! Local integer
212      !!
213      NAMELIST/nampisext/ln_co2int, atcco2, clname, nn_offset
214      !!----------------------------------------------------------------------
215      IF(lwp) THEN
216         WRITE(numout,*)
217         WRITE(numout,*) ' p4z_flx_init : atmospheric conditions for air-sea flux calculation'
218         WRITE(numout,*) ' ~~~~~~~~~~~~'
219      ENDIF
220      !
221      READ  ( numnatp_ref, nampisext, IOSTAT = ios, ERR = 901)
222901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampisext in reference namelist' )
223      READ  ( numnatp_cfg, nampisext, IOSTAT = ios, ERR = 902 )
224902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisext in configuration namelist' )
225      IF(lwm) WRITE ( numonp, nampisext )
226      !
227      IF(lwp) THEN                         ! control print
228         WRITE(numout,*) '   Namelist : nampisext --- parameters for air-sea exchange'
229         WRITE(numout,*) '      reading in the atm pCO2 file or constant value   ln_co2int =', ln_co2int
230      ENDIF
231      !
232      CALL p4z_patm( nit000 )
233      !
234      IF( .NOT.ln_co2int .AND. .NOT.ln_presatmco2 ) THEN
235         IF(lwp) THEN                         ! control print
236            WRITE(numout,*) '         Constant Atmospheric pCO2 value               atcco2    =', atcco2
237         ENDIF
238         satmco2(:,:)  = atcco2      ! Initialisation of atmospheric pco2
239      ELSEIF( ln_co2int .AND. .NOT.ln_presatmco2 ) THEN
240         IF(lwp)  THEN
241            WRITE(numout,*) '         Constant Atmospheric pCO2 value               atcco2    =', atcco2
242            WRITE(numout,*) '         Atmospheric pCO2 value  from file             clname    =', TRIM( clname )
243            WRITE(numout,*) '         Offset model-data start year                  nn_offset =', nn_offset
244         ENDIF
245         CALL ctl_opn( numco2, TRIM( clname) , 'OLD', 'FORMATTED', 'SEQUENTIAL', -1 , numout, lwp )
246         jm = 0                      ! Count the number of record in co2 file
247         DO
248           READ(numco2,*,END=100) 
249           jm = jm + 1
250         END DO
251 100     nmaxrec = jm - 1 
252         ALLOCATE( years  (nmaxrec) )   ;   years  (:) = 0._wp
253         ALLOCATE( atcco2h(nmaxrec) )   ;   atcco2h(:) = 0._wp
254         !
255         REWIND(numco2)
256         DO jm = 1, nmaxrec          ! get  xCO2 data
257            READ(numco2, *)  years(jm), atcco2h(jm)
258            IF(lwp) WRITE(numout, '(f6.0,f7.2)')  years(jm), atcco2h(jm)
259         END DO
260         CLOSE(numco2)
261      ELSEIF( .NOT.ln_co2int .AND. ln_presatmco2 ) THEN
262         IF(lwp) WRITE(numout,*) '    Spatialized Atmospheric pCO2 from an external file'
263      ELSE
264         IF(lwp) WRITE(numout,*) '    Spatialized Atmospheric pCO2 from an external file'
265      ENDIF
266      !
267      oce_co2(:,:)  = 0._wp                ! Initialization of Flux of Carbon
268      t_oce_co2_flx = 0._wp
269      t_atm_co2_flx = 0._wp
270      !
271   END SUBROUTINE p4z_flx_init
272
273
274   SUBROUTINE p4z_patm( kt )
275      !!----------------------------------------------------------------------
276      !!                  ***  ROUTINE p4z_atm  ***
277      !!
278      !! ** Purpose :   Read and interpolate the external atmospheric sea-level pressure
279      !! ** Method  :   Read the files and interpolate the appropriate variables
280      !!
281      !!----------------------------------------------------------------------
282      INTEGER, INTENT(in) ::   kt   ! ocean time step
283      !
284      INTEGER            ::   ierr, ios   ! Local integer
285      CHARACTER(len=100) ::   cn_dir      ! Root directory for location of ssr files
286      TYPE(FLD_N)        ::   sn_patm     ! informations about the fields to be read
287      TYPE(FLD_N)        ::   sn_atmco2   ! informations about the fields to be read
288      !!
289      NAMELIST/nampisatm/ ln_presatm, ln_presatmco2, sn_patm, sn_atmco2, cn_dir
290      !!----------------------------------------------------------------------
291      !
292      IF( kt == nit000 ) THEN    !==  First call kt=nittrc000  ==!
293         !
294         IF(lwp) THEN
295            WRITE(numout,*)
296            WRITE(numout,*) ' p4z_patm : sea-level atmospheric pressure'
297            WRITE(numout,*) ' ~~~~~~~~'
298         ENDIF
299         !
300         READ  ( numnatp_ref, nampisatm, IOSTAT = ios, ERR = 901)
301901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisatm in reference namelist' )
302         READ  ( numnatp_cfg, nampisatm, IOSTAT = ios, ERR = 902 )
303902      IF( ios >  0 )   CALL ctl_nam ( ios , 'nampisatm in configuration namelist' )
304         IF(lwm) WRITE ( numonp, nampisatm )
305         !
306         !
307         IF(lwp) THEN                                 !* control print
308            WRITE(numout,*) '   Namelist : nampisatm --- Atmospheric Pressure as external forcing'
309            WRITE(numout,*) '      constant atmopsheric pressure (F) or from a file (T)  ln_presatm    = ', ln_presatm
310            WRITE(numout,*) '      spatial atmopsheric CO2 for flux calcs                ln_presatmco2 = ', ln_presatmco2
311         ENDIF
312         !
313         IF( ln_presatm ) THEN
314            ALLOCATE( sf_patm(1), STAT=ierr )           !* allocate and fill sf_patm (forcing structure) with sn_patm
315            IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_flx: unable to allocate sf_patm structure' )
316            !
317            CALL fld_fill( sf_patm, (/ sn_patm /), cn_dir, 'p4z_flx', 'Atmospheric pressure ', 'nampisatm' )
318                                   ALLOCATE( sf_patm(1)%fnow(jpi,jpj,1)   )
319            IF( sn_patm%ln_tint )  ALLOCATE( sf_patm(1)%fdta(jpi,jpj,1,2) )
320         ENDIF
321         !                                         
322         IF( ln_presatmco2 ) THEN
323            ALLOCATE( sf_atmco2(1), STAT=ierr )           !* allocate and fill sf_atmco2 (forcing structure) with sn_atmco2
324            IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_flx: unable to allocate sf_atmco2 structure' )
325            !
326            CALL fld_fill( sf_atmco2, (/ sn_atmco2 /), cn_dir, 'p4z_flx', 'Atmospheric co2 partial pressure ', 'nampisatm' )
327                                   ALLOCATE( sf_atmco2(1)%fnow(jpi,jpj,1)   )
328            IF( sn_atmco2%ln_tint )  ALLOCATE( sf_atmco2(1)%fdta(jpi,jpj,1,2) )
329         ENDIF
330         !
331         IF( .NOT.ln_presatm )   patm(:,:) = 1._wp    ! Initialize patm if no reading from a file
332         !
333      ENDIF
334      !
335      IF( ln_presatm ) THEN
336         CALL fld_read( kt, 1, sf_patm )               !* input Patm provided at kt + 1/2
337         patm(:,:) = sf_patm(1)%fnow(:,:,1)                        ! atmospheric pressure
338      ENDIF
339      !
340      IF( ln_presatmco2 ) THEN
341         CALL fld_read( kt, 1, sf_atmco2 )               !* input atmco2 provided at kt + 1/2
342         satmco2(:,:) = sf_atmco2(1)%fnow(:,:,1)                        ! atmospheric pressure
343      ELSE
344         satmco2(:,:) = atcco2    ! Initialize atmco2 if no reading from a file
345      ENDIF
346      !
347   END SUBROUTINE p4z_patm
348
349
350   INTEGER FUNCTION p4z_flx_alloc()
351      !!----------------------------------------------------------------------
352      !!                     ***  ROUTINE p4z_flx_alloc  ***
353      !!----------------------------------------------------------------------
354      ALLOCATE( satmco2(jpi,jpj), patm(jpi,jpj), STAT=p4z_flx_alloc )
355      !
356      IF( p4z_flx_alloc /= 0 )   CALL ctl_stop( 'STOP', 'p4z_flx_alloc : failed to allocate arrays' )
357      !
358   END FUNCTION p4z_flx_alloc
359
360   !!======================================================================
361END MODULE p4zflx
Note: See TracBrowser for help on using the repository browser.