- Timestamp:
- 2019-11-29T16:59:07+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcdcy.F90
r10425 r12015 7 7 !! NEMO 2.0 ! 2006-02 (S. Masson, G. Madec) adaptation to NEMO 8 8 !! 3.1 ! 2009-07 (J.M. Molines) adaptation to v3.1 9 !! 4.* ! 2019-10 (L. Brodeau) nothing really new, but the routine 10 !! ! "sbc_dcy_param" has been extracted from old function "sbc_dcy" 11 !! ! => this allows the warm-layer param of COARE3* to know the time 12 !! ! of dawn and dusk even if "ln_dm2dc=.false." (rdawn_dcy & rdusk_dcy 13 !! ! are now public) 9 14 !!---------------------------------------------------------------------- 10 15 … … 22 27 IMPLICIT NONE 23 28 PRIVATE 24 29 25 30 INTEGER, PUBLIC :: nday_qsr !: day when parameters were computed 26 31 27 32 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: raa , rbb , rcc , rab ! diurnal cycle parameters 28 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: rtmd, rdawn, rdusk, rscal ! - - - 29 33 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: rtmd, rscal ! - - - 34 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: rdawn_dcy, rdusk_dcy ! - - - 35 30 36 PUBLIC sbc_dcy ! routine called by sbc 37 PUBLIC sbc_dcy_param ! routine used here and called by warm-layer parameterization (sbcblk_skin_coare*) 31 38 32 39 !!---------------------------------------------------------------------- 33 40 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 34 !! $Id$ 41 !! $Id$ 35 42 !! Software governed by the CeCILL license (see ./LICENSE) 36 43 !!---------------------------------------------------------------------- 37 44 CONTAINS 38 45 39 40 41 42 43 44 & rtmd(jpi,jpj) , rdawn(jpi,jpj) , rdusk(jpi,jpj) , rscal(jpi,jpj) , STAT=sbc_dcy_alloc )45 46 47 48 46 INTEGER FUNCTION sbc_dcy_alloc() 47 !!---------------------------------------------------------------------- 48 !! *** FUNCTION sbc_dcy_alloc *** 49 !!---------------------------------------------------------------------- 50 ALLOCATE( raa (jpi,jpj) , rbb (jpi,jpj) , rcc (jpi,jpj) , rab (jpi,jpj) , & 51 & rtmd(jpi,jpj) , rdawn_dcy(jpi,jpj) , rdusk_dcy(jpi,jpj) , rscal(jpi,jpj) , STAT=sbc_dcy_alloc ) 52 ! 53 CALL mpp_sum ( 'sbcdcy', sbc_dcy_alloc ) 54 IF( sbc_dcy_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sbc_dcy_alloc: failed to allocate arrays' ) 55 END FUNCTION sbc_dcy_alloc 49 56 50 57 … … 60 67 !! 61 68 !! reference : Bernie, DJ, E Guilyardi, G Madec, JM Slingo, and SJ Woolnough, 2007 62 !! Impact of resolving the diurnal cycle in an ocean--atmosphere GCM. 69 !! Impact of resolving the diurnal cycle in an ocean--atmosphere GCM. 63 70 !! Part 1: a diurnally forced OGCM. Climate Dynamics 29:6, 575-590. 64 71 !!---------------------------------------------------------------------- 65 72 LOGICAL , OPTIONAL , INTENT(in) :: l_mask ! use the routine for night mask computation 66 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqsrin ! input daily QSR flux 73 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqsrin ! input daily QSR flux 67 74 REAL(wp), DIMENSION(jpi,jpj) :: zqsrout ! output QSR flux with diurnal cycle 68 75 !! 69 76 INTEGER :: ji, jj ! dummy loop indices 70 77 INTEGER, DIMENSION(jpi,jpj) :: imask_night ! night mask 71 REAL(wp) :: ztwopi, zinvtwopi, zconvrad72 78 REAL(wp) :: zlo, zup, zlousd, zupusd 73 REAL(wp) :: zdsws, zdecrad, ztx, zsin, zcos 74 REAL(wp) :: ztmp, ztmp1, ztmp2, ztest 79 REAL(wp) :: ztmp, ztmp1, ztmp2 75 80 REAL(wp) :: ztmpm, ztmpm1, ztmpm2 76 !---------------------------statement functions------------------------77 REAL(wp) :: fintegral, pt1, pt2, paaa, pbbb, pccc ! dummy statement function arguments78 fintegral( pt1, pt2, paaa, pbbb, pccc ) = &79 & paaa * pt2 + zinvtwopi * pbbb * SIN(pccc + ztwopi * pt2) &80 & - paaa * pt1 - zinvtwopi * pbbb * SIN(pccc + ztwopi * pt1)81 81 !!--------------------------------------------------------------------- 82 82 ! 83 83 ! Initialization 84 84 ! -------------- 85 ztwopi = 2._wp * rpi86 zinvtwopi = 1._wp / ztwopi87 zconvrad = ztwopi / 360._wp88 89 85 ! When are we during the day (from 0 to 1) 90 86 zlo = ( REAL(nsec_day, wp) - 0.5_wp * rdt ) / rday 91 87 zup = zlo + ( REAL(nn_fsbc, wp) * rdt ) / rday 92 ! 93 IF( nday_qsr == -1 ) THEN ! first time step only 88 ! 89 IF( nday_qsr == -1 ) THEN ! first time step only 94 90 IF(lwp) THEN 95 91 WRITE(numout,*) … … 98 94 WRITE(numout,*) 99 95 ENDIF 96 ENDIF 97 98 ! Setting parameters for each new day: 99 CALL sbc_dcy_param() 100 101 !CALL iom_put( "rdusk_dcy", rdusk_dcy(:,:)*tmask(:,:,1) ) !LB 102 !CALL iom_put( "rdawn_dcy", rdawn_dcy(:,:)*tmask(:,:,1) ) !LB 103 !CALL iom_put( "rscal_dcy", rscal(:,:)*tmask(:,:,1) ) !LB 104 105 106 ! 3. update qsr with the diurnal cycle 107 ! ------------------------------------ 108 109 imask_night(:,:) = 0 110 DO jj = 1, jpj 111 DO ji = 1, jpi 112 ztmpm = 0._wp 113 IF( ABS(rab(ji,jj)) < 1. ) THEN ! day duration is less than 24h 114 ! 115 IF( rdawn_dcy(ji,jj) < rdusk_dcy(ji,jj) ) THEN ! day time in one part 116 zlousd = MAX(zlo, rdawn_dcy(ji,jj)) 117 zlousd = MIN(zlousd, zup) 118 zupusd = MIN(zup, rdusk_dcy(ji,jj)) 119 zupusd = MAX(zupusd, zlo) 120 ztmp = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 121 zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 122 ztmpm = zupusd - zlousd 123 IF( ztmpm .EQ. 0 ) imask_night(ji,jj) = 1 124 ! 125 ELSE ! day time in two parts 126 zlousd = MIN(zlo, rdusk_dcy(ji,jj)) 127 zupusd = MIN(zup, rdusk_dcy(ji,jj)) 128 ztmp1 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 129 ztmpm1=zupusd-zlousd 130 zlousd = MAX(zlo, rdawn_dcy(ji,jj)) 131 zupusd = MAX(zup, rdawn_dcy(ji,jj)) 132 ztmp2 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 133 ztmpm2 =zupusd-zlousd 134 ztmp = ztmp1 + ztmp2 135 ztmpm = ztmpm1 + ztmpm2 136 zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 137 IF(ztmpm .EQ. 0.) imask_night(ji,jj) = 1 138 ENDIF 139 ELSE ! 24h light or 24h night 140 ! 141 IF( raa(ji,jj) > rbb(ji,jj) ) THEN ! 24h day 142 ztmp = fintegral(zlo, zup, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 143 zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 144 imask_night(ji,jj) = 0 145 ! 146 ELSE ! No day 147 zqsrout(ji,jj) = 0.0_wp 148 imask_night(ji,jj) = 1 149 ENDIF 150 ENDIF 151 END DO 152 END DO 153 ! 154 IF( PRESENT(l_mask) .AND. l_mask ) THEN 155 zqsrout(:,:) = float(imask_night(:,:)) 156 ENDIF 157 ! 158 END FUNCTION sbc_dcy 159 160 161 SUBROUTINE sbc_dcy_param( ) 162 !! 163 INTEGER :: ji, jj ! dummy loop indices 164 !INTEGER, DIMENSION(jpi,jpj) :: imask_night ! night mask 165 REAL(wp) :: zdsws, zdecrad, ztx, zsin, zcos 166 REAL(wp) :: ztmp, ztest 167 !---------------------------statement functions------------------------ 168 ! 169 IF( nday_qsr == -1 ) THEN ! first time step only 100 170 ! allocate sbcdcy arrays 101 171 IF( sbc_dcy_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_dcy_alloc : unable to allocate arrays' ) 102 172 ! Compute rcc needed to compute the time integral of the diurnal cycle 103 rcc(:,:) = zconvrad * glamt(:,:) - rpi173 rcc(:,:) = rad * glamt(:,:) - rpi 104 174 ! time of midday 105 175 rtmd(:,:) = 0.5_wp - glamt(:,:) / 360._wp … … 107 177 ENDIF 108 178 109 ! If this is a new day, we have to update the dawn, dusk and scaling function 179 ! If this is a new day, we have to update the dawn, dusk and scaling function 110 180 !---------------------- 111 112 ! 2.1 dawn and dusk 113 114 ! nday is the number of days since the beginning of the current month 115 IF( nday_qsr /= nday ) THEN 181 182 ! 2.1 dawn and dusk 183 184 ! nday is the number of days since the beginning of the current month 185 IF( nday_qsr /= nday ) THEN 116 186 ! save the day of the year and the daily mean of qsr 117 nday_qsr = nday 118 ! number of days since the previous winter solstice (supposed to be always 21 December) 187 nday_qsr = nday 188 ! number of days since the previous winter solstice (supposed to be always 21 December) 119 189 zdsws = REAL(11 + nday_year, wp) 120 190 ! declination of the earths orbit 121 zdecrad = (-23.5_wp * zconvrad) * COS( zdsws * ztwopi / REAL(nyear_len(1),wp) )191 zdecrad = (-23.5_wp * rad) * COS( zdsws * 2._wp*rpi / REAL(nyear_len(1),wp) ) 122 192 ! Compute A and B needed to compute the time integral of the diurnal cycle 123 193 … … 125 195 DO jj = 1, jpj 126 196 DO ji = 1, jpi 127 ztmp = zconvrad * gphit(ji,jj)197 ztmp = rad * gphit(ji,jj) 128 198 raa(ji,jj) = SIN( ztmp ) * zsin 129 199 rbb(ji,jj) = COS( ztmp ) * zcos 130 END DO 131 END DO 200 END DO 201 END DO 132 202 ! Compute the time of dawn and dusk 133 203 134 ! rab to test if the day time is equal to 0, less than 24h of full day 204 ! rab to test if the day time is equal to 0, less than 24h of full day 135 205 rab(:,:) = -raa(:,:) / rbb(:,:) 136 206 DO jj = 1, jpj 137 207 DO ji = 1, jpi 138 IF 139 ! When is it night?140 ztx = zinvtwopi* (ACOS(rab(ji,jj)) - rcc(ji,jj))141 ztest = -rbb(ji,jj) * SIN( rcc(ji,jj) + ztwopi * ztx )142 ! is it dawn or dusk?143 IF 144 rdawn (ji,jj) = ztx145 rdusk (ji,jj) = rtmd(ji,jj) + ( rtmd(ji,jj) - rdawn(ji,jj) )208 IF( ABS(rab(ji,jj)) < 1._wp ) THEN ! day duration is less than 24h 209 ! When is it night? 210 ztx = 1._wp/(2._wp*rpi) * (ACOS(rab(ji,jj)) - rcc(ji,jj)) 211 ztest = -rbb(ji,jj) * SIN( rcc(ji,jj) + 2._wp*rpi * ztx ) 212 ! is it dawn or dusk? 213 IF( ztest > 0._wp ) THEN 214 rdawn_dcy(ji,jj) = ztx 215 rdusk_dcy(ji,jj) = rtmd(ji,jj) + ( rtmd(ji,jj) - rdawn_dcy(ji,jj) ) 146 216 ELSE 147 rdusk (ji,jj) = ztx148 rdawn (ji,jj) = rtmd(ji,jj) - ( rdusk(ji,jj) - rtmd(ji,jj) )217 rdusk_dcy(ji,jj) = ztx 218 rdawn_dcy(ji,jj) = rtmd(ji,jj) - ( rdusk_dcy(ji,jj) - rtmd(ji,jj) ) 149 219 ENDIF 150 220 ELSE 151 rdawn (ji,jj) = rtmd(ji,jj) + 0.5_wp152 rdusk (ji,jj) = rdawn(ji,jj)153 ENDIF 154 END DO155 END DO 156 rdawn (:,:) = MOD( (rdawn(:,:) + 1._wp), 1._wp )157 rdusk (:,:) = MOD( (rdusk(:,:) + 1._wp), 1._wp )221 rdawn_dcy(ji,jj) = rtmd(ji,jj) + 0.5_wp 222 rdusk_dcy(ji,jj) = rdawn_dcy(ji,jj) 223 ENDIF 224 END DO 225 END DO 226 rdawn_dcy(:,:) = MOD( (rdawn_dcy(:,:) + 1._wp), 1._wp ) 227 rdusk_dcy(:,:) = MOD( (rdusk_dcy(:,:) + 1._wp), 1._wp ) 158 228 ! 2.2 Compute the scaling function: 159 229 ! S* = the inverse of the time integral of the diurnal cycle from dawn to dusk … … 162 232 DO jj = 1, jpj 163 233 DO ji = 1, jpi 164 IF 234 IF( ABS(rab(ji,jj)) < 1._wp ) THEN ! day duration is less than 24h 165 235 rscal(ji,jj) = 0.0_wp 166 IF ( rdawn(ji,jj) < rdusk(ji,jj) ) THEN ! day time in one part167 IF( (rdusk (ji,jj) - rdawn(ji,jj) ) .ge. 0.001_wp ) THEN168 rscal(ji,jj) = fintegral(rdawn(ji,jj), rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj))169 rscal(ji,jj) = 1._wp / rscal(ji,jj)236 IF( rdawn_dcy(ji,jj) < rdusk_dcy(ji,jj) ) THEN ! day time in one part 237 IF( (rdusk_dcy(ji,jj) - rdawn_dcy(ji,jj) ) .ge. 0.001_wp ) THEN 238 rscal(ji,jj) = fintegral(rdawn_dcy(ji,jj), rdusk_dcy(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 239 rscal(ji,jj) = 1._wp / rscal(ji,jj) 170 240 ENDIF 171 241 ELSE ! day time in two parts 172 IF( (rdusk (ji,jj) + (1._wp - rdawn(ji,jj)) ) .ge. 0.001_wp ) THEN173 rscal(ji,jj) = fintegral(0._wp, rdusk(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) &174 & + fintegral(rdawn(ji,jj), 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj))175 rscal(ji,jj) = 1. / rscal(ji,jj)242 IF( (rdusk_dcy(ji,jj) + (1._wp - rdawn_dcy(ji,jj)) ) .ge. 0.001_wp ) THEN 243 rscal(ji,jj) = fintegral(0._wp, rdusk_dcy(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) & 244 & + fintegral(rdawn_dcy(ji,jj), 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 245 rscal(ji,jj) = 1. / rscal(ji,jj) 176 246 ENDIF 177 247 ENDIF 178 248 ELSE 179 IF 180 rscal(ji,jj) = fintegral(0._wp, 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 249 IF( raa(ji,jj) > rbb(ji,jj) ) THEN ! 24h day 250 rscal(ji,jj) = fintegral(0._wp, 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 181 251 rscal(ji,jj) = 1._wp / rscal(ji,jj) 182 252 ELSE ! No day … … 184 254 ENDIF 185 255 ENDIF 186 END DO 187 END DO 256 END DO 257 END DO 188 258 ! 189 259 ztmp = rday / ( rdt * REAL(nn_fsbc, wp) ) 190 260 rscal(:,:) = rscal(:,:) * ztmp 191 261 ! 192 ENDIF 193 ! 3. update qsr with the diurnal cycle 194 ! ------------------------------------ 195 196 imask_night(:,:) = 0 197 DO jj = 1, jpj 198 DO ji = 1, jpi 199 ztmpm = 0._wp 200 IF( ABS(rab(ji,jj)) < 1. ) THEN ! day duration is less than 24h 201 ! 202 IF( rdawn(ji,jj) < rdusk(ji,jj) ) THEN ! day time in one part 203 zlousd = MAX(zlo, rdawn(ji,jj)) 204 zlousd = MIN(zlousd, zup) 205 zupusd = MIN(zup, rdusk(ji,jj)) 206 zupusd = MAX(zupusd, zlo) 207 ztmp = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 208 zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 209 ztmpm = zupusd - zlousd 210 IF ( ztmpm .EQ. 0 ) imask_night(ji,jj) = 1 211 ! 212 ELSE ! day time in two parts 213 zlousd = MIN(zlo, rdusk(ji,jj)) 214 zupusd = MIN(zup, rdusk(ji,jj)) 215 ztmp1 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 216 ztmpm1=zupusd-zlousd 217 zlousd = MAX(zlo, rdawn(ji,jj)) 218 zupusd = MAX(zup, rdawn(ji,jj)) 219 ztmp2 = fintegral(zlousd, zupusd, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 220 ztmpm2 =zupusd-zlousd 221 ztmp = ztmp1 + ztmp2 222 ztmpm = ztmpm1 + ztmpm2 223 zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 224 IF (ztmpm .EQ. 0.) imask_night(ji,jj) = 1 225 ENDIF 226 ELSE ! 24h light or 24h night 227 ! 228 IF( raa(ji,jj) > rbb(ji,jj) ) THEN ! 24h day 229 ztmp = fintegral(zlo, zup, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 230 zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 231 imask_night(ji,jj) = 0 232 ! 233 ELSE ! No day 234 zqsrout(ji,jj) = 0.0_wp 235 imask_night(ji,jj) = 1 236 ENDIF 237 ENDIF 238 END DO 239 END DO 240 ! 241 IF( PRESENT(l_mask) .AND. l_mask ) THEN 242 zqsrout(:,:) = float(imask_night(:,:)) 243 ENDIF 244 ! 245 END FUNCTION sbc_dcy 262 ENDIF !IF( nday_qsr /= nday ) 263 ! 264 END SUBROUTINE sbc_dcy_param 265 266 267 FUNCTION fintegral( pt1, pt2, paaa, pbbb, pccc ) 268 REAL(wp), INTENT(in) :: pt1, pt2, paaa, pbbb, pccc 269 REAL(wp) :: fintegral 270 fintegral = paaa * pt2 + 1._wp/(2._wp*rpi) * pbbb * SIN(pccc + 2._wp*rpi*pt2) & 271 & - paaa * pt1 - 1._wp/(2._wp*rpi) * pbbb * SIN(pccc + 2._wp*rpi*pt1) 272 END FUNCTION fintegral 246 273 247 274 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
