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.
diu_bulk.F90 in NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/DIU – NEMO

source: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/DIU/diu_bulk.F90 @ 13900

Last change on this file since 13900 was 12377, checked in by acc, 4 years 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: 10.8 KB
Line 
1MODULE diu_bulk
2   !!======================================================================
3   !!                    ***  MODULE  diu_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   !! * Substitutions
39#  include "do_loop_substitute.h90"
40   
41   !!----------------------------------------------------------------------
42CONTAINS
43   
44   SUBROUTINE diurnal_sst_bulk_init
45      !!----------------------------------------------------------------------
46      !! *** ROUTINE diurnal_sst_init ***
47      !!
48      !! ** Purpose : Initialise the Takaya diurnal model   
49      !!----------------------------------------------------------------------     
50      INTEGER ::   ios   ! local integer
51      !!
52      NAMELIST /namdiu/ ln_diurnal, ln_diurnal_only
53      !!----------------------------------------------------------------------     
54
55      ! Read the namelist
56      READ  ( numnam_ref, namdiu, IOSTAT = ios, ERR = 901 )
57901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdiu in reference namelist' )
58      READ  ( numnam_cfg, namdiu, IOSTAT = ios, ERR = 902 )
59902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdiu in configuration namelist' )     
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_2D_11_11
133            IF(  ( x_solfrac(ji,jj) == 0._wp ) .AND. ( tmask(ji,jj,1) == 1._wp ) ) &
134               &   x_solfrac(ji,jj) = solfrac( zcoolthick(ji,jj),zthick(ji,jj) )
135         END_2D
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_2D_11_11
202         
203         ! Only calculate outside tmask
204         IF ( tmask(ji,jj,1) /= 1._wp ) THEN
205            t_imp(ji,jj) = 0._wp
206            CYCLE   
207         END IF
208         
209         IF (p_fvel(ji,jj) < pp_min_fvel) THEN
210            z_fvel = pp_min_fvel
211            WRITE(warn_string,*) "diurnal_sst_takaya step: "&
212            &//"friction velocity < minimum\n" &
213            &//"Setting friction velocity =",pp_min_fvel
214            CALL ctl_warn(warn_string)
215           
216         ELSE
217            z_fvel = p_fvel(ji,jj)
218         ENDIF
219             
220         ! Calculate the Obukhov length
221         IF ( (z_fvel < pp_veltol ) .AND. &
222         &    (p_dsst(ji,jj) > 0._wp) ) THEN
223            z_olength =  z_fvel  / &
224            &     SQRT( p_dsst(ji,jj) * vkarmn * grav * &
225            &             pp_alpha / ( 5._wp * pthick(ji,jj) ) )
226         ELSE
227            z_olength = &
228            &   ( prho(ji,jj) * rcp * z_fvel**3._wp ) / &
229            &   ( vkarmn * grav * pp_alpha *&
230            &     p_abflux(ji,jj) )
231         ENDIF
232         
233         ! Calculate the stability function         
234         z_sigma = pthick(ji,jj) / z_olength
235         z_sigma2 = z_sigma * z_sigma
236         IF ( z_sigma >= 0. ) THEN
237            z_stabfunc = 1._wp + &
238            &  ( ( 5._wp * z_sigma + 4._wp * z_sigma2 ) / &
239            &  ( 1._wp + 3._wp * z_sigma + 0.25_wp * &
240            &    z_sigma2 ) )
241         ELSE
242            z_stabfunc = 1._wp / &
243            &                   SQRT( 1._wp - 16._wp * z_sigma )
244         ENDIF
245
246         ! Calculate the T increment
247         z_term1 = ( p_abflux(ji,jj) * ( pmu(ji,jj) + 1._wp)  / &
248         & ( pmu(ji,jj) * pthick(ji,jj) * prho(ji,jj) * rcp ) )
249         
250         
251         z_term2 = -( ( pmu(ji,jj) + 1._wp) * &
252         &                       ( vkarmn * z_fvel * p_fla(ji,jj) ) / &
253         &      ( pthick(ji,jj) * z_stabfunc ) )     
254       
255         t_imp(ji,jj) = ( p_dsst(ji,jj) + p_rdt * z_term1 ) / &
256                        ( 1._wp - p_rdt * z_term2 )
257
258      END_2D
259     
260      END FUNCTION t_imp
261
262END MODULE diu_bulk
Note: See TracBrowser for help on using the repository browser.