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.
diurnal_bulk.F90 in branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/DIU – NEMO

source: branches/UKMO/dev_r5518_GO6_package_FOAMv14_dsst_CO2/NEMOGCM/NEMO/OPA_SRC/DIU/diurnal_bulk.F90 @ 15473

Last change on this file since 15473 was 15473, checked in by dford, 3 years ago

Use cool skin SST in air-sea CO2 flux calculations.
File size: 11.8 KB
Line 
1MODULE diurnal_bulk
2   !!======================================================================
3   !!                    ***  MODULE  diurnal_bulk  ***
4   !!     Takaya model of diurnal warming (Takaya, 2010)
5   !!=====================================================================
6   !! History :        !  11-10  (J. While)  Original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   diurnal_sst_bulk_init  : initalise diurnal model
11   !!   diurnal_sst_bulk_step  : timestep the diurnal model
12   !!----------------------------------------------------------------------
13   USE par_kind
14   USE phycst
15   USE dom_oce
16   USE lib_mpp
17   USE solfrac_mod
18   USE in_out_manager
19   
20   IMPLICIT NONE
21   
22   !Namelist parameters
23   LOGICAL, PUBLIC :: ln_diurnal = .FALSE.
24   LOGICAL, PUBLIC :: ln_diurnal_only = .FALSE.
25
26   !Parameters
27   REAL(wp), PRIVATE, PARAMETER :: pp_alpha = 2.0e-4_wp
28   REAL(wp), PRIVATE, PARAMETER :: pp_veltol = 0._wp
29   REAL(wp), PRIVATE, PARAMETER :: pp_min_fvel = 1.e-10_wp 
30   
31   !Key variables
32   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_dsst     !Delta SST
33   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: x_solfrac  !Fraction of
34                                                      !absorbed
35                                                      !radiation
36
37   PRIVATE
38   PUBLIC diurnal_sst_bulk_init, diurnal_sst_takaya_step
39   
40   CONTAINS
41   
42   SUBROUTINE diurnal_sst_bulk_init
43      !!----------------------------------------------------------------------
44      !! *** ROUTINE diurnal_sst_init ***
45      !!
46      !! ** Purpose : Initalise the Takaya diurnal model
47     
48      !!----------------------------------------------------------------------
49     
50      IMPLICIT NONE
51
52      INTEGER :: ios                 ! Local integer output status for namelist read
53     
54      NAMELIST /namdiu/ ln_diurnal, ln_diurnal_only
55     
56      !Read the namelist
57      !REWIND( numnam )
58      !READ( numnam, namdiu )
59
60      REWIND( numnam_ref )              ! Namelist namdiu in reference namelist
61      READ  ( numnam_ref, namdiu, IOSTAT = ios, ERR = 901)
62901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdiu in reference namelist', lwp )
63
64      REWIND( numnam_cfg )              ! Namelist namdiu in configuration namelist
65      READ  ( numnam_cfg, namdiu, IOSTAT = ios, ERR = 902 )
66902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdiu in configuration namelist', lwp )
67     
68      IF ( ln_diurnal_only .AND. ( .NOT. ln_diurnal ) ) THEN
69         CALL ctl_stop( "ln_diurnal_only set, but ln_diurnal = FALSE !" )
70      ENDIF
71     
72      IF ( ln_diurnal ) THEN     
73         
74         !Allocate arrays
75         ALLOCATE( x_dsst(jpi,jpj), x_solfrac(jpi,jpj) )
76         
77         !Initalise the solar fraction
78         x_solfrac=0._wp
79
80         IF ( ln_diurnal_only ) THEN
81            CALL ctl_warn( "ln_diurnal_only set; only the diurnal component of SST will be calculated" )
82         ENDIF
83      ENDIF
84     
85   END SUBROUTINE diurnal_sst_bulk_init
86   
87   SUBROUTINE diurnal_sst_takaya_step(psolflux, pqflux, ptauflux, prho, p_rdt,&
88            &                  pla, pthick, pcoolthick, pmu, ld_calcfrac, &
89            &                  p_fvel_bkginc, p_hflux_bkginc)
90      !!----------------------------------------------------------------------
91      !! *** ROUTINE diurnal_sst_takaya_step ***
92      !!
93      !! ** Purpose :   Timestep the Takaya diurnal model
94      !!
95      !! ** Method :    1) Calculate the Obukhov length
96      !!                2) Calculate the Simularity function
97      !!                2) Calculate the increment to dsst
98      !!                3) Apply the increment
99      !! ** Reference : Refinements to a prognostic scheme of skin sea surface
100      !!                temperature, Takaya et al, JGR, 2010
101      !!----------------------------------------------------------------------
102     
103      IMPLICIT NONE
104     
105      !Dummy variables
106      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: psolflux  !solar flux (Watts)
107      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: pqflux    !heat (non-solar)
108                                                            !flux (Watts)
109      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: ptauflux  !wind stress
110                                                            !(kg/ m s^2)
111      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: prho      !water density
112                                                            !(kg/m^3)
113      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: pLa   
114                                                     !Langmuir number
115      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: pthick 
116                                                     !warm layer thickness (m)
117      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: pcoolthick 
118                                                     !cool skin thickness (m)
119      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: pmu   
120                                                     !mu parameter
121      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: p_hflux_bkginc 
122                                                     !optional increment to the
123                                                     !heat flux
124      REAL(wp), OPTIONAL, INTENT(IN), DIMENSION(jpi,jpj) :: p_fvel_bkginc 
125                                                     !optional increment to the
126                                                     !friction velocity
127      REAL(wp), INTENT(IN) :: p_rdt              !timestep
128     
129      LOGICAL, OPTIONAL, &
130      &     INTENT(IN)  :: ld_calcfrac !Set TRUE to recalculate the
131                                               !solar fraction
132     
133      !Local variables
134      REAL(wp), DIMENSION(jpi,jpj) :: z_fvel              ! friction velocity     
135      REAL(wp), DIMENSION(jpi,jpj) :: zthick, zcoolthick, zmu, zla
136      REAL(wp), DIMENSION(jpi,jpj) :: z_abflux            ! absorbed flux           
137      REAL(wp), DIMENSION(jpi,jpj) :: z_fla               !Langmuir fuction value
138     
139      LOGICAL  :: ll_calcfrac
140     
141      INTEGER :: ji,jj
142
143      !Set optional arguments to their defualts
144      IF ( .NOT. PRESENT(pthick) ) THEN
145         zthick(:,:) = 3._wp
146      ELSE
147         zthick(:,:) = pthick(:,:)
148      ENDIF
149      IF ( .NOT. PRESENT(pcoolthick) ) THEN
150         zcoolthick(:,:) = 0._wp
151      ELSE
152         zcoolthick(:,:) = pcoolthick(:,:)
153      ENDIF
154      IF ( .NOT. PRESENT(pmu) ) THEN
155         zmu(:,:) = 0.3_wp
156      ELSE
157         zmu(:,:) = pmu(:,:)
158      ENDIF
159      IF ( .NOT. PRESENT(pla) ) THEN
160         zla(:,:) = 0.3_wp
161      ELSE
162         zla(:,:) = pla(:,:)
163      ENDIF
164      IF ( .NOT. PRESENT(ld_calcfrac) ) THEN
165         ll_calcfrac = .FALSE.
166      ELSE
167         ll_calcfrac = ld_calcfrac
168      ENDIF     
169     
170      !If not done already, calculate the solar fraction
171      IF (ll_calcfrac ) THEN
172         DO jj = 1,jpj
173            DO ji = 1, jpi
174               IF(  ( x_solfrac(ji,jj) == 0. ) .AND. ( tmask(ji,jj,1) == 1. ) ) &
175               &   x_solfrac(ji,jj) = solfrac( zcoolthick(ji,jj),zthick(ji,jj) )
176            END DO
177         END DO   
178      ENDIF   
179
180      !convert solar flux and heat flux to absorbed flux   
181      WHERE ( tmask(:,:,1) == 1.) 
182         z_abflux(:,:) = (  x_solfrac(:,:) *  psolflux (:,:)) + pqflux(:,:)     
183      ELSEWHERE
184         z_abflux(:,:) = 0._wp
185      ENDWHERE
186      IF( PRESENT(p_hflux_bkginc) ) z_abflux(:,:) = z_abflux(:,:) + p_hflux_bkginc !Optional increment
187      WHERE ( ABS( z_abflux(:,:) ) < rsmall )
188         z_abflux(:,:) = rsmall
189      ENDWHERE 
190     
191      !Calculate the friction velocity
192      WHERE ( (ptauflux /= 0) .AND. ( tmask(:,:,1) == 1.) )   
193         z_fvel(:,:) = SQRT( ptauflux(:,:) / prho(:,:) )
194      ELSEWHERE
195         z_fvel(:,:) = 0.
196      ENDWHERE
197      IF( PRESENT(p_hflux_bkginc) ) z_fvel(:,:) = z_fvel(:,:) + p_fvel_bkginc !Optional increment
198     
199       
200       
201      !Calculate the Langmuir function value
202      WHERE ( tmask(:,:,1) == 1.)
203         z_fla(:,:) = MAX( 1._wp, zla(:,:)**( -2._wp / 3._wp ) ) 
204      ENDWHERE     
205     
206      !!Increment the temperature using the implicit solution
207      x_dsst(:,:) = t_imp( x_dsst(:,:), p_rdt, z_abflux(:,:), z_fvel(:,:), &
208                                 z_fla(:,:), zmu(:,:), zthick(:,:), prho(:,:))
209     
210                               
211     
212   END SUBROUTINE diurnal_sst_takaya_step
213
214   
215   FUNCTION t_imp(p_dsst, p_rdt, p_abflux, p_fvel, &
216                          p_fla, pmu, pthick, prho )
217                         
218      IMPLICIT NONE
219     
220      !Function definition
221      REAL(wp), DIMENSION(jpi,jpj) :: t_imp
222      !Dummy variables
223      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_dsst    !Delta SST
224      REAL(wp), INTENT(IN) ::                 p_rdt           !Timestep
225      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_abflux  !Heat forcing
226      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fvel    !Friction velocity
227      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fla     !Langmuir number
228      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pmu       !Structure parameter
229      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pthick    !Layer thickness
230      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: prho      !Water density
231   
232      !Local variables
233      REAL(wp) :: z_olength          ! Obukhov length
234      REAL(wp) :: z_sigma, z_sigma2
235      REAL(wp) :: z_term1, z_term2     
236      REAL(wp) :: z_stabfunc         ! stability function value
237      REAL(wp) :: z_fvel     
238     
239      CHARACTER(200) :: warn_string
240       
241      INTEGER :: ji,jj
242                     
243      DO jj = 1, jpj
244         DO ji = 1, jpi   
245           
246            !Only calculate outside tmask
247            IF ( tmask(ji,jj,1) /= 1. ) THEN
248               t_imp(ji,jj) = 0._wp
249               CYCLE   
250            END IF
251           
252            IF (p_fvel(ji,jj) < pp_min_fvel) THEN
253               z_fvel = pp_min_fvel
254               WRITE(warn_string,*) "diurnal_sst_takaya step: "&
255               &//"friction velocity < minimum\n" &
256               &//"Setting friction velocity =",pp_min_fvel
257               CALL ctl_warn(warn_string)
258               
259            ELSE
260               z_fvel = p_fvel(ji,jj)
261            ENDIF
262                 
263            !Calculate the Obukhov length
264            IF ( (z_fvel < pp_veltol ) .AND. &
265            &    (p_dsst(ji,jj) > 0.) ) THEN
266               z_olength =  z_fvel  / &
267               &     SQRT( p_dsst(ji,jj) * vkarmn * grav * &
268               &             pp_alpha / ( 5._wp * pthick(ji,jj) ) )
269            ELSE
270               z_olength = &
271               &   ( prho(ji,jj) * rcp * z_fvel**3._wp ) / &
272               &   ( vkarmn * grav * pp_alpha *&
273               &     p_abflux(ji,jj) )
274            ENDIF
275             
276            !Calculate the stability function         
277            z_sigma = pthick(ji,jj) / z_olength
278            z_sigma2 = z_sigma * z_sigma
279 
280            IF ( z_sigma >= 0 ) THEN
281               z_stabfunc = 1._wp + &
282               &  ( ( 5._wp * z_sigma + 4._wp * z_sigma2 ) / &
283               &  ( 1._wp + 3._wp * z_sigma + 0.25_wp * &
284               &    z_sigma2 ) )
285            ELSE
286               z_stabfunc = 1._wp / &
287               &                   SQRT( 1._wp - 16._wp * z_sigma )
288            ENDIF
289
290            !Calculate the T increment
291            z_term1 = ( p_abflux(ji,jj) * ( pmu(ji,jj) + 1._wp)  / &
292            & ( pmu(ji,jj) * pthick(ji,jj) * prho(ji,jj) * rcp ) )
293           
294             
295            z_term2 = -( ( pmu(ji,jj) + 1._wp) * &
296            &                       ( vkarmn * z_fvel * p_fla(ji,jj) ) / &
297            &      ( pthick(ji,jj) * z_stabfunc ) )     
298         
299            t_imp(ji,jj) = ( p_dsst(ji,jj) + p_rdt * z_term1 ) / &
300                           ( 1._wp - p_rdt * z_term2 )
301
302         END DO
303      END DO
304     
305      END FUNCTION t_imp
306
307END MODULE diurnal_bulk
Note: See TracBrowser for help on using the repository browser.