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 NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DIU – NEMO

source: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DIU/diurnal_bulk.F90 @ 12202

Last change on this file since 12202 was 12202, checked in by cetlod, 4 years ago

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