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/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/DIU – NEMO

source: branches/2015/dev_r5656_Met_Office_3_diurnalSST/NEMOGCM/NEMO/OPA_SRC/DIU/diurnal_bulk.F90 @ 5676

Last change on this file since 5676 was 5676, checked in by jwhile, 9 years ago

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