Changeset 7168
- Timestamp:
- 2016-11-02T15:24:08+01:00 (7 years ago)
- Location:
- branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 3 edited
-
SBC/sbccpl.F90 (modified) (19 diffs)
-
SBC/sbcwave.F90 (modified) (8 diffs)
-
step.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
r7166 r7168 23 23 USE sbcapr 24 24 USE sbcdcy ! surface boundary condition: diurnal cycle 25 USE sbcwave ! surface boundary condition: waves 25 26 USE phycst ! physical constants 26 27 #if defined key_lim3 … … 105 106 INTEGER, PARAMETER :: jpr_e3t1st = 41 ! first T level thickness 106 107 INTEGER, PARAMETER :: jpr_fraqsr = 42 ! fraction of solar net radiation absorbed in the first ocean level 107 INTEGER, PARAMETER :: jprcv = 42 ! total number of fields received 108 INTEGER, PARAMETER :: jpr_mslp = 43 ! mean sea level pressure 109 INTEGER, PARAMETER :: jpr_hsig = 44 ! Hsig 110 INTEGER, PARAMETER :: jpr_phioc = 45 ! Wave=>ocean energy flux 111 INTEGER, PARAMETER :: jpr_sdrftx = 46 ! Stokes drift on grid 1 112 INTEGER, PARAMETER :: jpr_sdrfty = 47 ! Stokes drift on grid 2 113 INTEGER, PARAMETER :: jpr_wper = 48 ! Mean wave period 114 INTEGER, PARAMETER :: jpr_wnum = 49 ! Mean wavenumber 115 INTEGER, PARAMETER :: jpr_wstrf = 50 ! Stress fraction adsorbed by waves 116 INTEGER, PARAMETER :: jpr_wdrag = 51 ! Neutral surface drag coefficient 117 INTEGER, PARAMETER :: jprcv = 51 ! total number of fields received 108 118 109 119 INTEGER, PARAMETER :: jps_fice = 1 ! ice fraction sent to the atmosphere … … 135 145 INTEGER, PARAMETER :: jps_e3t1st = 27 ! first level depth (vvl) 136 146 INTEGER, PARAMETER :: jps_fraqsr = 28 ! fraction of solar net radiation absorbed in the first ocean level 137 INTEGER, PARAMETER :: jpsnd = 28 ! total number of fields sended 147 INTEGER, PARAMETER :: jps_ficet = 29 ! total ice fraction 148 INTEGER, PARAMETER :: jps_ocxw = 30 ! currents on grid 1 149 INTEGER, PARAMETER :: jps_ocyw = 31 ! currents on grid 2 150 INTEGER, PARAMETER :: jps_wlev = 32 ! water level 151 INTEGER, PARAMETER :: jpsnd = 32 ! total number of fields sent 138 152 139 153 ! !!** namelist namsbc_cpl ** … … 149 163 ! Received from the atmosphere ! 150 164 TYPE(FLD_C) :: sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf 151 TYPE(FLD_C) :: sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2 165 TYPE(FLD_C) :: sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp 166 ! Send to waves 167 TYPE(FLD_C) :: sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev 168 ! Received from waves 169 TYPE(FLD_C) :: sn_rcv_hsig,sn_rcv_phioc,sn_rcv_sdrfx,sn_rcv_sdrfy,sn_rcv_wper,sn_rcv_wnum,sn_rcv_wstrf,sn_rcv_wdrag 152 170 ! Other namelist parameters ! 153 171 INTEGER :: nn_cplmodel ! Maximum number of models to/from which NEMO is potentialy sending/receiving data … … 161 179 162 180 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: albedo_oce_mix ! ocean albedo sent to atmosphere (mix clear/overcast sky) 181 182 REAL(wp) :: rpref = 101000._wp ! reference atmospheric pressure[N/m2] 183 REAL(wp) :: r1_grau ! = 1.e0 / (grav * rau0) 163 184 164 185 INTEGER , ALLOCATABLE, SAVE, DIMENSION( :) :: nrcvinfo ! OASIS info argument … … 179 200 !! *** FUNCTION sbc_cpl_alloc *** 180 201 !!---------------------------------------------------------------------- 181 INTEGER :: ierr( 3)202 INTEGER :: ierr(4) 182 203 !!---------------------------------------------------------------------- 183 204 ierr(:) = 0 … … 190 211 ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) ) 191 212 ! 213 IF( .NOT. ln_apr_dyn ) ALLOCATE( ssh_ib(jpi,jpj), ssh_ibb(jpi,jpj), apr(jpi, jpj), STAT=ierr(4) ) 214 192 215 sbc_cpl_alloc = MAXVAL( ierr ) 193 216 IF( lk_mpp ) CALL mpp_sum ( sbc_cpl_alloc ) … … 216 239 REAL(wp), POINTER, DIMENSION(:,:) :: zacs, zaos 217 240 !! 218 NAMELIST/namsbc_cpl/ sn_snd_temp, sn_snd_alb , sn_snd_thick, sn_snd_crt , sn_snd_co2, & 219 & sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau , sn_rcv_dqnsdt, sn_rcv_qsr, & 220 & sn_rcv_qns , sn_rcv_emp , sn_rcv_rnf , sn_rcv_cal , sn_rcv_iceflx, & 221 & sn_rcv_co2 , nn_cplmodel , ln_usecplmask 241 NAMELIST/namsbc_cpl/ sn_snd_temp , sn_snd_alb , sn_snd_thick , sn_snd_crt , sn_snd_co2, & 242 & sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau , sn_rcv_dqnsdt, sn_rcv_qsr, & 243 & sn_snd_ifrac, sn_snd_crtw , sn_snd_wlev , sn_rcv_hsig , sn_rcv_phioc , & 244 & sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper , sn_rcv_wnum , sn_rcv_wstrf , & 245 & sn_rcv_wdrag, sn_rcv_qns , sn_rcv_emp , sn_rcv_rnf , sn_rcv_cal , & 246 & sn_rcv_iceflx,sn_rcv_co2 , nn_cplmodel , ln_usecplmask, sn_rcv_mslp 222 247 !!--------------------------------------------------------------------- 223 248 ! … … 260 285 WRITE(numout,*)' sea ice heat fluxes = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')' 261 286 WRITE(numout,*)' atm co2 = ', TRIM(sn_rcv_co2%cldes ), ' (', TRIM(sn_rcv_co2%clcat ), ')' 287 WRITE(numout,*)' significant wave heigth = ', TRIM(sn_rcv_hsig%cldes ), ' (', TRIM(sn_rcv_hsig%clcat ), ')' 288 WRITE(numout,*)' wave to oce energy flux = ', TRIM(sn_rcv_phioc%cldes ), ' (', TRIM(sn_rcv_phioc%clcat ), ')' 289 WRITE(numout,*)' Surface Stokes drift grid u = ', TRIM(sn_rcv_sdrfx%cldes ), ' (', TRIM(sn_rcv_sdrfx%clcat ), ')' 290 WRITE(numout,*)' Surface Stokes drift grid v = ', TRIM(sn_rcv_sdrfy%cldes ), ' (', TRIM(sn_rcv_sdrfy%clcat ), ')' 291 WRITE(numout,*)' Mean wave period = ', TRIM(sn_rcv_wper%cldes ), ' (', TRIM(sn_rcv_wper%clcat ), ')' 292 WRITE(numout,*)' Mean wave number = ', TRIM(sn_rcv_wnum%cldes ), ' (', TRIM(sn_rcv_wnum%clcat ), ')' 293 WRITE(numout,*)' Stress frac adsorbed by waves = ', TRIM(sn_rcv_wstrf%cldes ), ' (', TRIM(sn_rcv_wstrf%clcat ), ')' 294 WRITE(numout,*)' Neutral surf drag coefficient = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')' 262 295 WRITE(numout,*)' sent fields (multiple ice categories)' 263 296 WRITE(numout,*)' surface temperature = ', TRIM(sn_snd_temp%cldes ), ' (', TRIM(sn_snd_temp%clcat ), ')' 264 297 WRITE(numout,*)' albedo = ', TRIM(sn_snd_alb%cldes ), ' (', TRIM(sn_snd_alb%clcat ), ')' 265 298 WRITE(numout,*)' ice/snow thickness = ', TRIM(sn_snd_thick%cldes ), ' (', TRIM(sn_snd_thick%clcat ), ')' 299 WRITE(numout,*)' total ice fraction = ', TRIM(sn_snd_ifrac%cldes ), ' (', TRIM(sn_snd_ifrac%clcat ), ')' 266 300 WRITE(numout,*)' surface current = ', TRIM(sn_snd_crt%cldes ), ' (', TRIM(sn_snd_crt%clcat ), ')' 267 301 WRITE(numout,*)' - referential = ', sn_snd_crt%clvref … … 269 303 WRITE(numout,*)' - mesh = ', sn_snd_crt%clvgrd 270 304 WRITE(numout,*)' oce co2 flux = ', TRIM(sn_snd_co2%cldes ), ' (', TRIM(sn_snd_co2%clcat ), ')' 305 WRITE(numout,*)' water level = ', TRIM(sn_snd_wlev%cldes ), ' (', TRIM(sn_snd_wlev%clcat ), ')' 306 WRITE(numout,*)' mean sea level pressure = ', TRIM(sn_rcv_mslp%cldes ), ' (', TRIM(sn_rcv_mslp%clcat ), ')' 307 WRITE(numout,*)' surface current to waves = ', TRIM(sn_snd_crtw%cldes ), ' (', TRIM(sn_snd_crtw%clcat ), ')' 308 WRITE(numout,*)' - referential = ', sn_snd_crtw%clvref 309 WRITE(numout,*)' - orientation = ', sn_snd_crtw%clvor 310 WRITE(numout,*)' - mesh = ', sn_snd_crtw%clvgrd 271 311 WRITE(numout,*)' nn_cplmodel = ', nn_cplmodel 272 312 WRITE(numout,*)' ln_usecplmask = ', ln_usecplmask … … 470 510 ! ! ------------------------- ! 471 511 srcv(jpr_co2 )%clname = 'O_AtmCO2' ; IF( TRIM(sn_rcv_co2%cldes ) == 'coupled' ) srcv(jpr_co2 )%laction = .TRUE. 512 513 ! ! ------------------------- ! 514 ! ! Mean Sea Level Pressure ! 515 ! ! ------------------------- ! 516 srcv(jpr_mslp)%clname = 'O_MSLP' ; IF( TRIM(sn_rcv_mslp%cldes ) == 'coupled' ) srcv(jpr_mslp)%laction = .TRUE. 517 472 518 ! ! ------------------------- ! 473 519 ! ! topmelt and botmelt ! … … 483 529 srcv(jpr_topm:jpr_botm)%laction = .TRUE. 484 530 ENDIF 531 ! ! ------------------------- ! 532 ! ! Wave breaking ! 533 ! ! ------------------------- ! 534 srcv(jpr_hsig)%clname = 'O_Hsigwa' ! significant wave height 535 IF( TRIM(sn_rcv_hsig%cldes ) == 'coupled' ) THEN 536 srcv(jpr_hsig)%laction = .TRUE. 537 cpl_hsig = .TRUE. 538 ENDIF 539 srcv(jpr_phioc)%clname = 'O_PhiOce' ! wave to ocean energy 540 IF( TRIM(sn_rcv_phioc%cldes ) == 'coupled' ) THEN 541 srcv(jpr_phioc)%laction = .TRUE. 542 cpl_phioc = .TRUE. 543 ENDIF 544 srcv(jpr_sdrftx)%clname = 'O_Sdrfx' ! Stokes drift in the u direction 545 IF( TRIM(sn_rcv_sdrfx%cldes ) == 'coupled' ) THEN 546 srcv(jpr_sdrftx)%laction = .TRUE. 547 cpl_sdrftx = .TRUE. 548 ENDIF 549 srcv(jpr_sdrfty)%clname = 'O_Sdrfy' ! Stokes drift in the v direction 550 IF( TRIM(sn_rcv_sdrfy%cldes ) == 'coupled' ) THEN 551 srcv(jpr_sdrfty)%laction = .TRUE. 552 cpl_sdrfty = .TRUE. 553 ENDIF 554 srcv(jpr_wper)%clname = 'O_WPer' ! mean wave period 555 IF( TRIM(sn_rcv_wper%cldes ) == 'coupled' ) THEN 556 srcv(jpr_wper)%laction = .TRUE. 557 cpl_wper = .TRUE. 558 ENDIF 559 srcv(jpr_wnum)%clname = 'O_WNum' ! mean wave number 560 IF( TRIM(sn_rcv_wnum%cldes ) == 'coupled' ) THEN 561 srcv(jpr_wnum)%laction = .TRUE. 562 cpl_wnum = .TRUE. 563 ENDIF 564 srcv(jpr_wstrf)%clname = 'O_WStrf' ! stress fraction adsorbed by the wave 565 IF( TRIM(sn_rcv_wstrf%cldes ) == 'coupled' ) THEN 566 srcv(jpr_wstrf)%laction = .TRUE. 567 cpl_wstrf = .TRUE. 568 ENDIF 569 srcv(jpr_wdrag)%clname = 'O_WDrag' ! neutral surface drag coefficient 570 IF( TRIM(sn_rcv_wdrag%cldes ) == 'coupled' ) THEN 571 srcv(jpr_wdrag)%laction = .TRUE. 572 cpl_wdrag = .TRUE. 573 ENDIF 574 ! 485 575 ! ! ------------------------------- ! 486 576 ! ! OPA-SAS coupling - rcv by opa ! … … 637 727 ! ! ------------------------- ! 638 728 ssnd(jps_fice)%clname = 'OIceFrc' 729 ssnd(jps_ficet)%clname = 'OIceFrcT' 639 730 ssnd(jps_hice)%clname = 'OIceTck' 640 731 ssnd(jps_hsnw)%clname = 'OSnwTck' … … 645 736 ENDIF 646 737 738 IF (TRIM( sn_snd_ifrac%cldes ) == 'coupled') ssnd(jps_ficet)%laction = .TRUE. 739 647 740 SELECT CASE ( TRIM( sn_snd_thick%cldes ) ) 648 741 CASE( 'none' ) ! nothing to do … … 665 758 ssnd(jps_ocy1)%clname = 'O_OCury1' ; ssnd(jps_ivy1)%clname = 'O_IVely1' 666 759 ssnd(jps_ocz1)%clname = 'O_OCurz1' ; ssnd(jps_ivz1)%clname = 'O_IVelz1' 760 ssnd(jps_ocxw)%clname = 'O_OCurxw' 761 ssnd(jps_ocyw)%clname = 'O_OCuryw' 667 762 ! 668 763 ssnd(jps_ocx1:jps_ivz1)%nsgn = -1. ! vectors: change of the sign at the north fold … … 685 780 END SELECT 686 781 782 ssnd(jps_ocxw:jps_ocyw)%nsgn = -1. ! vectors: change of the sign at the north fold 783 784 IF( sn_snd_crtw%clvgrd == 'U,V' ) THEN 785 ssnd(jps_ocxw)%clgrid = 'U' ; ssnd(jps_ocyw)%clgrid = 'V' 786 ELSE IF( sn_snd_crtw%clvgrd /= 'T' ) THEN 787 CALL ctl_stop( 'sn_snd_crtw%clvgrd must be equal to T' ) 788 ENDIF 789 IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) ssnd(jps_ocxw:jps_ocyw)%nsgn = 1. 790 SELECT CASE( TRIM( sn_snd_crtw%cldes ) ) 791 CASE( 'none' ) ; ssnd(jps_ocxw:jps_ocyw)%laction = .FALSE. 792 CASE( 'oce only' ) ; ssnd(jps_ocxw:jps_ocyw)%laction = .TRUE. 793 CASE( 'weighted oce and ice' ) ! nothing to do 794 CASE( 'mixed oce-ice' ) ; ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE. 795 CASE default ; CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_crtw%cldes' ) 796 END SELECT 797 687 798 ! ! ------------------------- ! 688 799 ! ! CO2 flux ! 689 800 ! ! ------------------------- ! 690 801 ssnd(jps_co2)%clname = 'O_CO2FLX' ; IF( TRIM(sn_snd_co2%cldes) == 'coupled' ) ssnd(jps_co2 )%laction = .TRUE. 802 803 ! ! ------------------------- ! 804 ! ! Sea surface height ! 805 ! ! ------------------------- ! 806 ssnd(jps_wlev)%clname = 'O_Wlevel' ; IF( TRIM(sn_snd_wlev%cldes) == 'coupled' ) ssnd(jps_wlev)%laction = .TRUE. 691 807 692 808 ! ! ------------------------------- ! … … 837 953 !! emp upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case) 838 954 !!---------------------------------------------------------------------- 955 USE zdf_oce, ONLY : ln_zdfqiao 956 957 IMPLICIT NONE 958 839 959 INTEGER, INTENT(in) :: kt ! ocean model time step index 840 960 INTEGER, INTENT(in) :: k_fsbc ! frequency of sbc (-> ice model) computation … … 992 1112 IF( srcv(jpr_co2)%laction ) atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1) 993 1113 #endif 1114 ! 1115 ! ! ========================= ! 1116 ! ! Mean Sea Level Pressure ! (taum) 1117 ! ! ========================= ! 1118 ! 1119 IF( srcv(jpr_mslp)%laction ) THEN ! UKMO SHELF effect of atmospheric pressure on SSH 1120 IF( kt /= nit000 ) ssh_ibb(:,:) = ssh_ib(:,:) !* Swap of ssh_ib fields 1121 1122 r1_grau = 1.e0 / (grav * rau0) !* constant for optimization 1123 ssh_ib(:,:) = - ( frcv(jpr_mslp)%z3(:,:,1) - rpref ) * r1_grau ! equivalent ssh (inverse barometer) 1124 apr (:,:) = frcv(jpr_mslp)%z3(:,:,1) !atmospheric pressure 1125 1126 IF( kt == nit000 ) ssh_ibb(:,:) = ssh_ib(:,:) ! correct this later (read from restart if possible) 1127 END IF 1128 ! 1129 IF( ln_sdw ) THEN ! Stokes Drift correction activated 1130 ! ! ========================= ! 1131 ! ! Stokes drift u ! 1132 ! ! ========================= ! 1133 IF( srcv(jpr_sdrftx)%laction ) zusd2dt(:,:) = frcv(jpr_sdrftx)%z3(:,:,1) 1134 ! 1135 ! ! ========================= ! 1136 ! ! Stokes drift v ! 1137 ! ! ========================= ! 1138 IF( srcv(jpr_sdrfty)%laction ) zvsd2dt(:,:) = frcv(jpr_sdrfty)%z3(:,:,1) 1139 ! 1140 ! ! ========================= ! 1141 ! ! Wave mean period ! 1142 ! ! ========================= ! 1143 IF( srcv(jpr_wper)%laction ) wmp(:,:) = frcv(jpr_wper)%z3(:,:,1) 1144 ! 1145 ! ! ========================= ! 1146 ! ! Significant wave height ! 1147 ! ! ========================= ! 1148 IF( srcv(jpr_hsig)%laction ) swh(:,:) = frcv(jpr_hsig)%z3(:,:,1) 1149 ! 1150 ! ! ========================= ! 1151 ! ! Vertical mixing Qiao ! 1152 ! ! ========================= ! 1153 IF( srcv(jpr_wnum)%laction .AND. ln_zdfqiao ) wnum(:,:) = frcv(jpr_wnum)%z3(:,:,1) 1154 1155 ! Calculate the 3D Stokes drift both in coupled and not fully uncoupled mode 1156 IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. srcv(jpr_wper)%laction & 1157 .OR. srcv(jpr_hsig)%laction ) THEN 1158 CALL sbc_stokes() 1159 IF( ln_zdfqiao .AND. .NOT. srcv(jpr_wnum)%laction ) CALL sbc_qiao() 1160 ENDIF 1161 IF( ln_zdfqiao .AND. srcv(jpr_wnum)%laction ) CALL sbc_qiao() 1162 ENDIF 1163 ! ! ========================= ! 1164 ! ! Stress adsorbed by waves ! 1165 ! ! ========================= ! 1166 IF( srcv(jpr_wstrf)%laction .AND. ln_tauoc ) tauoc_wave(:,:) = frcv(jpr_wstrf)%z3(:,:,1) 1167 1168 ! ! ========================= ! 1169 ! ! Wave drag coefficient ! 1170 ! ! ========================= ! 1171 IF( srcv(jpr_wdrag)%laction .AND. ln_cdgw ) cdn_wave(:,:) = frcv(jpr_wdrag)%z3(:,:,1) 994 1172 995 1173 ! Fields received by SAS when OASIS coupling … … 1984 2162 ENDIF 1985 2163 ! 2164 ! ! ------------------------- ! 2165 ! ! Surface current to waves ! 2166 ! ! ------------------------- ! 2167 IF( ssnd(jps_ocxw)%laction .OR. ssnd(jps_ocyw)%laction ) THEN 2168 ! 2169 ! j+1 j -----V---F 2170 ! surface velocity always sent from T point ! | 2171 ! j | T U 2172 ! | | 2173 ! j j-1 -I-------| 2174 ! (for I) | | 2175 ! i-1 i i 2176 ! i i+1 (for I) 2177 SELECT CASE( TRIM( sn_snd_crtw%cldes ) ) 2178 CASE( 'oce only' ) ! C-grid ==> T 2179 DO jj = 2, jpjm1 2180 DO ji = fs_2, fs_jpim1 ! vector opt. 2181 zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj ,1) ) 2182 zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji , jj-1,1) ) 2183 END DO 2184 END DO 2185 CASE( 'weighted oce and ice' ) 2186 SELECT CASE ( cp_ice_msh ) 2187 CASE( 'C' ) ! Ocean and Ice on C-grid ==> T 2188 DO jj = 2, jpjm1 2189 DO ji = fs_2, fs_jpim1 ! vector opt. 2190 zotx1(ji,jj) = 0.5 * ( un (ji,jj,1) + un (ji-1,jj ,1) ) * zfr_l(ji,jj) 2191 zoty1(ji,jj) = 0.5 * ( vn (ji,jj,1) + vn (ji ,jj-1,1) ) * zfr_l(ji,jj) 2192 zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * fr_i(ji,jj) 2193 zity1(ji,jj) = 0.5 * ( v_ice(ji,jj ) + v_ice(ji ,jj-1 ) ) * fr_i(ji,jj) 2194 END DO 2195 END DO 2196 CASE( 'I' ) ! Ocean on C grid, Ice on I-point (B-grid) ==> T 2197 DO jj = 2, jpjm1 2198 DO ji = 2, jpim1 ! NO vector opt. 2199 zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj ,1) ) * zfr_l(ji,jj) 2200 zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji ,jj-1,1) ) * zfr_l(ji,jj) 2201 zitx1(ji,jj) = 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1) & 2202 & + u_ice(ji+1,jj ) + u_ice(ji,jj ) ) * fr_i(ji,jj) 2203 zity1(ji,jj) = 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1) & 2204 & + v_ice(ji+1,jj ) + v_ice(ji,jj ) ) * fr_i(ji,jj) 2205 END DO 2206 END DO 2207 CASE( 'F' ) ! Ocean on C grid, Ice on F-point (B-grid) ==> T 2208 DO jj = 2, jpjm1 2209 DO ji = 2, jpim1 ! NO vector opt. 2210 zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj ,1) ) * zfr_l(ji,jj) 2211 zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji ,jj-1,1) ) * zfr_l(ji,jj) 2212 zitx1(ji,jj) = 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1) & 2213 & + u_ice(ji-1,jj ) + u_ice(ji,jj ) ) * fr_i(ji,jj) 2214 zity1(ji,jj) = 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1) & 2215 & + v_ice(ji-1,jj ) + v_ice(ji,jj ) ) * fr_i(ji,jj) 2216 END DO 2217 END DO 2218 END SELECT 2219 CALL lbc_lnk( zitx1, 'T', -1. ) ; CALL lbc_lnk( zity1, 'T', -1. ) 2220 CASE( 'mixed oce-ice' ) 2221 SELECT CASE ( cp_ice_msh ) 2222 CASE( 'C' ) ! Ocean and Ice on C-grid ==> T 2223 DO jj = 2, jpjm1 2224 DO ji = fs_2, fs_jpim1 ! vector opt. 2225 zotx1(ji,jj) = 0.5 * ( un (ji,jj,1) + un (ji-1,jj ,1) ) * zfr_l(ji,jj) & 2226 & + 0.5 * ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * fr_i(ji,jj) 2227 zoty1(ji,jj) = 0.5 * ( vn (ji,jj,1) + vn (ji ,jj-1,1) ) * zfr_l(ji,jj) & 2228 & + 0.5 * ( v_ice(ji,jj ) + v_ice(ji ,jj-1 ) ) * fr_i(ji,jj) 2229 END DO 2230 END DO 2231 CASE( 'I' ) ! Ocean on C grid, Ice on I-point (B-grid) ==> T 2232 DO jj = 2, jpjm1 2233 DO ji = 2, jpim1 ! NO vector opt. 2234 zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj ,1) ) * zfr_l(ji,jj) & 2235 & + 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1) & 2236 & + u_ice(ji+1,jj ) + u_ice(ji,jj ) ) * fr_i(ji,jj) 2237 zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji ,jj-1,1) ) * zfr_l(ji,jj) & 2238 & + 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1) & 2239 & + v_ice(ji+1,jj ) + v_ice(ji,jj ) ) * fr_i(ji,jj) 2240 END DO 2241 END DO 2242 CASE( 'F' ) ! Ocean on C grid, Ice on F-point (B-grid) ==> T 2243 DO jj = 2, jpjm1 2244 DO ji = 2, jpim1 ! NO vector opt. 2245 zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj ,1) ) * zfr_l(ji,jj) & 2246 & + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1) & 2247 & + u_ice(ji-1,jj ) + u_ice(ji,jj ) ) * fr_i(ji,jj) 2248 zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji ,jj-1,1) ) * zfr_l(ji,jj) & 2249 & + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1) & 2250 & + v_ice(ji-1,jj ) + v_ice(ji,jj ) ) * fr_i(ji,jj) 2251 END DO 2252 END DO 2253 END SELECT 2254 END SELECT 2255 CALL lbc_lnk( zotx1, ssnd(jps_ocxw)%clgrid, -1. ) ; CALL lbc_lnk( zoty1, ssnd(jps_ocyw)%clgrid, -1. ) 2256 ! 2257 ! 2258 IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) THEN ! Rotation of the components 2259 ! ! Ocean component 2260 CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->e', ztmp1 ) ! 1st component 2261 CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->n', ztmp2 ) ! 2nd component 2262 zotx1(:,:) = ztmp1(:,:) ! overwrite the components 2263 zoty1(:,:) = ztmp2(:,:) 2264 IF( ssnd(jps_ivx1)%laction ) THEN ! Ice component 2265 CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 ) ! 1st component 2266 CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 ) ! 2nd component 2267 zitx1(:,:) = ztmp1(:,:) ! overwrite the components 2268 zity1(:,:) = ztmp2(:,:) 2269 ENDIF 2270 ENDIF 2271 ! 2272 ! ! spherical coordinates to cartesian -> 2 components to 3 components 2273 ! IF( TRIM( sn_snd_crtw%clvref ) == 'cartesian' ) THEN 2274 ! ztmp1(:,:) = zotx1(:,:) ! ocean currents 2275 ! ztmp2(:,:) = zoty1(:,:) 2276 ! CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 ) 2277 ! ! 2278 ! IF( ssnd(jps_ivx1)%laction ) THEN ! ice velocities 2279 ! ztmp1(:,:) = zitx1(:,:) 2280 ! ztmp1(:,:) = zity1(:,:) 2281 ! CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 ) 2282 ! ENDIF 2283 ! ENDIF 2284 ! 2285 IF( ssnd(jps_ocxw)%laction ) CALL cpl_snd( jps_ocxw, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info ) ! ocean x current 1st grid 2286 IF( ssnd(jps_ocyw)%laction ) CALL cpl_snd( jps_ocyw, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info ) ! ocean y current 1st grid 2287 ! 2288 ENDIF 2289 ! 2290 IF( ssnd(jps_ficet)%laction ) THEN 2291 CALL cpl_snd( jps_ficet, isec, RESHAPE ( fr_i, (/jpi,jpj,1/) ), info ) 2292 END IF 2293 ! ! ------------------------- ! 2294 ! ! Water levels to waves ! 2295 ! ! ------------------------- ! 2296 IF( ssnd(jps_wlev)%laction ) THEN 2297 IF( ln_apr_dyn ) THEN 2298 IF( kt /= nit000 ) THEN 2299 ztmp1(:,:) = sshb(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) 2300 ELSE 2301 ztmp1(:,:) = sshb(:,:) 2302 ENDIF 2303 ELSE 2304 ztmp1(:,:) = sshn(:,:) 2305 ENDIF 2306 CALL cpl_snd( jps_wlev , isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info ) 2307 END IF 1986 2308 ! 1987 2309 ! Fields sent by OPA to SAS when doing OPA<->SAS coupling -
branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
r7167 r7168 27 27 PRIVATE 28 28 29 PUBLIC sbc_stokes, sbc_qiao ! routines called in sbccpl 29 30 PUBLIC sbc_wave ! routine called in sbcmod 30 31 31 INTEGER , PARAMETER :: jpfld = 4 ! number of files to read for stokes drift 32 INTEGER , PARAMETER :: jp_usd = 1 ! index of stokes drift (i-component) (m/s) at T-point 33 INTEGER , PARAMETER :: jp_vsd = 2 ! index of stokes drift (j-component) (m/s) at T-point 34 INTEGER , PARAMETER :: jp_swh = 3 ! index of significant wave hight (m) at T-point 35 INTEGER , PARAMETER :: jp_wmp = 4 ! index of mean wave period (s) at T-point 32 ! Variables checking if the wave parameters are coupled (if not, they are read from file) 33 LOGICAL, PUBLIC :: cpl_hsig=.FALSE. 34 LOGICAL, PUBLIC :: cpl_phioc=.FALSE. 35 LOGICAL, PUBLIC :: cpl_sdrftx=.FALSE. 36 LOGICAL, PUBLIC :: cpl_sdrfty=.FALSE. 37 LOGICAL, PUBLIC :: cpl_wper=.FALSE. 38 LOGICAL, PUBLIC :: cpl_wnum=.FALSE. 39 LOGICAL, PUBLIC :: cpl_wstrf=.FALSE. 40 LOGICAL, PUBLIC :: cpl_wdrag=.FALSE. 41 42 INTEGER :: jpfld ! number of files to read for stokes drift 43 INTEGER :: jp_usd ! index of stokes drift (i-component) (m/s) at T-point 44 INTEGER :: jp_vsd ! index of stokes drift (j-component) (m/s) at T-point 45 INTEGER :: jp_swh ! index of significant wave hight (m) at T-point 46 INTEGER :: jp_wmp ! index of mean wave period (s) at T-point 36 47 37 48 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_cd ! structure of input fields (file informations, fields read) Drag Coefficient … … 57 68 CONTAINS 58 69 70 SUBROUTINE sbc_stokes( ) 71 !!--------------------------------------------------------------------- 72 !! *** ROUTINE sbc_stokes *** 73 !! 74 !! ** Purpose : compute the 3d Stokes Drift according to Breivik et al., 75 !! 2014 (DOI: 10.1175/JPO-D-14-0020.1) 76 !! 77 !! ** Method : - Calculate Stokes transport speed 78 !! - Calculate horizontal divergence 79 !! - Integrate the horizontal divergenze from the bottom 80 !! ** action 81 !!--------------------------------------------------------------------- 82 INTEGER :: jj,ji,jk 83 REAL(wp) :: ztransp, zsp0, zk, zus, zvs 84 REAL(wp), DIMENSION(jpi,jpj) :: zfac 85 REAL(wp), DIMENSION(:,:,:), POINTER :: ze3hdiv ! 3D workspace 86 87 88 CALL wrk_alloc( jpi,jpj,jpk, ze3hdiv ) 89 DO jk = 1, jpk 90 DO jj = 1, jpj 91 DO ji = 1, jpi 92 ! On T grid 93 ! Stokes transport speed estimated from Hs and Tmean 94 ztransp = 2.0_wp*rpi*swh(ji,jj)**2.0_wp/(16.0_wp*MAX(wmp(ji,jj),0.0000001_wp)) 95 ! Stokes surface speed 96 zsp0 = SQRT( sf_sd(jp_usd)%fnow(ji,jj,1)**2 + sf_sd(jp_vsd)%fnow(ji,jj,1)**2) 97 ! Wavenumber scale 98 zk = ABS(zsp0)/MAX(ABS(5.97_wp*ztransp),0.0000001_wp) 99 ! Depth attenuation 100 zfac(ji,jj) = EXP(-2.0_wp*zk*fsdept(ji,jj,jk))/(1.0_wp+8.0_wp*zk*fsdept(ji,jj,jk)) 101 END DO 102 END DO 103 104 DO jj = 1, jpjm1 105 DO ji = 1, jpim1 106 ! Into the U and V Grid 107 zus = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zfac(ji,jj) * tmask(ji,jj,1) & 108 & + zfac(ji+1,jj) * tmask(ji+1,jj,1) ) 109 zvs = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zfac(ji,jj) * tmask(ji,jj,1) & 110 & + zfac(ji,jj+1) * tmask(ji,jj+1,1) ) 111 usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zusd2dt(ji,jj) * tmask(ji,jj,1) & 112 & + zusd2dt(ji+1,jj) * tmask(ji+1,jj,1) ) 113 vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zvsd2dt(ji,jj) * tmask(ji,jj,1) & 114 & + zvsd2dt(ji,jj+1) * tmask(ji,jj+1,1) ) 115 usd3d(ji,jj,jk) = usd2d(ji,jj) * zus 116 vsd3d(ji,jj,jk) = vsd2d(ji,jj) * zvs 117 END DO 118 END DO 119 END DO 120 ! 121 CALL lbc_lnk( usd3d(:,:,:), 'U', -1. ) 122 CALL lbc_lnk( vsd3d(:,:,:), 'V', -1. ) 123 ! 124 DO jk = 1, jpkm1 ! e3t * Horizontal divergence 125 DO jj = 2, jpjm1 126 DO ji = fs_2, fs_jpim1 ! vector opt. 127 ze3hdiv(ji,jj,jk) = ( e2u(ji ,jj) * fse3u_n(ji ,jj,jk) * usd3d(ji ,jj,jk) & 128 & - e2u(ji-1,jj) * fse3u_n(ji-1,jj,jk) * usd3d(ji-1,jj,jk) & 129 & + e1v(ji,jj ) * fse3v_n(ji,jj ,jk) * vsd3d(ji,jj ,jk) & 130 & - e1v(ji,jj-1) * fse3v_n(ji,jj-1,jk) * vsd3d(ji,jj-1,jk) ) * r1_e1e2t(ji,jj) 131 END DO 132 END DO 133 IF( .NOT. AGRIF_Root() ) THEN 134 IF( nbondi == 1 .OR. nbondi == 2 ) ze3hdiv(nlci-1, : ,jk) = 0._wp ! east 135 IF( nbondi == -1 .OR. nbondi == 2 ) ze3hdiv( 2 , : ,jk) = 0._wp ! west 136 IF( nbondj == 1 .OR. nbondj == 2 ) ze3hdiv( : ,nlcj-1,jk) = 0._wp ! north 137 IF( nbondj == -1 .OR. nbondj == 2 ) ze3hdiv( : , 2 ,jk) = 0._wp ! south 138 ENDIF 139 END DO 140 CALL lbc_lnk( ze3hdiv, 'T', 1. ) 141 ! 142 DO jk = jpkm1, 1, -1 !* integrate from the bottom the e3t * hor. divergence 143 wsd3d(:,:,jk) = wsd3d(:,:,jk+1) - ze3hdiv(:,:,jk) 144 END DO 145 #if defined key_bdy 146 IF( lk_bdy ) THEN 147 DO jk = 1, jpkm1 148 wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:) 149 END DO 150 ENDIF 151 #endif 152 CALL wrk_dealloc( jpi,jpj,jpk, ze3hdiv ) 153 ! 154 END SUBROUTINE sbc_stokes 155 156 SUBROUTINE sbc_qiao( ) 157 !!--------------------------------------------------------------------- 158 !! *** ROUTINE sbc_qiao *** 159 !! 160 !! ** Purpose : Qiao formulation for wave enhanced turbulence 161 !! 2010 (DOI: 10.1007/s10236-010-0326) 162 !! 163 !! ** Method : - 164 !! ** action 165 !!--------------------------------------------------------------------- 166 INTEGER :: jj,ji 167 168 ! Calculate the module of the stokes drift on T grid 169 !------------------------------------------------- 170 DO jj = 1, jpj 171 DO ji = 1, jpi 172 tsd2d(ji,jj) = ((zusd2dt(ji,jj) * tmask(ji,jj,1))**2.0 + & 173 & (zvsd2dt(ji,jj) * tmask(ji,jj,1))**2.0)**0.5 174 END DO 175 END DO 176 ! 177 END SUBROUTINE sbc_qiao 178 59 179 SUBROUTINE sbc_wave( kt ) 60 180 !!--------------------------------------------------------------------- … … 71 191 !! ** action 72 192 !!--------------------------------------------------------------------- 73 USE zdf_oce, ONLY : ln_zdfqiao 193 USE zdf_oce , ONLY : ln_zdfqiao 194 195 IMPLICIT NONE 74 196 75 197 INTEGER, INTENT( in ) :: kt ! ocean time step 76 198 ! 77 199 INTEGER :: ierror ! return error code 78 INTEGER :: ifpr , jj,ji,jk200 INTEGER :: ifpr 79 201 INTEGER :: ios ! Local integer output status for namelist read 80 202 ! 81 203 CHARACTER(len=100) :: cn_dir ! Root directory for location of drag coefficient files 82 REAL(wp) :: ztransp, zsp0, zk, zus, zvs 83 REAL(wp), DIMENSION(jpi,jpj) :: zfac 84 REAL(wp), DIMENSION(:,:,:), POINTER :: ze3hdiv ! 3D workspace 85 TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read 204 TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: slf_i ! array of namelist informations on the fields to read 86 205 TYPE(FLD_N) :: sn_cdg, sn_usd, sn_vsd, & 87 206 & sn_swh, sn_wmp, sn_wnum, sn_tauoc ! informations about the fields to be read … … 103 222 ! 104 223 IF ( ln_cdgw ) THEN 105 ALLOCATE( sf_cd(1), STAT=ierror ) !* allocate and fill sf_wave with sn_cdg 106 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 107 ! 108 ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1) ) 109 IF( sn_cdg%ln_tint ) ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 110 CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 224 IF ( .NOT. cpl_wdrag ) THEN 225 ALLOCATE( sf_cd(1), STAT=ierror ) !* allocate and fill sf_wave with sn_cdg 226 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 227 ! 228 ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1) ) 229 IF( sn_cdg%ln_tint ) ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 230 CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 231 ENDIF 111 232 ALLOCATE( cdn_wave(jpi,jpj) ) 112 233 cdn_wave(:,:) = 0.0 … … 114 235 115 236 IF ( ln_tauoc ) THEN 116 ALLOCATE( sf_tauoc(1), STAT=ierror ) !* allocate and fill sf_wave with sn_tauoc 117 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 118 ! 119 ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1) ) 120 IF( sn_tauoc%ln_tint ) ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) ) 121 CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 237 IF ( .NOT. cpl_wstrf ) THEN 238 ALLOCATE( sf_tauoc(1), STAT=ierror ) !* allocate and fill sf_wave with sn_tauoc 239 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 240 ! 241 ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1) ) 242 IF( sn_tauoc%ln_tint ) ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) ) 243 CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 244 ENDIF 122 245 ALLOCATE( tauoc_wave(jpi,jpj) ) 123 246 tauoc_wave(:,:) = 0.0 124 ENDIF247 ENDIF 125 248 126 249 IF ( ln_sdw ) THEN 127 slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd; 128 slf_i(jp_swh) = sn_swh ; slf_i(jp_wmp) = sn_wmp; 129 ALLOCATE( sf_sd(jpfld), STAT=ierror ) !* allocate and fill sf_sd with stokes drift 130 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 131 ! 132 DO ifpr= 1, jpfld 133 ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 134 IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 135 END DO 136 137 CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 250 ! Find out how many fields have to be read from file if not coupled 251 jpfld=0 252 jp_usd=0; jp_vsd=0; jp_swh=0; jp_wmp=0 253 IF( .NOT. cpl_sdrftx ) THEN 254 jpfld=jpfld+1 255 jp_usd=jpfld 256 ENDIF 257 IF( .NOT. cpl_sdrfty ) THEN 258 jpfld=jpfld+1 259 jp_vsd=jpfld 260 ENDIF 261 IF( .NOT. cpl_hsig ) THEN 262 jpfld=jpfld+1 263 jp_swh=jpfld 264 ENDIF 265 IF( .NOT. cpl_wper ) THEN 266 jpfld=jpfld+1 267 jp_wmp=jpfld 268 ENDIF 269 270 ! Read from file only the non-coupled fields 271 IF( jpfld > 0 ) THEN 272 ALLOCATE( slf_i(jpfld) ) 273 IF( jp_usd > 0 ) slf_i(jp_usd) = sn_usd 274 IF( jp_vsd > 0 ) slf_i(jp_vsd) = sn_vsd 275 IF( jp_swh > 0 ) slf_i(jp_swh) = sn_swh 276 IF( jp_wmp > 0 ) slf_i(jp_wmp) = sn_wmp 277 ALLOCATE( sf_sd(jpfld), STAT=ierror ) !* allocate and fill sf_sd with stokes drift 278 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 279 ! 280 DO ifpr= 1, jpfld 281 ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 282 IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 283 END DO 284 285 CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 286 ENDIF 138 287 ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj) ) 139 288 ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) ) … … 145 294 swh (:,:) = 0._wp ; wmp (:,:) = 0._wp ; 146 295 IF ( ln_zdfqiao ) THEN !== Vertical mixing enhancement using Qiao,2010 ==! 147 ALLOCATE( sf_wn(1), STAT=ierror ) !* allocate and fill sf_wave with sn_wnum 148 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave: unable toallocate sf_wave structure' ) 149 ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1) ) 150 IF( sn_wnum%ln_tint ) ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) 151 CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 296 IF ( .NOT. cpl_wnum ) THEN 297 ALLOCATE( sf_wn(1), STAT=ierror ) !* allocate and fill sf_wave with sn_wnum 298 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave: unable toallocate sf_wave structure' ) 299 ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1) ) 300 IF( sn_wnum%ln_tint ) ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) 301 CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' ) 302 ENDIF 152 303 ALLOCATE( wnum(jpi,jpj),tsd2d(jpi,jpj) ) 153 304 wnum(:,:) = 0._wp ; tsd2d(:,:) = 0._wp … … 156 307 ENDIF 157 308 ! 158 IF ( ln_cdgw ) THEN !== Neutral drag coefficient ==!309 IF ( ln_cdgw .AND. .NOT. cpl_wdrag ) THEN !== Neutral drag coefficient ==! 159 310 CALL fld_read( kt, nn_fsbc, sf_cd ) ! read from external forcing 160 311 cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 161 312 ENDIF 162 313 163 IF ( ln_tauoc ) THEN !== Wave induced stress ==!314 IF ( ln_tauoc .AND. .NOT. cpl_wstrf ) THEN !== Wave induced stress ==! 164 315 CALL fld_read( kt, nn_fsbc, sf_tauoc ) !* read wave norm stress from external forcing 165 316 tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1) … … 168 319 IF ( ln_sdw ) THEN !== Computation of the 3d Stokes Drift ==! 169 320 ! 170 CALL fld_read( kt, nn_fsbc, sf_sd ) !* read wave parameters from external forcing 171 swh(:,:) = sf_sd(jp_swh)%fnow(:,:,1) ! significant wave height 172 wmp(:,:) = sf_sd(jp_wmp)%fnow(:,:,1) ! wave mean period 173 zusd2dt(:,:) = sf_sd(jp_usd)%fnow(:,:,1) ! 2D zonal Stokes Drift at T point 174 zvsd2dt(:,:) = sf_sd(jp_vsd)%fnow(:,:,1) ! 2D meridional Stokes Drift at T point 321 ! Read from file only if the field is not coupled 322 IF( jpfld > 0 ) THEN 323 CALL fld_read( kt, nn_fsbc, sf_sd ) !* read wave parameters from external forcing 324 IF( jp_swh > 0 ) swh(:,:) = sf_sd(jp_swh)%fnow(:,:,1) ! significant wave height 325 IF( jp_wmp > 0 ) wmp(:,:) = sf_sd(jp_wmp)%fnow(:,:,1) ! wave mean period 326 IF( jp_usd > 0 ) zusd2dt(:,:) = sf_sd(jp_usd)%fnow(:,:,1) ! 2D zonal Stokes Drift at T point 327 IF( jp_vsd > 0 ) zvsd2dt(:,:) = sf_sd(jp_vsd)%fnow(:,:,1) ! 2D meridional Stokes Drift at T point 328 ENDIF 175 329 ! 330 ! Read also wave number if needed, so that it is available in coupling routines 331 IF ( ln_zdfqiao .AND. .NOT. cpl_wnum ) THEN 332 CALL fld_read( kt, nn_fsbc, sf_wn ) !* read wave parameters from external forcing 333 wnum(:,:) = sf_wn(1)%fnow(:,:,1) 334 ENDIF 335 176 336 !== Computation of the 3d Stokes Drift according to Breivik et al.,2014 177 337 !(DOI: 10.1175/JPO-D-14-0020.1)==! 178 338 ! 179 CALL wrk_alloc( jpi,jpj,jpk, ze3hdiv ) 180 DO jk = 1, jpk 181 DO jj = 1, jpj 182 DO ji = 1, jpi 183 ! On T grid 184 ! Stokes transport speed estimated from Hs and Tmean 185 ztransp = 2.0_wp*rpi*swh(ji,jj)**2.0_wp/(16.0_wp*MAX(wmp(ji,jj),0.0000001_wp)) 186 ! Stokes surface speed 187 zsp0 = SQRT( sf_sd(jp_usd)%fnow(ji,jj,1)**2 + sf_sd(jp_vsd)%fnow(ji,jj,1)**2) 188 ! Wavenumber scale 189 zk = ABS(zsp0)/MAX(ABS(5.97_wp*ztransp),0.0000001_wp) 190 ! Depth attenuation 191 zfac(ji,jj) = EXP(-2.0_wp*zk*fsdept(ji,jj,jk))/(1.0_wp+8.0_wp*zk*fsdept(ji,jj,jk)) 192 END DO 193 END DO 194 195 DO jj = 1, jpjm1 196 DO ji = 1, jpim1 197 ! Into the U and V Grid 198 zus = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zfac(ji,jj) * tmask(ji,jj,1) & 199 & + zfac(ji+1,jj) * tmask(ji+1,jj,1) ) 200 201 zvs = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zfac(ji,jj) * tmask(ji,jj,1) & 202 & + zfac(ji,jj+1) * tmask(ji,jj+1,1) ) 203 204 usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zusd2dt(ji,jj) * tmask(ji,jj,1) & 205 & + zusd2dt(ji+1,jj) * tmask(ji+1,jj,1) ) 206 207 vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zvsd2dt(ji,jj) * tmask(ji,jj,1) & 208 & + zvsd2dt(ji,jj+1) * tmask(ji,jj+1,1) ) 209 210 usd3d(ji,jj,jk) = usd2d(ji,jj) * zus 211 vsd3d(ji,jj,jk) = vsd2d(ji,jj) * zvs 212 END DO 213 END DO 214 END DO 215 ! 216 CALL lbc_lnk( usd3d(:,:,:), 'U', -1. ) 217 CALL lbc_lnk( vsd3d(:,:,:), 'V', -1. ) 218 ! 219 DO jk = 1, jpkm1 ! e3t * Horizontal divergence 220 DO jj = 2, jpjm1 221 DO ji = fs_2, fs_jpim1 ! vector opt. 222 ze3hdiv(ji,jj,jk) = ( e2u(ji ,jj) * fse3u_n(ji ,jj,jk) * usd3d(ji ,jj,jk) & 223 & - e2u(ji-1,jj) * fse3u_n(ji-1,jj,jk) * usd3d(ji-1,jj,jk) & 224 & + e1v(ji,jj ) * fse3v_n(ji,jj ,jk) * vsd3d(ji,jj ,jk) & 225 & - e1v(ji,jj-1) * fse3v_n(ji,jj-1,jk) * vsd3d(ji,jj-1,jk) ) * r1_e1e2t(ji,jj) 226 END DO 227 END DO 228 IF( .NOT. AGRIF_Root() ) THEN 229 IF( nbondi == 1 .OR. nbondi == 2 ) ze3hdiv(nlci-1, : ,jk) = 0._wp ! east 230 IF( nbondi == -1 .OR. nbondi == 2 ) ze3hdiv( 2 , : ,jk) = 0._wp ! west 231 IF( nbondj == 1 .OR. nbondj == 2 ) ze3hdiv( : ,nlcj-1,jk) = 0._wp ! north 232 IF( nbondj == -1 .OR. nbondj == 2 ) ze3hdiv( : , 2 ,jk) = 0._wp ! south 233 ENDIF 234 END DO 235 CALL lbc_lnk( ze3hdiv, 'T', 1. ) 236 ! 237 DO jk = jpkm1, 1, -1 !* integrate from the bottom the e3t * hor. divergence 238 wsd3d(:,:,jk) = wsd3d(:,:,jk+1) - ze3hdiv(:,:,jk) 239 END DO 240 #if defined key_bdy 241 IF( lk_bdy ) THEN 242 DO jk = 1, jpkm1 243 wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:) 244 END DO 245 ENDIF 246 #endif 247 CALL wrk_dealloc( jpi,jpj,jpk, ze3hdiv ) 248 249 IF ( ln_zdfqiao ) THEN 250 wnum(:,:) = sf_wn(1)%fnow(:,:,1) 251 252 ! Calculate the module of the stokes drift on T grid 253 !------------------------------------------------- 254 DO jj = 1, jpj 255 DO ji = 1, jpi 256 tsd2d(ji,jj) = ((sf_sd(jp_usd)%fnow(ji,jj,1) * tmask(ji,jj,1))**2.0 + & 257 & (sf_sd(jp_vsd)%fnow(ji,jj,1) * tmask(ji,jj,1))**2.0)**0.5 258 END DO 259 END DO 260 ENDIF 261 ! 339 ! Calculate only if no necessary fields are coupled, if not calculate later after coupling 340 IF( jpfld == 4 ) THEN 341 CALL sbc_stokes() 342 IF ( ln_zdfqiao .AND. .NOT. cpl_wnum ) THEN 343 CALL sbc_qiao() 344 ENDIF 345 ENDIF 262 346 ENDIF 263 347 ! -
branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/step.F90
r7167 r7168 231 231 CALL dyn_adv ( kstp ) ! advection (vector or flux form) 232 232 CALL dyn_vor ( kstp ) ! vorticity term including Coriolis 233 IF( ln_stcor ) CALL dyn_stcor ( kstp ) ! Stokes-Coriolis forcing 233 234 CALL dyn_ldf ( kstp ) ! lateral mixing 234 IF( ln_stcor ) CALL dyn_stcor ( kstp ) ! Stokes-Coriolis forcing235 235 CALL dyn_hpg ( kstp ) ! horizontal gradient of Hydrostatic pressure 236 236 CALL dyn_spg ( kstp ) ! surface pressure gradient
Note: See TracChangeset
for help on using the changeset viewer.
