Changeset 840
- Timestamp:
- 2008-03-11T14:51:35+01:00 (16 years ago)
- Location:
- branches/dev_001_SBC/NEMO/OPA_SRC/SBC
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/dev_001_SBC/NEMO/OPA_SRC/SBC/sbcblk_clio.F90
r710 r840 1 1 MODULE sbcblk_clio 2 2 !!====================================================================== 3 !! 4 !! Ocean forcing: momentum, heat and freshwater flux formulation3 !! *** MODULE sbcblk_clio *** 4 !! Ocean forcing: bulk thermohaline forcing of the ocean (or ice) 5 5 !!===================================================================== 6 !! History : 8.0 ! 01-04 (Louvain-La-Neuve) Original code 7 !! 8.5 ! 02-09 (C. Ethe , G. Madec ) F90: Free form and module 8 !! 9.0 ! 06-06 (G. Madec) surface module 6 !! History : OPA ! 1997-06 (Louvain-La-Neuve) Original code 7 !! ! 2001-04 (C. Ethe) add flx_blk_declin 8 !! NEMO 2.0 ! 2002-08 (C. Ethe, G. Madec) F90: Free form and module 9 !! 3.0 ! 2008-03 (C. Talandier, G. Madec) surface module + LIM3 9 10 !!---------------------------------------------------------------------- 10 11 !!---------------------------------------------------------------------- 12 !! sbc_blk_clio : bulk formulation as ocean surface boundary condition 13 !! (forced mode, CORE bulk formulea) 14 !! blk_oce_clio : ocean: computes momentum, heat and freshwater fluxes 15 !! blk_ice_clio : ice : computes momentum, heat and freshwater fluxes 16 !!---------------------------------------------------------------------- 17 !! flx_blk_declin : Computation of the solar declination 11 !! sbc_blk_clio : CLIO bulk formulation: read and update required input fields 12 !! blk_clio_oce : ocean CLIO bulk formulea: compute momentum, heat and freswater fluxes for the ocean 13 !! blk_ice_clio : ice CLIO bulk formulea: compute momentum, heat and freswater fluxes for the sea-ice 14 !! blk_clio_qsr_oce : shortwave radiation for ocean computed from the cloud cover 15 !! blk_clio_qsr_ice : shortwave radiation for ice computed from the cloud cover 16 !! flx_blk_declin : solar declinaison 18 17 !!---------------------------------------------------------------------- 19 18 USE oce ! ocean dynamics and tracers 20 19 USE dom_oce ! ocean space and time domain 21 USE cpl_oce ! ???22 20 USE phycst ! physical constants 23 USE daymod 24 21 USE daymod ! calendar 22 USE ocfzpt ! ocean freezing point 23 USE fldread ! read input fields 25 24 USE sbc_oce ! Surface boundary condition: ocean fields 26 USE sbc_ice ! Surface boundary condition: ocean fields 27 28 USE fldread ! read input fields 29 30 USE ocfzpt ! ??? 31 32 USE iom 33 USE in_out_manager 34 USE lbclnk 25 USE iom ! I/O manager library 26 USE in_out_manager ! I/O manager 27 USE lib_mpp ! distribued memory computing library 28 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 29 30 USE albedo 35 31 USE prtctl ! Print control 36 32 #if defined key_lim3 33 USE par_ice 34 #elif defined key_lim2 35 USE ice_2 36 #endif 37 37 IMPLICIT NONE 38 38 PRIVATE 39 39 40 PUBLIC sbc_blk_clio ! routine called in sbcmod 41 PUBLIC blk_ice_clio ! routine called in sbcice_lim module 42 43 INTEGER , PARAMETER :: & 44 jpfld = 7, & ! number of files to read 45 jp_utau = 1, & ! index of wind stress (i-component) (m/s) at U-point 46 jp_vtau = 2, & ! index of wind stress (j-component) (m/s) at V-point 47 jp_wspd = 3, & ! index of XXm wind module (m/s) at T-point 48 jp_tair = 4, & ! index of 2m air temperature (Celcius) 49 jp_humi = 5, & ! index of specific humidity ( % ) 50 jp_cldc = 6, & ! Cloud cover ( % ) 51 jp_prec = 7 ! index of total precipitation (kg/m2/s) 52 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf ! structure of input fields (informations on files, fields read) 53 54 !! * CLIO bulk parameters 40 PUBLIC sbc_blk_clio ! routine called by flx.F90 41 PUBLIC blk_ice_clio ! routine called by flx.F90 42 43 INTEGER , PARAMETER :: jpfld = 7 ! maximum number of files to read 44 INTEGER , PARAMETER :: jp_wndi = 1 ! index of 10m wind velocity (i-component) (m/s) at U-point 45 INTEGER , PARAMETER :: jp_wndj = 2 ! index of 10m wind velocity (j-component) (m/s) at V-point 46 INTEGER , PARAMETER :: jp_wndm = 3 ! index of 10m wind module (m/s) at T-point 47 INTEGER , PARAMETER :: jp_humi = 4 ! index of specific humidity ( % ) 48 INTEGER , PARAMETER :: jp_ccov = 5 ! index of cloud cover ( % ) 49 INTEGER , PARAMETER :: jp_tair = 6 ! index of 10m air temperature (Kelvin) 50 INTEGER , PARAMETER :: jp_prec = 7 ! index of total precipitation (rain+snow) (Kg/m2/s) 51 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf ! structure of input fields (file informations, fields read) 52 53 54 !! 55 !!!!gm to be moved 56 INTEGER, PARAMETER :: jpl = 1 ! number of layer in the ice 57 !!!!gm to be moved 58 59 55 60 INTEGER, PARAMETER :: jpintsr = 24 ! number of time step between sunrise and sunset 56 ! ! uses for heat flux computation 57 LOGICAL :: lbulk_init = .TRUE. ! flag, bulk initialization done or not) 58 59 REAL(wp), DIMENSION(jpi,jpj) :: stauc ! cloud optical depth 60 REAL(wp), DIMENSION(jpi,jpj) :: sbudyko ! ??? 61 62 !! * constants for bulk computation (flx_blk) 63 REAL(wp), DIMENSION(19) :: budyko ! BUDYKO's coefficient 64 ! ! BUDYKO's coefficient (cloudiness effect on LW radiation): 65 DATA budyko / 1.00, 0.98, 0.95, 0.92, 0.89, 0.86, 0.83, 0.80, 0.78, 0.75, & 61 ! ! uses for heat flux computation 62 LOGICAL :: lbulk_init = .TRUE. ! flag, bulk initialization done or not) 63 64 REAL(wp) :: cai = 1.40e-3 ! best estimate of atm drag in order to get correct FS export in ORCA2-LIM 65 REAL(wp) :: cao = 1.00e-3 ! chosen by default ==> should depends on many things... !!gmto be updated 66 67 REAL(wp) :: yearday !: number of days per year 68 REAL(wp) :: rdtbs2 !: number of days per year 69 70 REAL(wp), DIMENSION(19) :: budyko ! BUDYKO's coefficient (cloudiness effect on LW radiation) 71 DATA budyko / 1.00, 0.98, 0.95, 0.92, 0.89, 0.86, 0.83, 0.80, 0.78, 0.75, & 66 72 & 0.72, 0.69, 0.67, 0.64, 0.61, 0.58, 0.56, 0.53, 0.50 / 67 REAL(wp), DIMENSION(20) :: tauco ! cloud optical depth coefficient 68 ! ! Cloud optical depth coefficient 73 REAL(wp), DIMENSION(20) :: tauco ! cloud optical depth coefficient 69 74 DATA tauco / 6.6, 6.6, 7.0, 7.2, 7.1, 6.8, 6.5, 6.6, 7.1, 7.6, & 70 75 & 6.6, 6.1, 5.6, 5.5, 5.8, 5.8, 5.6, 5.6, 5.6, 5.6 / 71 72 REAL(wp) :: zeps = 1.e-20 , & ! constant values 73 & zeps0 = 1.e-13 , & 74 & zeps1 = 1.e-06 , & 75 & zzero = 0.e0 , & 76 & zone = 1.0 77 78 REAL(wp) :: yearday ! 79 76 !! 77 REAL(wp), DIMENSION(jpi,jpj) :: sbudyko ! cloudiness effect on LW radiation 78 REAL(wp), DIMENSION(jpi,jpj) :: stauc ! cloud optical depth 79 80 REAL(wp) :: zeps = 1.e-20 ! constant values 81 REAL(wp) :: zeps0 = 1.e-13 82 83 # include "vectopt_loop_substitute.h90" 80 84 !!---------------------------------------------------------------------- 81 !! OPA 9.0 , LOCEAN-IPSL (2006)82 !! $ Header: $85 !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 86 !! $Id:$ 83 87 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 84 88 !!---------------------------------------------------------------------- … … 86 90 CONTAINS 87 91 92 88 93 SUBROUTINE sbc_blk_clio( kt ) 89 94 !!--------------------------------------------------------------------- 90 !! *** ROUTINE sbc_blk_c lio***91 !! 95 !! *** ROUTINE sbc_blk_core *** 96 !! 92 97 !! ** Purpose : provide at each time step the surface ocean fluxes 93 !! (momentum, heat, freshwater and runoff) 98 !! (momentum, heat, freshwater and runoff) 94 99 !! 95 100 !! ** Method : READ each fluxes in NetCDF files … … 118 123 INTEGER, INTENT( in ) :: kt ! ocean time step 119 124 !! 120 INTEGER :: jf 125 INTEGER :: jf ! dummy indices 121 126 INTEGER :: ierror ! return error code 122 127 !! 123 CHARACTER(len=100) :: cn_dir ! Root directory for location of clio files 124 TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read 125 TYPE(FLD_N) :: sn_utau, sn_vtau, sn_wspd, sn_tair, & ! informations about the fields to be read 126 & sn_humi, sn_cldc, sn_prec 127 NAMELIST/namsbc_clio/ cn_dir, sn_utau, sn_vtau, sn_wspd, sn_tair, & 128 & sn_humi, sn_cldc, sn_prec 128 CHARACTER(len=100) :: cn_dir ! Root directory for location of CLIO files 129 TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read 130 TYPE(FLD_N) :: sn_wndi, sn_wndj, sn_wndm, sn_tair ! informations about the fields to be read 131 TYPE(FLD_N) :: sn_humi, sn_ccov, sn_prec ! " " 132 !! 133 NAMELIST/namsbc_core/ cn_dir, sn_wndi, sn_wndj, sn_wndm, sn_humi, & 134 & sn_ccov, sn_tair, sn_prec 129 135 !!--------------------------------------------------------------------- 130 136 … … 134 140 ! set file information (default values) 135 141 cn_dir = './' ! directory in which the model is executed 142 136 143 ! (NB: frequency positive => hours, negative => months) 137 ! ! file! frequency ! variable ! time intep ! clim ! starting !138 ! ! name! (hours) ! name ! (T/F) ! (0/1) ! record !139 sn_ utau = FLD_N( 'utau' , 24. , 'utau' , .FALSE. , 0 , 0 )140 sn_ vtau = FLD_N( 'vtau' , 24. , 'vtau' , .FALSE. , 0 , 0 )141 sn_w spd = FLD_N( 'wspd' , 24. , 'wspd' , .FALSE. , 0 , 0 )142 sn_tair = FLD_N( 'tair ' , 24. , 'Tair' , .FALSE. , 0 , 0 )143 sn_humi = FLD_N( 'humi ' , -12. , 'humi' , .FALSE. , 0 , 0 )144 sn_c ldc = FLD_N( 'cloud', -12. , 'cloud' , .FALSE. , 0 , 0 )145 sn_prec = FLD_N( ' rain' , -12. , 'precip' , .FALSE. , 0 , 0 )146 147 REWIND ( numnam ) ! ... read in namlist namsbc_clio148 READ ( numnam, namsbc_clio)144 ! ! file ! frequency ! variable ! time intep ! clim ! starting ! 145 ! ! name ! (hours) ! name ! (T/F) ! (0/1) ! record ! 146 sn_wndi = FLD_N( 'uwnd10m' , 24. , 'u_10' , .true. , 0 , 0 ) 147 sn_wndj = FLD_N( 'vwnd10m' , 24. , 'v_10' , .true. , 0 , 0 ) 148 sn_wndm = FLD_N( 'mwnd10m' , 24. , 'm_10' , .true. , 0 , 0 ) 149 sn_tair = FLD_N( 'tair10m' , 24. , 't_10' , .FALSE. , 0 , 0 ) 150 sn_humi = FLD_N( 'humi10m' , 24. , 'q_10' , .FALSE. , 0 , 0 ) 151 sn_ccov = FLD_N( 'ccover' , -12. , 'cloud' , .TRUE. , 0 , 0 ) 152 sn_prec = FLD_N( 'precip' , -12. , 'precip' , .TRUE. , 0 , 0 ) 153 154 REWIND( numnam ) ! ... read in namlist namsbc_core 155 READ ( numnam, namsbc_core ) 149 156 150 157 ! store namelist information in an array 151 slf_i(jp_utau) = sn_utau ; slf_i(jp_vtau) = sn_vtau 152 slf_i(jp_wspd) = sn_wspd ; slf_i(jp_tair) = sn_tair 153 slf_i(jp_humi) = sn_humi ; slf_i(jp_cldc) = sn_cldc 154 slf_i(jp_prec) = sn_prec 155 158 slf_i(jp_wndi) = sn_wndi ; slf_i(jp_wndj) = sn_wndj ; slf_i(jp_wndm) = sn_wndm 159 slf_i(jp_tair) = sn_tair ; slf_i(jp_humi) = sn_humi 160 slf_i(jp_ccov) = sn_ccov ; slf_i(jp_prec) = sn_prec 161 156 162 ! set sf structure 157 163 ALLOCATE( sf(jpfld), STAT=ierror ) 158 164 IF( ierror > 0 ) THEN 159 CALL ctl_stop( 'sbc_blk_clio: unable to allocate sf _sststructure' ) ; RETURN165 CALL ctl_stop( 'sbc_blk_clio: unable to allocate sf structure' ) ; RETURN 160 166 ENDIF 161 167 ! 162 168 DO jf = 1, jpfld 163 169 WRITE(sf(jf)%clrootname,'(a,a)' ) TRIM( cn_dir ), TRIM( slf_i(jf)%clname ) … … 170 176 171 177 IF(lwp) THEN ! control print 172 WRITE(numout,*) 173 WRITE(numout,*) 'sbc_blk_clio : CLIO bulk formulation for ocean surface boundary condition'178 WRITE(numout,*) 179 WRITE(numout,*) 'sbc_blk_clio : flux formulattion for ocean surface boundary condition' 174 180 WRITE(numout,*) '~~~~~~~~~~~~ ' 175 WRITE(numout,*) ' namsbc_c lioNamelist'181 WRITE(numout,*) ' namsbc_core Namelist' 176 182 WRITE(numout,*) ' list of files and frequency (>0: in hours ; <0 in months)' 177 183 DO jf = 1, jpfld 178 WRITE(numout,*) ' root filename: ' , trim( sf(jf)%clrootname ), &179 & ' variable name: ' , trim( sf(jf)%clvar )184 WRITE(numout,*) ' file root name: ' , TRIM( sf(jf)%clrootname ), & 185 & ' variable name: ' , TRIM( sf(jf)%clvar ) 180 186 WRITE(numout,*) ' frequency: ' , sf(jf)%freqh , & 181 187 & ' time interp: ' , sf(jf)%ln_tint , & … … 186 192 ! 187 193 ENDIF 188 189 CALL fld_read( kt, nn_fsbc, sf ) ! Read input fields and provides the 190 ! ! input fields at the current time-step 191 !CDIR COLLAPSE 192 utau(:,:) = sf(jp_utau)%fnow(:,:) ! set surface ocean stresses directly 193 !CDIR COLLAPSE 194 vtau(:,:) = sf(jp_vtau)%fnow(:,:) ! from the input values 195 196 CALL blk_oce_clio( ) ! set the ocean surface heat and freshwater fluxes 197 ! ! using CLIO bulk formulea 198 199 ! temporary staff : set fluxes to zero.... 200 qns (:,:)= 0.e0 201 qsr (:,:)= 0.e0 202 emp (:,:)= 0.e0 203 emps(:,:)= 0.e0 204 205 ! control print (if less than 100 time-step asked) 206 !!! IF( nitend-nit000 <= 100 .AND. lwp ) THEN 207 IF( kt == nit000 .AND. lwp ) THEN 208 WRITE(numout,*) 209 WRITE(numout,*) ' CLIO bulk fields at nit000' 210 DO jf = 1, jpfld 211 WRITE(numout,*) 212 WRITE(numout,*) ' day: ', ndastp , TRIM(sf(jf)%clvar) 213 CALL prihre( sf(jf)%fnow, jpi, jpj, 1, jpi, 20, 1, jpj, 10, 1., numout ) 214 END DO 215 CALL FLUSH(numout) 194 ! ! ====================== ! 195 ! ! At each time-step ! 196 ! ! ====================== ! 197 ! 198 CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step 199 ! 200 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 201 CALL blk_oce_clio( sst_m, ssu_m, ssv_m ) ! compute the surface ocean fluxes using CLIO bulk formulea 202 ENDIF ! 203 ! 204 END SUBROUTINE sbc_blk_clio 205 206 207 SUBROUTINE blk_oce_clio( pst, pu, pv ) 208 !!--------------------------------------------------------------------------- 209 !! *** ROUTINE blk_oce_clio *** 210 !! 211 !! ** Purpose : Compute momentum, heat and freshwater fluxes at ocean surface 212 !! using CLIO bulk formulea 213 !! 214 !! ** Method : The flux of heat at the ocean surfaces are derived 215 !! from semi-empirical ( or bulk ) formulae which relate the flux to 216 !! the properties of the surface and of the lower atmosphere. Here, we 217 !! follow the work of Oberhuber, 1988 218 !! - momentum flux (stresses) directly read in files at U- and V-points 219 !! - compute ocean and ice albedos (call flx_blk_albedo) 220 !! - compute shortwave radiation for ocean (call blk_clio_qsr_oce) 221 !! - compute long-wave radiation for the ocean 222 !! - compute the turbulent heat fluxes over the ocean 223 !! - deduce the evaporation over the ocean 224 !! ** Action : Fluxes over the ocean: 225 !! - utau, vtau i- and j-component of the wind stress 226 !! - qns, qsr non-slor and solar heat flux 227 !! - emp, emps evaporation minus precipitation 228 !!---------------------------------------------------------------------- 229 REAL(wp), INTENT(in), DIMENSION(jpi,jpj) :: pst ! surface temperature [Celcius] 230 REAL(wp), INTENT(in), DIMENSION(jpi,jpj) :: pu ! surface current at U-point (i-component) [m/s] 231 REAL(wp), INTENT(in), DIMENSION(jpi,jpj) :: pv ! surface current at V-point (j-component) [m/s] 232 !! 233 INTEGER :: ji, jj ! dummy loop indices 234 !! 235 REAL(wp) :: zrhova, zcsho, zcleo, zcldeff ! temporary scalars 236 REAL(wp) :: zqsato, zdteta, zdeltaq, ztvmoy, zobouks ! - - 237 REAL(wp) :: zpsims, zpsihs, zpsils, zobouku, zxins, zpsimu ! - - 238 REAL(wp) :: zpsihu, zpsilu, zstab,zpsim, zpsih, zpsil ! - - 239 REAL(wp) :: zvatmg, zcmn, zchn, zcln, zcmcmn, zdenum ! - - 240 REAL(wp) :: zdtetar, ztvmoyr, zlxins, zchcm, zclcm ! - - 241 REAL(wp) :: zmt1, zmt2, zmt3, ztatm3, ztamr, ztaevbk ! - - 242 REAL(wp) :: zcoef, zsst, ztatm, zcco1, zpatm ! - - 243 REAL(wp) :: zrhoa, zev, zes, zeso, zqatm, zevsqr ! - - 244 !! 245 REAL(wp), DIMENSION(jpi,jpj) :: zqlw ! long-wave heat flux over ocean 246 REAL(wp), DIMENSION(jpi,jpj) :: zqla ! latent heat flux over ocean 247 REAL(wp), DIMENSION(jpi,jpj) :: zqsb ! sensible heat flux over ocean 248 !!--------------------------------------------------------------------- 249 250 zpatm = 101000. ! atmospheric pressure (assumed constant here) 251 252 !------------------------------------! 253 ! momentum fluxes (utau, vtau ) ! 254 !------------------------------------! 255 256 ! Change from wind speed to wind stress over OCEAN (cao is used) 257 ! note: 1.3 = air density) 258 !!gm the module is wrong as U and V points are not the same 259 !!gm the surface currents are not taken into account... 260 !CDIR COLLAPSE 261 DO jj = 1 , jpj 262 DO ji = 1, jpi 263 zcoef = 1.3 * cao * SQRT( sf(jp_wndi)%fnow(ji,jj)*sf(jp_wndi)%fnow(ji,jj) & 264 & + sf(jp_wndj)%fnow(ji,jj)*sf(jp_wndj)%fnow(ji,jj) ) 265 utau(ji,jj) = 1.3 * cao * zcoef * sf(jp_wndi)%fnow(ji,jj) 266 vtau(ji,jj) = 1.3 * cao * zcoef * sf(jp_wndj)%fnow(ji,jj) 267 END DO 268 END DO 269 !!gm better coding: use the wind module and averaging to U and V points 270 ! zcoef = 1.3 * cao * 0.5 271 ! DO jj = 1 , jpjm1 272 ! DO ji = 1, fs_jpim1 273 ! utau(ji,jj) = zcoef * ( sf(jp_wndm)%fnow(ji+1,jj) + sf(jp_wndm)%fnow(ji,jj) ) * sf(jp_wndi)%fnow(ji,jj) 274 ! vtau(ji,jj) = zcoef * ( sf(jp_wndm)%fnow(ji,jj+1) + sf(jp_wndm)%fnow(ji,jj) ) * sf(jp_wndj)%fnow(ji,jj) 275 ! END DO 276 ! END DO 277 ! CALL lbc_lnk( utau(:,:), 'U', -1. ) 278 ! CALL lbc_lnk( vtau(:,:), 'V', -1. ) 279 !!gm end coding 280 281 282 !------------------------------------------------! 283 ! Shortwave radiation for ocean and snow/ice ! 284 !------------------------------------------------! 285 286 CALL blk_clio_qsr_oce( qsr ) 287 288 289 !------------------------! 290 ! Other ocean fluxes ! 291 !------------------------! 292 !CDIR NOVERRCHK 293 !CDIR COLLAPSE 294 DO jj = 1, jpj 295 !CDIR NOVERRCHK 296 DO ji = 1, jpi 297 ! 298 zsst = pst(ji,jj) + rt0 ! converte Celcius to Kelvin the SST and Tair 299 ztatm = sf(jp_tair)%fnow(ji,jj) + rt0 ! and set minimum value far above 0 K (=rt0 over land) 300 zcco1 = 1.0 - sf(jp_ccov)%fnow(ji,jj) ! fraction of clear sky ( 1 - cloud cover) 301 zrhoa = zpatm / ( 287.04 * ztatm ) ! air density (equation of state for dry air) 302 ztamr = ztatm - rtt ! Saturation water vapour 303 zmt1 = SIGN( 17.269, ztamr ) ! || 304 zmt2 = SIGN( 21.875, ztamr ) ! \ / 305 zmt3 = SIGN( 28.200, -ztamr ) ! \/ 306 zes = 611.0 * EXP( ABS( ztamr ) * MIN ( zmt1, zmt2 ) / ( ztatm - 35.86 + MAX( 0.e0, zmt3 ) ) ) 307 zev = sf(jp_humi)%fnow(ji,jj) * zes ! vapour pressure 308 zevsqr = SQRT( zev * 0.01 ) ! square-root of vapour pressure 309 zqatm = 0.622 * zev / ( zpatm - 0.378 * zev ) ! specific humidity 310 311 !--------------------------------------! 312 ! long-wave radiation over the ocean ! ( Berliand 1952 ; all latitudes ) 313 !--------------------------------------! 314 ztatm3 = ztatm * ztatm * ztatm 315 zcldeff = 1.0 - sbudyko(ji,jj) * sf(jp_ccov)%fnow(ji,jj) * sf(jp_ccov)%fnow(ji,jj) 316 ztaevbk = ztatm * ztatm3 * zcldeff * ( 0.39 - 0.05 * zevsqr ) 317 ! 318 zqlw(ji,jj) = - emic * stefan * ( ztaevbk + 4. * ztatm3 * ( zsst - ztatm ) ) 319 320 !-------------------------------------------------- 321 ! Latent and sensible heat fluxes over the ocean 322 !-------------------------------------------------- 323 ! ! vapour pressure at saturation of ocean 324 zeso = 611.0 * EXP ( 17.2693884 * ( zsst - rtt ) * tmask(ji,jj,1) / ( zsst - 35.86 ) ) 325 326 zqsato = ( 0.622 * zeso ) / ( zpatm - 0.378 * zeso ) ! humidity close to the ocean surface (at saturation) 327 328 ! Drag coefficients from Large and Pond (1981,1982) 329 ! ! Stability parameters 330 zdteta = zsst - ztatm 331 zdeltaq = zqatm - zqsato 332 ztvmoy = ztatm * ( 1. + 2.2e-3 * ztatm * zqatm ) 333 zdenum = MAX( sf(jp_wndm)%fnow(ji,jj) * sf(jp_wndm)%fnow(ji,jj) * ztvmoy, zeps ) 334 zdtetar = zdteta / zdenum 335 ztvmoyr = ztvmoy * ztvmoy * zdeltaq / zdenum 336 ! ! case of stable atmospheric conditions 337 zobouks = -70.0 * 10. * ( zdtetar + 3.2e-3 * ztvmoyr ) 338 zobouks = MAX( 0.e0, zobouks ) 339 zpsims = -7.0 * zobouks 340 zpsihs = zpsims 341 zpsils = zpsims 342 ! ! case of unstable atmospheric conditions 343 zobouku = MIN( 0.e0, -100.0 * 10.0 * ( zdtetar + 2.2e-3 * ztvmoyr ) ) 344 zxins = ( 1. - 16. * zobouku )**0.25 345 zlxins = LOG( ( 1. + zxins * zxins ) / 2. ) 346 zpsimu = 2. * LOG( ( 1 + zxins ) * 0.5 ) + zlxins - 2. * ATAN( zxins ) + rpi * 0.5 347 zpsihu = 2. * zlxins 348 zpsilu = zpsihu 349 ! ! intermediate values 350 zstab = MAX( 0.e0, SIGN( 1.e0, zdteta ) ) 351 zpsim = zstab * zpsimu + ( 1.0 - zstab ) * zpsims 352 zpsih = zstab * zpsihu + ( 1.0 - zstab ) * zpsihs 353 zpsil = zpsih 354 355 zvatmg = MAX( 0.032 * 1.5e-3 * sf(jp_wndm)%fnow(ji,jj) * sf(jp_wndm)%fnow(ji,jj) / grav, zeps ) 356 zcmn = vkarmn / LOG ( 10. / zvatmg ) 357 zchn = 0.0327 * zcmn 358 zcln = 0.0346 * zcmn 359 zcmcmn = 1. / ( 1. - zcmn * zpsim / vkarmn ) 360 zchcm = zcmcmn / ( 1. - zchn * zpsih / ( vkarmn * zcmn ) ) 361 zclcm = zchcm 362 ! ! transfert coef. (Large and Pond 1981,1982) 363 zcsho = zchn * zchcm 364 zcleo = zcln * zclcm 365 366 zrhova = zrhoa * sf(jp_wndm)%fnow(ji,jj) 367 368 ! sensible heat flux 369 zqsb(ji,jj) = zrhova * zcsho * 1004.0 * ( zsst - ztatm ) 370 371 ! latent heat flux (bounded by zero) 372 zqla(ji,jj) = MAX( 0.e0, zrhova * zcleo * 2.5e+06 * ( zqsato - zqatm ) ) 373 ! 374 END DO 375 END DO 376 377 ! ----------------------------------------------------------------------------- ! 378 ! III Total FLUXES ! 379 ! ----------------------------------------------------------------------------- ! 380 !CDIR COLLAPSE 381 qns (:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) ! Downward Non Solar flux 382 !CDIR COLLAPSE 383 emp (:,:) = zqla(:,:) / cevap - sf(jp_prec)%fnow(:,:) / rday * tmask(:,:,1) 384 !CDIR COLLAPSE 385 emps(:,:) = emp(:,:) 386 ! 387 END SUBROUTINE blk_oce_clio 388 389 390 SUBROUTINE blk_ice_clio( pst , pui , pvi , palb_cs, palb_os , & 391 & p_taui, p_tauj, p_qns , p_qsr, & 392 & p_qla , p_dqns, p_dqla, & 393 & p_tpr , p_spr , & 394 & p_fr1 , p_fr2 , cd_grid ) 395 !!--------------------------------------------------------------------------- 396 !! *** ROUTINE blk_ice_clio *** 397 !! 398 !! ** Purpose : Computation of the heat fluxes at ocean and snow/ice 399 !! surface the solar heat at ocean and snow/ice surfaces and the 400 !! sensitivity of total heat fluxes to the SST variations 401 !! 402 !! ** Method : The flux of heat at the ice and ocean surfaces are derived 403 !! from semi-empirical ( or bulk ) formulae which relate the flux to 404 !! the properties of the surface and of the lower atmosphere. Here, we 405 !! follow the work of Oberhuber, 1988 406 !! 407 !! ** Action : call flx_blk_albedo to compute ocean and ice albedo 408 !! computation of snow precipitation 409 !! computation of solar flux at the ocean and ice surfaces 410 !! computation of the long-wave radiation for the ocean and sea/ice 411 !! computation of turbulent heat fluxes over water and ice 412 !! computation of evaporation over water 413 !! computation of total heat fluxes sensitivity over ice (dQ/dT) 414 !! computation of latent heat flux sensitivity over ice (dQla/dT) 415 !! 416 !!---------------------------------------------------------------------- 417 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpl) :: pst ! ice surface temperature [Kelvin] 418 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pui ! ice surface velocity (i-component, I-point) [m/s] 419 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pvi ! ice surface velocity (j-component, I-point) [m/s] 420 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpl) :: palb_cs ! ice albedo (clear sky) (alb_ice_cs) [%] 421 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpl) :: palb_os ! ice albedo (overcast sky) (alb_ice_cs) [%] 422 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: p_taui ! surface ice stress at I-point (i-component) [N/m2] 423 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: p_tauj ! surface ice stress at I-point (j-component) [N/m2] 424 REAL(wp), INTENT( out), DIMENSION(jpi,jpj,jpl) :: p_qns ! non solar heat flux over ice (T-point) [W/m2] 425 REAL(wp), INTENT( out), DIMENSION(jpi,jpj,jpl) :: p_qsr ! solar heat flux over ice (T-point) [W/m2] 426 REAL(wp), INTENT( out), DIMENSION(jpi,jpj,jpl) :: p_qla ! latent heat flux over ice (T-point) [W/m2] 427 REAL(wp), INTENT( out), DIMENSION(jpi,jpj,jpl) :: p_dqns ! non solar heat sensistivity (T-point) [W/m2] 428 REAL(wp), INTENT( out), DIMENSION(jpi,jpj,jpl) :: p_dqla ! latent heat sensistivity (T-point) [W/m2] 429 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: p_tpr ! total precipitation (T-point) [Kg/m2/s] 430 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: p_spr ! solid precipitation (T-point) [Kg/m2/s] 431 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: p_fr1 ! 1sr fraction of qsr penetration in ice [%] 432 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: p_fr2 ! 2nd fraction of qsr penetration in ice [%] 433 CHARACTER(len=1), INTENT(in ) :: cd_grid ! type of sea-ice grid ("C" or "B" grid) 434 !! 435 INTEGER :: ji, jj, jl ! dummy loop indices 436 !! 437 REAL(wp) :: zcoef, zmt1, zmt2, zmt3, ztatm3 ! temporary scalars 438 REAL(wp) :: ztaevbk, zind1, zind2, zind3, ztamr ! - - 439 REAL(wp) :: zesi, zqsati, zdesidt ! - - 440 REAL(wp) :: zdqla, zcldeff, zev, zes, zpatm, zrhova ! - - 441 REAL(wp) :: zcshi, zclei, zrhovaclei, zrhovacshi ! - - 442 REAL(wp) :: ztice3, zticemb, zticemb2, zdqlw, zdqsb ! - - 443 !! 444 REAL(wp), DIMENSION(jpi,jpj) :: ztatm ! Tair in Kelvin 445 REAL(wp), DIMENSION(jpi,jpj) :: zqatm ! specific humidity 446 REAL(wp), DIMENSION(jpi,jpj) :: zevsqr ! vapour pressure square-root 447 REAL(wp), DIMENSION(jpi,jpj) :: zrhoa ! air density 448 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qlw, z_qsb 449 !!--------------------------------------------------------------------- 450 451 zpatm = 101000. ! atmospheric pressure (assumed constant here) 452 453 !------------------------------------! 454 ! momentum fluxes (utau, vtau ) ! 455 !------------------------------------! 456 457 SELECT CASE( cd_grid ) 458 CASE( 'C' ) ! C-grid ice dynamics 459 ! Change from wind speed to wind stress over OCEAN (cao is used) 460 zcoef = cai / cao 461 !CDIR COLLAPSE 462 DO jj = 1 , jpj 463 DO ji = 1, jpi 464 p_taui(ji,jj) = zcoef * utau(ji,jj) 465 p_tauj(ji,jj) = zcoef * vtau(ji,jj) 466 END DO 467 END DO 468 CASE( 'B' ) ! B-grid ice dynamics 469 ! stress from ocean U- and V-points to ice U,V point 470 zcoef = cai / cao * 0.5 471 !CDIR COLLAPSE 472 DO jj = 2, jpj 473 DO ji = fs_2, jpi ! vector opt. 474 p_taui(ji,jj) = 0.5 * ( utau(ji-1,jj ) + utau(ji-1,jj-1) ) 475 p_tauj(ji,jj) = 0.5 * ( vtau(ji ,jj-1) + vtau(ji-1,jj-1) ) 476 END DO 477 END DO 478 CALL lbc_lnk( p_taui(:,:), 'I', -1. ) ! I-point (i.e. ice U-V point) 479 CALL lbc_lnk( p_tauj(:,:), 'I', -1. ) ! I-point (i.e. ice U-V point) 480 END SELECT 481 482 483 ! Determine cloud optical depths as a function of latitude (Chou et al., 1981). 484 ! and the correction factor for taking into account the effect of clouds 485 !------------------------------------------------------ 486 !CDIR NOVERRCHK 487 !CDIR COLLAPSE 488 DO jj = 1, jpj 489 !CDIR NOVERRCHK 490 DO ji = 1, jpi 491 ztatm (ji,jj) = sf(jp_tair)%fnow(ji,jj) + rt0 ! air temperature in Kelvins 492 493 zrhoa(ji,jj) = zpatm / ( 287.04 * ztatm(ji,jj) ) ! air density (equation of state for dry air) 494 495 ztamr = ztatm(ji,jj) - rtt ! Saturation water vapour 496 zmt1 = SIGN( 17.269, ztamr ) 497 zmt2 = SIGN( 21.875, ztamr ) 498 zmt3 = SIGN( 28.200, -ztamr ) 499 zes = 611.0 * EXP( ABS( ztamr ) * MIN ( zmt1, zmt2 ) & 500 & / ( ztatm(ji,jj) - 35.86 + MAX( 0.e0, zmt3 ) ) ) 501 502 zev = sf(jp_humi)%fnow(ji,jj) * zes ! vapour pressure 503 zevsqr(ji,jj) = SQRT( zev * 0.01 ) ! square-root of vapour pressure 504 zqatm(ji,jj) = 0.622 * zev / ( zpatm - 0.378 * zev ) ! specific humidity 505 506 !---------------------------------------------------- 507 ! Computation of snow precipitation (Ledley, 1985) | 508 !---------------------------------------------------- 509 zmt1 = 253.0 - ztatm(ji,jj) ; zind1 = MAX( 0.e0, SIGN( 1.e0, zmt1 ) ) 510 zmt2 = ( 272.0 - ztatm(ji,jj) ) / 38.0 ; zind2 = MAX( 0.e0, SIGN( 1.e0, zmt2 ) ) 511 zmt3 = ( 281.0 - ztatm(ji,jj) ) / 18.0 ; zind3 = MAX( 0.e0, SIGN( 1.e0, zmt3 ) ) 512 p_spr(ji,jj) = sf(jp_prec)%fnow(ji,jj) / rday & ! rday = converte mm/day to kg/m2/s 513 & * ( zind1 & ! solid (snow) precipitation [kg/m2/s] 514 & + ( 1.0 - zind1 ) * ( zind2 * ( 0.5 + zmt2 ) & 515 & + ( 1.0 - zind2 ) * zind3 * zmt3 ) ) 516 517 !----------------------------------------------------! 518 ! fraction of net penetrative shortwave radiation ! 519 !----------------------------------------------------! 520 ! fraction of qsr_ice which is NOT absorbed in the thin surface layer 521 ! and thus which penetrates inside the ice cover ( Maykut and Untersteiner, 1971 ; Elbert anbd Curry, 1993 ) 522 p_fr1(ji,jj) = 0.18 * ( 1.e0 - sf(jp_ccov)%fnow(ji,jj) ) + 0.35 * sf(jp_ccov)%fnow(ji,jj) 523 p_fr2(ji,jj) = 0.82 * ( 1.e0 - sf(jp_ccov)%fnow(ji,jj) ) + 0.65 * sf(jp_ccov)%fnow(ji,jj) 524 END DO 525 END DO 526 527 !-----------------------------------------------------------! 528 ! snow/ice Shortwave radiation (abedo already computed) ! 529 !-----------------------------------------------------------! 530 CALL blk_clio_qsr_ice( palb_cs, palb_os, p_qsr ) 531 532 ! ! ========================== ! 533 DO jl = 1, jpl ! Loop over ice categories ! 534 ! ! ========================== ! 535 !CDIR NOVERRCHK 536 !CDIR COLLAPSE 537 DO jj = 1 , jpj 538 !CDIR NOVERRCHK 539 DO ji = 1, jpi 540 !-------------------------------------------! 541 ! long-wave radiation over ice categories ! ( Berliand 1952 ; all latitudes ) 542 !-------------------------------------------! 543 ztatm3 = ztatm(ji,jj) * ztatm(ji,jj) * ztatm(ji,jj) 544 zcldeff = 1.0 - sbudyko(ji,jj) * sf(jp_ccov)%fnow(ji,jj) * sf(jp_ccov)%fnow(ji,jj) 545 ztaevbk = ztatm3 * ztatm(ji,jj) * zcldeff * ( 0.39 - 0.05 * zevsqr(ji,jj) ) 546 ! 547 z_qlw(ji,jj,jl) = - emic * stefan * ( ztaevbk + 4. * ztatm3 * ( pst(ji,jj,jl) - ztatm(ji,jj) ) ) 548 549 !---------------------------------------- 550 ! Turbulent heat fluxes over snow/ice ( Latent and sensible ) 551 !---------------------------------------- 552 553 ! vapour pressure at saturation of ice (tmask to avoid overflow in the exponential) 554 zesi = 611.0 * EXP( 21.8745587 * tmask(ji,jj,1) * ( pst(ji,jj,jl) - rtt )/ ( pst(ji,jj,jl) - 7.66 ) ) 555 ! humidity close to the ice surface (at saturation) 556 zqsati = ( 0.622 * zesi ) / ( zpatm - 0.378 * zesi ) 557 558 ! computation of intermediate values 559 zticemb = pst(ji,jj,jl) - 7.66 560 zticemb2 = zticemb * zticemb 561 ztice3 = pst(ji,jj,jl) * pst(ji,jj,jl) * pst(ji,jj,jl) 562 zdesidt = zesi * ( 9.5 * LOG( 10.0 ) * ( rtt - 7.66 ) / zticemb2 ) 563 564 ! Transfer cofficients assumed to be constant (Parkinson 1979 ; Maykut 1982) 565 zcshi = 1.75e-03 566 zclei = zcshi 567 568 ! sensible and latent fluxes over ice 569 zrhova = zrhoa(ji,jj) * sf(jp_wndm)%fnow(ji,jj) ! computation of intermediate values 570 zrhovaclei = zrhova * zcshi * 2.834e+06 571 zrhovacshi = zrhova * zclei * 1004.0 572 573 ! sensible heat flux 574 z_qsb(ji,jj,jl) = zrhovacshi * ( pst(ji,jj,jl) - ztatm(ji,jj) ) 575 576 ! latent heat flux 577 p_qla(ji,jj,jl) = MAX( 0.e0, zrhovaclei * ( zqsati - zqatm(ji,jj) ) ) 578 579 ! sensitivity of non solar fluxes (dQ/dT) (long-wave, sensible and latent fluxes) 580 zdqlw = 4.0 * emic * stefan * ztice3 581 zdqsb = zrhovacshi 582 zdqla = zrhovaclei * ( zdesidt * ( zqsati * zqsati / ( zesi * zesi ) ) * ( zpatm / 0.622 ) ) 583 ! 584 p_dqla(ji,jj,jl) = zdqla ! latent flux sensitivity 585 p_dqns(ji,jj,jl) = -( zdqlw + zdqsb + zdqla ) ! total non solar sensitivity 586 END DO 587 END DO 588 END DO 589 590 ! 591 ! ----------------------------------------------------------------------------- ! 592 ! III Total FLUXES ! 593 ! ----------------------------------------------------------------------------- ! 594 ! 595 !CDIR COLLAPSE 596 p_qns(:,:,:) = z_qlw (:,:,:) - z_qsb (:,:,:) - p_qla (:,:,:) ! Downward Non Solar flux 597 !CDIR COLLAPSE 598 p_tpr(:,:) = sf(jp_prec)%fnow(:,:) / rday ! total precipitation [kg/m2/s] 599 ! 600 !!gm : not necessary as all input data are lbc_lnk... 601 CALL lbc_lnk( p_fr1 (:,:) , 'T', 1. ) 602 CALL lbc_lnk( p_fr2 (:,:) , 'T', 1. ) 603 DO jl = 1, jpl 604 CALL lbc_lnk( p_qns (:,:,jl) , 'T', 1. ) 605 CALL lbc_lnk( p_dqns(:,:,jl) , 'T', 1. ) 606 CALL lbc_lnk( p_qla (:,:,jl) , 'T', 1. ) 607 CALL lbc_lnk( p_dqla(:,:,jl) , 'T', 1. ) 608 END DO 609 610 !!gm : mask is not required on forcing 611 DO jl = 1, jpl 612 p_qns (:,:,jl) = p_qns (:,:,jl) * tmask(:,:,1) 613 p_qla (:,:,jl) = p_qla (:,:,jl) * tmask(:,:,1) 614 p_dqns(:,:,jl) = p_dqns(:,:,jl) * tmask(:,:,1) 615 p_dqla(:,:,jl) = p_dqla(:,:,jl) * tmask(:,:,1) 616 END DO 617 END SUBROUTINE blk_ice_clio 618 619 620 SUBROUTINE blk_clio_qsr_oce( pqsr_oce ) 621 !!--------------------------------------------------------------------------- 622 !! *** ROUTINE blk_clio_qsr_oce *** 623 !! 624 !! ** Purpose : Computation of the shortwave radiation at the ocean and the 625 !! snow/ice surfaces. 626 !! 627 !! ** Method : - computed qsr from the cloud cover for both ice and ocean 628 !! - also initialise sbudyko and stauc once for all 629 !!---------------------------------------------------------------------- 630 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: pqsr_oce ! shortwave radiation over the ocean 631 !! 632 INTEGER, PARAMETER :: jp24 = 24 ! sampling of the daylight period (sunrise to sunset) into 24 equal parts 633 !! 634 INTEGER :: ji, jj, jt ! dummy loop indices 635 INTEGER :: indaet ! = -1, 0, 1 for odd, normal and leap years resp. 636 INTEGER :: iday ! integer part of day 637 INTEGER :: indxb, indxc ! index for cloud depth coefficient 638 639 REAL(wp) :: zalat , zclat, zcmue, zcmue2 ! local scalars 640 REAL(wp) :: zmt1, zmt2, zmt3 ! 641 REAL(wp) :: zdecl, zsdecl , zcdecl ! 642 REAL(wp) :: za_oce, ztamr, zinda ! 643 644 REAL(wp) :: zdl, zlha ! local scalars 645 REAL(wp) :: zlmunoon, zcldcor, zdaycor ! 646 REAL(wp) :: zxday, zdist, zcoef, zcoef1 ! 647 REAL(wp) :: zes 648 !! 649 REAL(wp), DIMENSION(jpi,jpj) :: zev ! vapour pressure 650 REAL(wp), DIMENSION(jpi,jpj) :: zdlha, zlsrise, zlsset ! 2D workspace 651 652 REAL(wp), DIMENSION(jpi,jpj) :: zps, zpc ! sine (cosine) of latitude per sine (cosine) of solar declination 653 !!--------------------------------------------------------------------- 654 655 656 IF( lbulk_init ) THEN ! Initilization at first time step only 657 rdtbs2 = nn_fsbc * rdt * 0.5 658 ! cloud optical depths as a function of latitude (Chou et al., 1981). 659 ! and the correction factor for taking into account the effect of clouds 660 DO jj = 1, jpj 661 DO ji = 1 , jpi 662 zalat = ( 90.e0 - ABS( gphit(ji,jj) ) ) / 5.e0 663 zclat = ( 95.e0 - gphit(ji,jj) ) / 10.e0 664 indxb = 1 + INT( zalat ) 665 indxc = 1 + INT( zclat ) 666 zdl = zclat - INT( zclat ) 667 ! correction factor to account for the effect of clouds 668 sbudyko(ji,jj) = budyko(indxb) 669 stauc (ji,jj) = ( 1.e0 - zdl ) * tauco( indxc ) + zdl * tauco( indxc + 1 ) 670 END DO 671 END DO 672 IF ( nleapy == 1 ) THEN ; yearday = 366.e0 673 ELSEIF( nleapy == 0 ) THEN ; yearday = 365.e0 674 ELSEIF( nleapy == 30) THEN ; yearday = 360.e0 675 ENDIF 676 lbulk_init = .FALSE. 216 677 ENDIF 217 678 218 END SUBROUTINE sbc_blk_clio 219 220 221 SUBROUTINE blk_oce_clio( ) 222 END SUBROUTINE blk_oce_clio 223 224 225 SUBROUTINE blk_ice_clio 226 END SUBROUTINE blk_ice_clio 227 228 229 SUBROUTINE blk_qsr_clio( ) 230 END SUBROUTINE blk_qsr_clio 679 680 ! Saturated water vapour and vapour pressure 681 ! ------------------------------------------ 682 !CDIR NOVERRCHK 683 !CDIR COLLAPSE 684 DO jj = 1, jpj 685 !CDIR NOVERRCHK 686 DO ji = 1, jpi 687 ztamr = sf(jp_tair)%fnow(ji,jj) + rt0 - rtt 688 zmt1 = SIGN( 17.269, ztamr ) 689 zmt2 = SIGN( 21.875, ztamr ) 690 zmt3 = SIGN( 28.200, -ztamr ) 691 zes = 611.0 * EXP( ABS( ztamr ) * MIN ( zmt1, zmt2 ) & ! Saturation water vapour 692 & / ( sf(jp_tair)%fnow(ji,jj) + rt0 - 35.86 + MAX( 0.e0, zmt3 ) ) ) 693 zev(ji,jj) = sf(jp_humi)%fnow(ji,jj) * zes * 1.0e-05 ! vapour pressure 694 END DO 695 END DO 696 697 !-----------------------------------! 698 ! Computation of solar irradiance ! 699 !-----------------------------------! 700 !!gm : hard coded leap year ??? 701 indaet = 1 ! = -1, 0, 1 for odd, normal and leap years resp. 702 zxday = nday_year + rdtbs2 / rday ! day of the year at which the fluxes are calculated 703 iday = INT( zxday ) ! (centred at the middle of the ice time step) 704 CALL flx_blk_declin( indaet, iday, zdecl ) ! solar declination of the current day 705 zsdecl = SIN( zdecl * rad ) ! its sine 706 zcdecl = COS( zdecl * rad ) ! its cosine 707 708 709 ! correction factor added for computation of shortwave flux to take into account the variation of 710 ! the distance between the sun and the earth during the year (Oberhuber 1988) 711 zdist = zxday * 2. * rpi / yearday 712 zdaycor = 1.0 + 0.0013 * SIN( zdist ) + 0.0342 * COS( zdist ) 713 714 !CDIR NOVERRCHK 715 DO jj = 1, jpj 716 !CDIR NOVERRCHK 717 DO ji = 1, jpi 718 ! product of sine (cosine) of latitude and sine (cosine) of solar declination 719 zps(ji,jj) = SIN( gphit(ji,jj) * rad ) * zsdecl 720 zpc(ji,jj) = COS( gphit(ji,jj) * rad ) * zcdecl 721 ! computation of the both local time of sunrise and sunset 722 zlsrise(ji,jj) = ACOS( - SIGN( 1.e0, zps(ji,jj) ) & 723 & * MIN( 1.e0, SIGN( 1.e0, zps(ji,jj) ) * ( zps(ji,jj) / zpc(ji,jj) ) ) ) 724 zlsset (ji,jj) = - zlsrise(ji,jj) 725 ! dividing the solar day into jp24 segments of length zdlha 726 zdlha (ji,jj) = ( zlsrise(ji,jj) - zlsset(ji,jj) ) / REAL( jp24 ) 727 END DO 728 END DO 729 730 731 !---------------------------------------------! 732 ! shortwave radiation absorbed by the ocean ! 733 !---------------------------------------------! 734 pqsr_oce(:,:) = 0.e0 ! set ocean qsr to zero 735 736 ! compute and sum ocean qsr over the daylight (i.e. between sunrise and sunset) 737 !CDIR NOVERRCHK 738 DO jt = 1, jp24 739 zcoef = FLOAT( jt ) - 0.5 740 !CDIR NOVERRCHK 741 !CDIR COLLAPSE 742 DO jj = 1, jpj 743 !CDIR NOVERRCHK 744 DO ji = 1, jpi 745 zlha = COS( zlsrise(ji,jj) - zcoef * zdlha(ji,jj) ) ! local hour angle 746 zcmue = MAX( 0.e0 , zps(ji,jj) + zpc(ji,jj) * zlha ) ! cos of local solar altitude 747 zcmue2 = 1368.0 * zcmue * zcmue 748 749 ! ocean albedo depending on the cloud cover (Payne, 1972) 750 za_oce = ( 1.0 - sf(jp_ccov)%fnow(ji,jj) ) * 0.05 / ( 1.1 * zcmue**1.4 + 0.15 ) & ! clear sky 751 & + sf(jp_ccov)%fnow(ji,jj) * 0.06 ! overcast 752 753 ! solar heat flux absorbed by the ocean (Zillman, 1972) 754 pqsr_oce(ji,jj) = pqsr_oce(ji,jj) & 755 & + ( 1.0 - za_oce ) * zdlha(ji,jj) * zcmue2 & 756 & / ( ( zcmue + 2.7 ) * zev(ji,jj) + 1.085 * zcmue + 0.10 ) 757 END DO 758 END DO 759 END DO 760 ! Taking into account the ellipsity of the earth orbit, the clouds AND masked if sea-ice cover > 0% 761 !!gm : bug zinda is always 0 si ice.... 762 zcoef1 = srgamma * zdaycor / ( 2. * rpi ) 763 !CDIR COLLAPSE 764 DO jj = 1, jpj 765 DO ji = 1, jpi 766 zlmunoon = ASIN( zps(ji,jj) + zpc(ji,jj) ) / rad ! local noon solar altitude 767 zcldcor = MIN( 1.e0, ( 1.e0 - 0.62 * sf(jp_ccov)%fnow(ji,jj) & ! cloud correction (Reed 1977) 768 & + 0.0019 * zlmunoon ) ) 769 zinda = MAX( 0.e0, SIGN( 1.e0, -( -1.5 + freeze(ji,jj) ) ) ) ! 0 if more than 0% of ice 770 pqsr_oce(ji,jj) = zcoef1 * zcldcor * zinda * pqsr_oce(ji,jj) * tmask(ji,jj,1) ! and zcoef1: ellipsity 771 END DO 772 END DO 773 774 END SUBROUTINE blk_clio_qsr_oce 775 776 777 SUBROUTINE blk_clio_qsr_ice( pa_ice_cs, pa_ice_os, pqsr_ice ) 778 !!--------------------------------------------------------------------------- 779 !! *** ROUTINE blk_clio_qsr_ice *** 780 !! 781 !! ** Purpose : Computation of the shortwave radiation at the ocean and the 782 !! snow/ice surfaces. 783 !! 784 !! ** Method : - computed qsr from the cloud cover for both ice and ocean 785 !! - also initialise sbudyko and stauc once for all 786 !!---------------------------------------------------------------------- 787 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpl) :: pa_ice_cs ! albedo of ice under clear sky 788 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpl) :: pa_ice_os ! albedo of ice under overcast sky 789 REAL(wp), INTENT( out), DIMENSION(jpi,jpj,jpl) :: pqsr_ice ! shortwave radiation over the ice/snow 790 !! 791 INTEGER, PARAMETER :: jp24 = 24 ! sampling of the daylight period (sunrise to sunset) into 24 equal parts 792 !! 793 INTEGER :: ji, jj, jl, jt ! dummy loop indices 794 INTEGER :: indaet ! = -1, 0, 1 for odd, normal and leap years resp. 795 INTEGER :: iday ! integer part of day 796 797 REAL(wp) :: zcmue, zcmue2 ! local scalars 798 REAL(wp) :: zmt1, zmt2, zmt3 ! 799 REAL(wp) :: zdecl, zsdecl , zcdecl ! 800 REAL(wp) :: ztamr ! 801 802 REAL(wp) :: zlha ! local scalars 803 REAL(wp) :: zdaycor ! 804 REAL(wp) :: zxday, zdist, zcoef, zcoef1 ! 805 REAL(wp) :: zes 806 REAL(wp) :: zqsr_ice_cs, zqsr_ice_os 807 !! 808 REAL(wp), DIMENSION(jpi,jpj) :: zev ! vapour pressure 809 REAL(wp), DIMENSION(jpi,jpj) :: zdlha, zlsrise, zlsset ! 2D workspace 810 811 REAL(wp), DIMENSION(jpi,jpj) :: zps, zpc ! sine (cosine) of latitude per sine (cosine) of solar declination 812 !!--------------------------------------------------------------------- 813 814 ! Saturated water vapour and vapour pressure 815 ! ------------------------------------------ 816 !CDIR NOVERRCHK 817 !CDIR COLLAPSE 818 DO jj = 1, jpj 819 !CDIR NOVERRCHK 820 DO ji = 1, jpi 821 ztamr = sf(jp_tair)%fnow(ji,jj) + rt0 - rtt 822 zmt1 = SIGN( 17.269, ztamr ) 823 zmt2 = SIGN( 21.875, ztamr ) 824 zmt3 = SIGN( 28.200, -ztamr ) 825 zes = 611.0 * EXP( ABS( ztamr ) * MIN ( zmt1, zmt2 ) & ! Saturation water vapour 826 & / ( sf(jp_tair)%fnow(ji,jj) + rt0 - 35.86 + MAX( 0.e0, zmt3 ) ) ) 827 zev(ji,jj) = sf(jp_humi)%fnow(ji,jj) * zes * 1.0e-05 ! vapour pressure 828 END DO 829 END DO 830 831 !-----------------------------------! 832 ! Computation of solar irradiance ! 833 !-----------------------------------! 834 !!gm : hard coded leap year ??? 835 indaet = 1 ! = -1, 0, 1 for odd, normal and leap years resp. 836 zxday = nday_year + rdtbs2 / rday ! day of the year at which the fluxes are calculated 837 iday = INT( zxday ) ! (centred at the middle of the ice time step) 838 CALL flx_blk_declin( indaet, iday, zdecl ) ! solar declination of the current day 839 zsdecl = SIN( zdecl * rad ) ! its sine 840 zcdecl = COS( zdecl * rad ) ! its cosine 841 842 843 ! correction factor added for computation of shortwave flux to take into account the variation of 844 ! the distance between the sun and the earth during the year (Oberhuber 1988) 845 zdist = zxday * 2. * rpi / yearday 846 zdaycor = 1.0 + 0.0013 * SIN( zdist ) + 0.0342 * COS( zdist ) 847 848 !CDIR NOVERRCHK 849 DO jj = 1, jpj 850 !CDIR NOVERRCHK 851 DO ji = 1, jpi 852 ! product of sine (cosine) of latitude and sine (cosine) of solar declination 853 zps(ji,jj) = SIN( gphit(ji,jj) * rad ) * zsdecl 854 zpc(ji,jj) = COS( gphit(ji,jj) * rad ) * zcdecl 855 ! computation of the both local time of sunrise and sunset 856 zlsrise(ji,jj) = ACOS( - SIGN( 1.e0, zps(ji,jj) ) & 857 & * MIN( 1.e0, SIGN( 1.e0, zps(ji,jj) ) * ( zps(ji,jj) / zpc(ji,jj) ) ) ) 858 zlsset (ji,jj) = - zlsrise(ji,jj) 859 ! dividing the solar day into jp24 segments of length zdlha 860 zdlha (ji,jj) = ( zlsrise(ji,jj) - zlsset(ji,jj) ) / REAL( jp24 ) 861 END DO 862 END DO 863 864 865 !---------------------------------------------! 866 ! shortwave radiation absorbed by the ice ! 867 !---------------------------------------------! 868 ! compute and sum ice qsr over the daylight for each ice categories 869 pqsr_ice(:,:,:) = 0.e0 870 zcoef1 = zdaycor / ( 2. * rpi ) 871 872 ! !----------------------------! 873 DO jl = 1, jpl ! loop over ice categories ! 874 ! !----------------------------! 875 !CDIR NOVERRCHK 876 DO jt = 1, jp24 877 zcoef = FLOAT( jt ) - 0.5 878 !CDIR NOVERRCHK 879 !CDIR COLLAPSE 880 DO jj = 1, jpj 881 !CDIR NOVERRCHK 882 DO ji = 1, jpi 883 zlha = COS( zlsrise(ji,jj) - zcoef * zdlha(ji,jj) ) ! local hour angle 884 zcmue = MAX( 0.e0 , zps(ji,jj) + zpc(ji,jj) * zlha ) ! cos of local solar altitude 885 zcmue2 = 1368.0 * zcmue * zcmue 886 887 ! solar heat flux absorbed by the ice/snow system (Shine and Crane 1984 adapted to high albedo) 888 zqsr_ice_cs = ( 1.0 - pa_ice_cs(ji,jj,jl) ) * zdlha(ji,jj) * zcmue2 & ! clear sky 889 & / ( ( 1.0 + zcmue ) * zev(ji,jj) + 1.2 * zcmue + 0.0455 ) 890 zqsr_ice_os = zdlha(ji,jj) * SQRT( zcmue ) & ! overcast sky 891 & * ( 53.5 + 1274.5 * zcmue ) * ( 1.0 - 0.996 * pa_ice_os(ji,jj,jl) ) & 892 & / ( 1.0 + 0.139 * stauc(ji,jj) * ( 1.0 - 0.9435 * pa_ice_os(ji,jj,jl) ) ) 893 894 pqsr_ice(ji,jj,jl) = pqsr_ice(ji,jj,jl) + ( ( 1.0 - sf(jp_ccov)%fnow(ji,jj) ) * zqsr_ice_cs & 895 & + sf(jp_ccov)%fnow(ji,jj) * zqsr_ice_os ) 896 END DO 897 END DO 898 END DO 899 ! 900 ! Correction : Taking into account the ellipsity of the earth orbit 901 pqsr_ice(:,:,jl) = pqsr_ice(:,:,jl) * zcoef1 * tmask(:,:,1) 902 ! 903 ! !--------------------------------! 904 END DO ! end loop over ice categories ! 905 ! !--------------------------------! 906 907 908 !!gm : this should be suppress as input data have been passed through lbc_lnk 909 DO jl = 1, jpl 910 CALL lbc_lnk( pqsr_ice(:,:,jl) , 'T', 1. ) 911 END DO 912 ! 913 END SUBROUTINE blk_clio_qsr_ice 231 914 232 915 … … 236 919 !! 237 920 !! ** Purpose : Computation of the solar declination for the day 238 !! kday ( in decimal degrees ).239 921 !! 240 !! ** Method : 922 !! ** Method : ??? 241 923 !!--------------------------------------------------------------------- 242 924 INTEGER , INTENT(in ) :: ky ! = -1, 0, 1 for odd, normal and leap years resp. 243 925 INTEGER , INTENT(in ) :: kday ! day of the year ( kday = 1 on january 1) 244 926 REAL(wp), INTENT( out) :: pdecl ! solar declination 245 246 REAL(wp) :: zday , & ! corresponding day of type year (cf. ky) 247 & zp1, zp2, zp3, zp4 ! temporary scalars 248 REAL(wp) :: a0 = 0.39507671 , & ! constants used in solar 249 & a1 = 22.85684301 , & ! declinaison computation 250 & a2 = -0.38637317 , & 251 & a3 = 0.15096535 , & 252 & a4 = -0.00961411 , & 253 & b1 = -4.29692073 , & 254 & b2 = 0.05702074 , & 255 & b3 = -0.09028607 , & 256 & b4 = 0.00592797 927 !! 928 REAL(wp) :: a0 = 0.39507671 ! coefficients for solar declinaison computation 929 REAL(wp) :: a1 = 22.85684301 ! " "" " 930 REAL(wp) :: a2 = -0.38637317 ! " "" " 931 REAL(wp) :: a3 = 0.15096535 ! " "" " 932 REAL(wp) :: a4 = -0.00961411 ! " "" " 933 REAL(wp) :: b1 = -4.29692073 ! " "" " 934 REAL(wp) :: b2 = 0.05702074 ! " "" " 935 REAL(wp) :: b3 = -0.09028607 ! " "" " 936 REAL(wp) :: b4 = 0.00592797 937 !! 938 REAL(wp) :: zday ! corresponding day of type year (cf. ky) 939 REAL(wp) :: zp ! temporary scalars 257 940 !!--------------------------------------------------------------------- 258 !259 SELECT CASE ( ky )260 CASE ( 1 )261 zday = REAL( kday, wp ) - 0.5262 CASE ( 3 )263 zday = REAL( kday, wp ) - 1.0264 CASE DEFAULT265 zday = REAL( kday, wp )266 END SELECT267 941 268 zp1 = rpi * ( 2.0 * zday - 367.0 ) / yearday 269 zp2 = 2. * zp1 270 zp3 = 3. * zp1 271 zp4 = 4. * zp1 942 IF ( ky == 1 ) THEN ; zday = REAL( kday ) - 0.5 943 ELSEIF( ky == 3 ) THEN ; zday = REAL( kday ) - 1. 944 ELSE ; zday = REAL( kday ) 945 ENDIF 946 947 zp = rpi * ( 2.0 * zday - 367.0 ) / yearday 272 948 273 949 pdecl = a0 & 274 & + a1 * COS( zp 1 ) + a2 * COS( zp2 ) + a3 * COS( zp3 ) + a4 * COS( zp4) &275 & + b1 * SIN( zp 1 ) + b2 * SIN( zp2 ) + b3 * SIN( zp3 ) + b4 * SIN( zp4)950 & + a1 * COS( zp ) + a2 * COS( 2. * zp ) + a3 * COS( 3. * zp ) + a4 * COS( 4. * zp ) & 951 & + b1 * SIN( zp ) + b2 * SIN( 2. * zp ) + b3 * SIN( 3. * zp ) + b4 * SIN( 4. * zp ) 276 952 ! 277 953 END SUBROUTINE flx_blk_declin -
branches/dev_001_SBC/NEMO/OPA_SRC/SBC/sbcice_lim.F90
r751 r840 125 125 SELECT CASE( kblk ) 126 126 CASE( 3 ) ! CLIO bulk formulation 127 CALL blk_ice_clio() 127 CALL blk_ice_clio( sist , ui_ice , vi_ice , alb_ice_cs, alb_ice_os, & 128 & utaui_ice, vtaui_ice , qns_ice , qsr_ice, & 129 & qla_ice , dqns_ice , dqla_ice , & 130 & tprecip , sprecip , & 131 & fr1_i0 , fr2_i0 , 'B' ) 128 132 CASE( 4 ) ! CORE bulk formulation 129 133 CALL blk_ice_core( sist , ui_ice , vi_ice , alb_ice_cs, &
Note: See TracChangeset
for help on using the changeset viewer.