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/UKMO/NEMO_4.0_mirror_SI3_GPU/src/OCE/DIU – NEMO

source: NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_GPU/src/OCE/DIU/diurnal_bulk.F90 @ 13307

Last change on this file since 13307 was 13307, checked in by andmirek, 4 years ago

Ticket #2482: Few additional corrections

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