- Timestamp:
- 2015-11-20T09:39:06+01:00 (9 years ago)
- Location:
- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC
- Files:
-
- 16 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90
r5620 r5901 4 4 !! Coupled O/A : coupled ocean-atmosphere case using OASIS3-MCT 5 5 !!===================================================================== 6 !! History : 7 !! 9.0 ! 04-06 (R. Redler, NEC Laboratories Europe, Germany) Original code 8 !! " " ! 04-11 (R. Redler, NEC Laboratories Europe; N. Keenlyside, W. Park, IFM-GEOMAR, Germany) revision 9 !! " " ! 04-11 (V. Gayler, MPI M&D) Grid writing 10 !! " " ! 05-08 (R. Redler, W. Park) frld initialization, paral(2) revision 11 !! " " ! 05-09 (R. Redler) extended to allow for communication over root only 12 !! " " ! 06-01 (W. Park) modification of physical part 13 !! " " ! 06-02 (R. Redler, W. Park) buffer array fix for root exchange 14 !! 3.4 ! 11-11 (C. Harris) Changes to allow mutiple category fields 15 !!---------------------------------------------------------------------- 6 !! History : 1.0 ! 2004-06 (R. Redler, NEC Laboratories Europe, Germany) Original code 7 !! - ! 2004-11 (R. Redler, NEC Laboratories Europe; N. Keenlyside, W. Park, IFM-GEOMAR, Germany) revision 8 !! - ! 2004-11 (V. Gayler, MPI M&D) Grid writing 9 !! 2.0 ! 2005-08 (R. Redler, W. Park) frld initialization, paral(2) revision 10 !! - ! 2005-09 (R. Redler) extended to allow for communication over root only 11 !! - ! 2006-01 (W. Park) modification of physical part 12 !! - ! 2006-02 (R. Redler, W. Park) buffer array fix for root exchange 13 !! 3.4 ! 2011-11 (C. Harris) Changes to allow mutiple category fields 14 !! 3.6 ! 2014-11 (S. Masson) OASIS3-MCT 15 !!---------------------------------------------------------------------- 16 16 17 !!---------------------------------------------------------------------- 17 18 !! 'key_oasis3' coupled Ocean/Atmosphere via OASIS3-MCT … … 20 21 !! cpl_init : initialization of coupled mode communication 21 22 !! cpl_define : definition of grid and fields 22 !! cpl_snd : snd out fields in coupled mode23 !! cpl_rcv : receive fields in coupled mode23 !! cpl_snd : snd out fields in coupled mode 24 !! cpl_rcv : receive fields in coupled mode 24 25 !! cpl_finalize : finalize the coupled mode communication 25 26 !!---------------------------------------------------------------------- … … 99 100 !! ** Method : OASIS3 MPI communication 100 101 !!-------------------------------------------------------------------- 101 CHARACTER(len = *), INTENT(in ) :: cd_modname ! model name as set in namcouple file102 INTEGER , INTENT(out) :: kl_comm ! local communicator of the model102 CHARACTER(len = *), INTENT(in ) :: cd_modname ! model name as set in namcouple file 103 INTEGER , INTENT( out) :: kl_comm ! local communicator of the model 103 104 !!-------------------------------------------------------------------- 104 105 … … 163 164 CALL oasis_abort ( ncomp_id, 'cpl_define', 'nsnd is larger than nmaxfld, increase nmaxfld') ; RETURN 164 165 ENDIF 165 166 166 ! 167 167 ! ... Define the shape for the area that excludes the halo -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90
r4162 r5901 195 195 196 196 DO jj = 2, jpjm1 197 !CDIR NOVERRCHK198 197 DO ji = fs_2, jpi ! vector opt. 199 198 -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
r5620 r5901 80 80 INTEGER , PUBLIC, PARAMETER :: jp_mfs = 6 !: MFS bulk formulation 81 81 INTEGER , PUBLIC, PARAMETER :: jp_none = 7 !: for OPA when doing coupling via SAS module 82 INTEGER , PUBLIC, PARAMETER :: jp_esopa = -1 !: esopa test, ALL formulations83 82 84 83 !!---------------------------------------------------------------------- … … 200 199 !!--------------------------------------------------------------------- 201 200 zcoef = 0.5 / ( zrhoa * zcdrag ) 202 !CDIR NOVERRCHK203 201 DO jj = 2, jpjm1 204 !CDIR NOVERRCHK205 202 DO ji = fs_2, fs_jpim1 ! vect. opt. 206 203 ztx = utau(ji-1,jj ) + utau(ji,jj) -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90
r4792 r5901 279 279 ! module of wind stress and wind speed at T-point 280 280 zcoef = 1. / ( zrhoa * zcdrag ) 281 !CDIR NOVERRCHK282 281 DO jj = 2, jpjm1 283 !CDIR NOVERRCHK284 282 DO ji = fs_2, fs_jpim1 ! vect. opt. 285 283 ztx = utau(ji-1,jj ) + utau(ji,jj) -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90
r5620 r5901 62 62 !!--------------------------------------------------------------------- 63 63 INTEGER, INTENT(in):: kt ! ocean time step 64 ! !64 ! 65 65 INTEGER :: ierror ! local integer 66 66 INTEGER :: ios ! Local integer output status for namelist read … … 71 71 NAMELIST/namsbc_apr/ cn_dir, sn_apr, ln_ref_apr, rn_pref, ln_apr_obc 72 72 !!---------------------------------------------------------------------- 73 !74 73 ! 75 74 ! ! -------------------- ! -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90
r5620 r5901 243 243 ! momentum fluxes (utau, vtau ) ! 244 244 !------------------------------------! 245 !CDIR COLLAPSE246 245 utau(:,:) = sf(jp_utau)%fnow(:,:,1) 247 !CDIR COLLAPSE248 246 vtau(:,:) = sf(jp_vtau)%fnow(:,:,1) 249 247 … … 251 249 ! wind stress module (taum ) ! 252 250 !------------------------------------! 253 !CDIR NOVERRCHK254 251 DO jj = 2, jpjm1 255 !CDIR NOVERRCHK256 252 DO ji = fs_2, fs_jpim1 ! vector opt. 257 253 ztx2 = utau(ji-1,jj ) + utau(ji,jj) … … 268 264 ! store the wind speed (wndm ) ! 269 265 !------------------------------------! 270 !CDIR COLLAPSE271 266 wndm(:,:) = sf(jp_wndm)%fnow(:,:,1) 272 267 wndm(:,:) = wndm(:,:) * tmask(:,:,1) … … 281 276 ! Other ocean fluxes ! 282 277 !------------------------! 283 !CDIR NOVERRCHK284 !CDIR COLLAPSE285 278 DO jj = 1, jpj 286 !CDIR NOVERRCHK287 279 DO ji = 1, jpi 288 280 ! … … 375 367 zcprec = rcp / rday ! convert prec ( mm/day ==> m/s) ==> W/m2 376 368 377 !CDIR COLLAPSE378 369 emp(:,:) = zqla(:,:) / cevap & ! freshwater flux 379 370 & - sf(jp_prec)%fnow(:,:,1) / rday * tmask(:,:,1) 380 371 ! 381 !CDIR COLLAPSE382 372 qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) & ! Downward Non Solar flux 383 373 & - zqla(:,:) * pst(:,:) * zcevap & ! remove evap. heat content at SST in Celcius … … 415 405 416 406 # if defined key_lim2 || defined key_lim3 407 417 408 SUBROUTINE blk_ice_clio_tau 418 409 !!--------------------------------------------------------------------------- … … 429 420 ! 430 421 IF( nn_timing == 1 ) CALL timing_start('blk_ice_clio_tau') 431 422 ! 432 423 SELECT CASE( cp_ice_msh ) 433 424 ! 434 425 CASE( 'C' ) ! C-grid ice dynamics 435 426 ! 436 427 zcoef = cai / cao ! Change from air-sea stress to air-ice stress 437 428 utau_ice(:,:) = zcoef * utau(:,:) 438 429 vtau_ice(:,:) = zcoef * vtau(:,:) 439 430 ! 440 431 CASE( 'I' ) ! I-grid ice dynamics: I-point (i.e. F-point lower-left corner) 441 432 ! 442 433 zcoef = 0.5_wp * cai / cao ! Change from air-sea stress to air-ice stress 443 434 DO jj = 2, jpj ! stress from ocean U- and V-points to ice U,V point … … 447 438 END DO 448 439 END DO 449 440 ! 450 441 CALL lbc_lnk( utau_ice(:,:), 'I', -1. ) ; CALL lbc_lnk( vtau_ice(:,:), 'I', -1. ) ! I-point 451 442 ! 452 443 END SELECT 453 444 ! 454 445 IF(ln_ctl) THEN 455 446 CALL prt_ctl(tab2d_1=utau_ice , clinfo1=' blk_ice_clio: utau_ice : ', tab2d_2=vtau_ice , clinfo2=' vtau_ice : ') 456 447 ENDIF 457 448 ! 458 449 IF( nn_timing == 1 ) CALL timing_stop('blk_ice_clio_tau') 459 450 ! 460 451 END SUBROUTINE blk_ice_clio_tau 452 461 453 #endif 462 454 463 455 # if defined key_lim2 || defined key_lim3 456 464 457 SUBROUTINE blk_ice_clio_flx( ptsu , palb_cs, palb_os, palb ) 465 458 !!--------------------------------------------------------------------------- … … 520 513 !-------------------------------------------------------------------------------- 521 514 522 !CDIR NOVERRCHK523 !CDIR COLLAPSE524 515 DO jj = 1, jpj 525 !CDIR NOVERRCHK526 516 DO ji = 1, jpi 527 517 ztatm (ji,jj) = sf(jp_tair)%fnow(ji,jj,1) ! air temperature in Kelvins … … 573 563 574 564 ! ! ========================== ! 575 DO jl = 1, jpl ! Loop over ice categories !565 DO jl = 1, jpl ! Loop over ice categories ! 576 566 ! ! ========================== ! 577 !CDIR NOVERRCHK578 !CDIR COLLAPSE579 567 DO jj = 1 , jpj 580 !CDIR NOVERRCHK581 568 DO ji = 1, jpi 582 569 !-------------------------------------------! … … 636 623 ! ----------------------------------------------------------------------------- ! 637 624 ! 638 !CDIR COLLAPSE639 625 qns_ice(:,:,:) = z_qlw (:,:,:) - z_qsb (:,:,:) - qla_ice (:,:,:) ! Downward Non Solar flux 640 !CDIR COLLAPSE641 626 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) / rday ! total precipitation [kg/m2/s] 642 627 ! … … 644 629 ! Correct the OCEAN non solar flux with the existence of solid precipitation ! 645 630 ! ---------------=====--------------------------------------------------------- ! 646 !CDIR COLLAPSE647 631 qns(:,:) = qns(:,:) & ! update the non-solar heat flux with: 648 632 & - sprecip(:,:) * lfus & ! remove melting solid precip … … 782 766 ! Saturated water vapour and vapour pressure 783 767 ! ------------------------------------------ 784 !CDIR NOVERRCHK785 !CDIR COLLAPSE786 768 DO jj = 1, jpj 787 !CDIR NOVERRCHK788 769 DO ji = 1, jpi 789 770 ztamr = sf(jp_tair)%fnow(ji,jj,1) - rtt … … 814 795 zdaycor = 1.0 + 0.0013 * SIN( zdist ) + 0.0342 * COS( zdist ) 815 796 816 !CDIR NOVERRCHK817 797 DO jj = 1, jpj 818 !CDIR NOVERRCHK819 798 DO ji = 1, jpi 820 799 ! product of sine (cosine) of latitude and sine (cosine) of solar declination … … 837 816 838 817 ! compute and sum ocean qsr over the daylight (i.e. between sunrise and sunset) 839 !CDIR NOVERRCHK840 818 DO jt = 1, jp24 841 819 zcoef = FLOAT( jt ) - 0.5 842 !CDIR NOVERRCHK843 !CDIR COLLAPSE844 820 DO jj = 1, jpj 845 !CDIR NOVERRCHK846 821 DO ji = 1, jpi 847 822 zlha = COS( zlsrise(ji,jj) - zcoef * zdlha(ji,jj) ) ! local hour angle … … 862 837 ! Taking into account the ellipsity of the earth orbit, the clouds AND masked if sea-ice cover > 0% 863 838 zcoef1 = srgamma * zdaycor / ( 2. * rpi ) 864 !CDIR COLLAPSE865 839 DO jj = 1, jpj 866 840 DO ji = 1, jpi … … 920 894 ! Saturated water vapour and vapour pressure 921 895 ! ------------------------------------------ 922 !CDIR NOVERRCHK923 !CDIR COLLAPSE924 896 DO jj = 1, jpj 925 !CDIR NOVERRCHK926 897 DO ji = 1, jpi 927 898 ztamr = sf(jp_tair)%fnow(ji,jj,1) - rtt … … 952 923 zdaycor = 1.0 + 0.0013 * SIN( zdist ) + 0.0342 * COS( zdist ) 953 924 954 !CDIR NOVERRCHK955 925 DO jj = 1, jpj 956 !CDIR NOVERRCHK957 926 DO ji = 1, jpi 958 927 ! product of sine (cosine) of latitude and sine (cosine) of solar declination … … 979 948 DO jl = 1, ijpl ! loop over ice categories ! 980 949 ! !----------------------------! 981 !CDIR NOVERRCHK982 950 DO jt = 1, jp24 983 951 zcoef = FLOAT( jt ) - 0.5 984 !CDIR NOVERRCHK985 !CDIR COLLAPSE986 952 DO jj = 1, jpj 987 !CDIR NOVERRCHK988 953 DO ji = 1, jpi 989 954 zlha = COS( zlsrise(ji,jj) - zcoef * zdlha(ji,jj) ) ! local hour angle -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90
r5620 r5901 251 251 ! for basin budget and cooerence 252 252 !-------------------------------------------------- 253 !CDIR COLLAPSE 254 emp (:,:) = evap(:,:) - sf(jp_prec)%fnow(:,:,1) * tmask(:,:,1) 255 !CDIR COLLAPSE 253 emp(:,:) = evap(:,:) - sf(jp_prec)%fnow(:,:,1) * tmask(:,:,1) 256 254 257 255 CALL iom_put( "qlw_oce", qbw ) ! output downward longwave heat over the ocean -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
r5620 r5901 931 931 ! => need to be done only when otx1 was changed 932 932 IF( llnewtx ) THEN 933 !CDIR NOVERRCHK934 933 DO jj = 2, jpjm1 935 !CDIR NOVERRCHK936 934 DO ji = fs_2, fs_jpim1 ! vect. opt. 937 935 zzx = frcv(jpr_otx1)%z3(ji-1,jj ,1) + frcv(jpr_otx1)%z3(ji,jj,1) … … 961 959 IF( llnewtau ) THEN 962 960 zcoef = 1. / ( zrhoa * zcdrag ) 963 !CDIR NOVERRCHK964 961 DO jj = 1, jpj 965 !CDIR NOVERRCHK966 962 DO ji = 1, jpi 967 963 frcv(jpr_w10m)%z3(ji,jj,1) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef ) -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90
r5038 r5901 131 131 ELSE ; qsr(:,:) = sf(jp_qsr)%fnow(:,:,1) 132 132 ENDIF 133 !CDIR COLLAPSE134 133 DO jj = 1, jpj ! set the ocean fluxes from read fields 135 134 DO ji = 1, jpi … … 145 144 ! ! module of wind stress and wind speed at T-point 146 145 zcoef = 1. / ( zrhoa * zcdrag ) 147 !CDIR NOVERRCHK148 146 DO jj = 2, jpjm1 149 !CDIR NOVERRCHK150 147 DO ji = fs_2, fs_jpim1 ! vect. opt. 151 148 ztx = utau(ji-1,jj ) + utau(ji,jj) -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90
r5620 r5901 108 108 ! 109 109 IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN 110 z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) -snwice_fmass(:,:) ) ) / area ! sum over the global domain110 z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) - snwice_fmass(:,:) ) ) / area ! sum over the global domain 111 111 zcoef = z_fwf * rcp 112 112 emp(:,:) = emp(:,:) - z_fwf * tmask(:,:,1) … … 162 162 zsurf_pos = glob_sum( e1e2t(:,:)*ztmsk_pos(:,:) ) 163 163 ! ! fwf global mean (excluding ocean to ice/snow exchanges) 164 z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + rdivisf *fwfisf(:,:) - snwice_fmass(:,:) ) ) / area164 z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) - snwice_fmass(:,:) ) ) / area 165 165 ! 166 166 IF( z_fwf < 0._wp ) THEN ! spread out over >0 erp area to increase evaporation -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90
r5620 r5901 67 67 PRIVATE 68 68 69 !! * Routine accessibility70 69 PUBLIC cice_sbc_init ! routine called by sbc_init 71 70 PUBLIC cice_sbc_final ! routine called by sbc_final … … 95 94 !! * Substitutions 96 95 # include "domzgr_substitute.h90" 97 96 !!---------------------------------------------------------------------- 97 !! NEMO/OPA 3.7 , NEMO-consortium (2015) 98 98 !! $Id$ 99 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 100 !!---------------------------------------------------------------------- 99 101 CONTAINS 100 102 … … 154 156 END SUBROUTINE sbc_ice_cice 155 157 156 SUBROUTINE cice_sbc_init (ksbc) 158 159 SUBROUTINE cice_sbc_init( ksbc ) 157 160 !!--------------------------------------------------------------------- 158 161 !! *** ROUTINE cice_sbc_init *** 159 162 !! ** Purpose: Initialise ice related fields for NEMO and coupling 160 163 !! 164 !!--------------------------------------------------------------------- 161 165 INTEGER, INTENT( in ) :: ksbc ! surface forcing type 162 166 REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2 … … 289 293 290 294 291 SUBROUTINE cice_sbc_in (kt, ksbc)295 SUBROUTINE cice_sbc_in( kt, ksbc ) 292 296 !!--------------------------------------------------------------------- 293 297 !! *** ROUTINE cice_sbc_in *** … … 296 300 INTEGER, INTENT(in ) :: kt ! ocean time step 297 301 INTEGER, INTENT(in ) :: ksbc ! surface forcing type 298 302 ! 299 303 INTEGER :: ji, jj, jl ! dummy loop indices 300 304 REAL(wp), DIMENSION(:,:), POINTER :: ztmp, zpice … … 490 494 ! x comp and y comp of sea surface slope (on F points) 491 495 ! T point to F point 492 DO jj=1,jpjm1 493 DO ji=1,jpim1 494 ztmp(ji,jj)=0.5 * ( (zpice(ji+1,jj )-zpice(ji,jj ))/e1u(ji,jj ) & 495 + (zpice(ji+1,jj+1)-zpice(ji,jj+1))/e1u(ji,jj+1) ) & 496 * fmask(ji,jj,1) 497 ENDDO 498 ENDDO 499 CALL nemo2cice(ztmp,ss_tltx,'F', -1. ) 496 DO jj = 1, jpjm1 497 DO ji = 1, jpim1 498 ztmp(ji,jj)=0.5 * ( (zpice(ji+1,jj )-zpice(ji,jj )) * r1_e1u(ji,jj ) & 499 & + (zpice(ji+1,jj+1)-zpice(ji,jj+1)) * r1_e1u(ji,jj+1) ) * fmask(ji,jj,1) 500 END DO 501 END DO 502 CALL nemo2cice( ztmp,ss_tltx,'F', -1. ) 500 503 501 504 ! T point to F point 502 DO jj=1,jpjm1 503 DO ji=1,jpim1 504 ztmp(ji,jj)=0.5 * ( (zpice(ji ,jj+1)-zpice(ji ,jj))/e2v(ji ,jj) & 505 + (zpice(ji+1,jj+1)-zpice(ji+1,jj))/e2v(ji+1,jj) ) & 506 * fmask(ji,jj,1) 507 ENDDO 508 ENDDO 505 DO jj = 1, jpjm1 506 DO ji = 1, jpim1 507 ztmp(ji,jj)=0.5 * ( (zpice(ji ,jj+1)-zpice(ji ,jj)) * r1_e2v(ji ,jj) & 508 & + (zpice(ji+1,jj+1)-zpice(ji+1,jj)) * r1_e2v(ji+1,jj) ) * fmask(ji,jj,1) 509 END DO 510 END DO 509 511 CALL nemo2cice(ztmp,ss_tlty,'F', -1. ) 510 512 … … 517 519 518 520 519 SUBROUTINE cice_sbc_out (kt,ksbc)521 SUBROUTINE cice_sbc_out( kt, ksbc ) 520 522 !!--------------------------------------------------------------------- 521 523 !! *** ROUTINE cice_sbc_out *** … … 575 577 ! Update taum with modulus of ice-ocean stress 576 578 ! strocnxT and strocnyT are not weighted by ice fraction in CICE so must be done here 577 taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1* *2. + ztmp2**2.)579 taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2) 578 580 579 581 ! Freshwater fluxes … … 888 890 #endif 889 891 !!--------------------------------------------------------------------- 890 891 892 CHARACTER(len=1), INTENT( in ) :: & 892 893 cd_type ! nature of pn grid-point … … 908 909 909 910 INTEGER :: ji, jj, jn ! dummy loop indices 911 !!--------------------------------------------------------------------- 910 912 911 913 ! A. Ensure all haloes are filled in NEMO field (pn) … … 1096 1098 !! Default option Dummy module NO CICE sea-ice model 1097 1099 !!---------------------------------------------------------------------- 1098 !! $Id$1099 1100 CONTAINS 1100 1101 -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
r5620 r5901 18 18 USE eosbn2 ! equation of state 19 19 USE sbc_oce ! surface boundary condition: ocean fields 20 USE zdfbfr ! 21 ! 22 USE in_out_manager ! I/O manager 23 USE iom ! I/O manager library 24 USE fldread ! read input field at current time step 20 25 USE lbclnk ! 21 USE iom ! I/O manager library22 USE in_out_manager ! I/O manager23 26 USE wrk_nemo ! Memory allocation 24 27 USE timing ! Timing 25 28 USE lib_fortran ! glob_sum 26 USE zdfbfr27 USE fldread ! read input field at current time step28 29 30 29 31 30 IMPLICIT NONE 32 31 PRIVATE 33 32 34 PUBLIC sbc_isf, sbc_isf_div, sbc_isf_alloc ! routine called in sbcmod and div cur33 PUBLIC sbc_isf, sbc_isf_div, sbc_isf_alloc ! routine called in sbcmod and divhor 35 34 36 35 ! public in order to be able to output then … … 53 52 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: risfLeff !:effective length (Leff) BG03 nn_isf==2 54 53 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: ttbl, stbl, utbl, vtbl !:top boundary layer variable at T point 55 #if defined key_agrif56 ! AGRIF can not handle these arrays as integers. The reason is a mystery but problems avoided by declaring them as reals57 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: misfkt, misfkb !:Level of ice shelf base58 !: (first wet level and last level include in the tbl)59 #else60 54 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: misfkt, misfkb !:Level of ice shelf base 61 #endif62 63 55 64 56 REAL(wp), PUBLIC, SAVE :: rcpi = 2000.0_wp ! phycst ? … … 79 71 # include "domzgr_substitute.h90" 80 72 !!---------------------------------------------------------------------- 81 !! NEMO/OPA 3. 0 , LOCEAN-IPSL (2008)73 !! NEMO/OPA 3.7 , LOCEAN-IPSL (2015) 82 74 !! $Id$ 83 75 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 84 76 !!---------------------------------------------------------------------- 85 86 77 CONTAINS 87 78 88 SUBROUTINE sbc_isf(kt) 89 INTEGER, INTENT(in) :: kt ! ocean time step 90 INTEGER :: ji, jj, jk, ijkmin, inum, ierror 91 INTEGER :: ikt, ikb ! top and bottom level of the isf boundary layer 92 REAL(wp) :: rmin 93 REAL(wp) :: zhk 94 CHARACTER(len=256) :: cfisf, cvarzisf, cvarhisf ! name for isf file 95 CHARACTER(LEN=256) :: cnameis ! name of iceshelf file 96 CHARACTER (LEN=32) :: cvarLeff ! variable name for efficient Length scale 97 INTEGER :: ios ! Local integer output status for namelist read 98 ! 79 SUBROUTINE sbc_isf(kt) 99 80 !!--------------------------------------------------------------------- 81 !! *** ROUTINE sbc_isf *** 82 !!--------------------------------------------------------------------- 83 INTEGER, INTENT(in) :: kt ! ocean time step 84 ! 85 INTEGER :: ji, jj, jk, ijkmin, inum, ierror 86 INTEGER :: ikt, ikb ! top and bottom level of the isf boundary layer 87 REAL(wp) :: rmin 88 REAL(wp) :: zhk 89 REAL(wp) :: zt_frz, zpress 90 CHARACTER(len=256) :: cfisf , cvarzisf, cvarhisf ! name for isf file 91 CHARACTER(LEN=256) :: cnameis ! name of iceshelf file 92 CHARACTER (LEN=32) :: cvarLeff ! variable name for efficient Length scale 93 INTEGER :: ios ! Local integer output status for namelist read 94 !! 100 95 NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, ln_divisf, ln_conserve, rn_gammat0, rn_gammas0, nn_gammablk, & 101 &sn_fwfisf, sn_qisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf102 ! 96 & sn_fwfisf, sn_qisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf 97 !!--------------------------------------------------------------------- 103 98 ! 104 99 ! ! ====================== ! … … 113 108 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp ) 114 109 IF(lwm) WRITE ( numond, namsbc_isf ) 115 116 110 117 111 IF ( lwp ) WRITE(numout,*) … … 194 188 END IF 195 189 190 ! compute bottom level of isf tbl and thickness of tbl below the ice shelf 196 191 rhisf_tbl_0(:,:) = rhisf_tbl(:,:) 197 198 ! compute bottom level of isf tbl and thickness of tbl below the ice shelf199 192 DO jj = 1,jpj 200 193 DO ji = 1,jpi … … 217 210 END DO 218 211 END DO 219 212 ! 220 213 END IF 221 214 … … 270 263 END IF 271 264 ! compute tsc due to isf 272 ! WARNING water add at temp = 0C, correction term is added in trasbc, maybe better here but need a 3D variable). 273 risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp ! 265 ! WARNING water add at temp = 0C, correction term is added, maybe better here but need a 3D variable). 266 ! zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04 267 zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress ) 268 risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - rdivisf * fwfisf(:,:) * zt_frz * r1_rau0 ! 274 269 275 270 ! salt effect already take into account in vertical advection 276 271 risf_tsc(:,:,jp_sal) = (1.0_wp-rdivisf) * fwfisf(:,:) * stbl(:,:) * r1_rau0 277 272 273 ! output 274 IF( iom_use('qisf' ) ) CALL iom_put('qisf' , qisf) 275 IF( iom_use('fwfisf') ) CALL iom_put('fwfisf', fwfisf * stbl(:,:) / soce ) 276 277 ! if apply only on the trend and not as a volume flux (rdivisf = 0), fwfisf have to be set to 0 now 278 fwfisf(:,:) = rdivisf * fwfisf(:,:) 279 278 280 ! lbclnk 279 281 CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.) … … 295 297 ENDIF 296 298 ! 297 ! output298 CALL iom_put('qisf' , qisf)299 IF( iom_use('fwfisf') ) CALL iom_put('fwfisf', fwfisf * stbl(:,:) / soce )300 299 END IF 301 300 ! 302 301 END SUBROUTINE sbc_isf 302 303 303 304 304 INTEGER FUNCTION sbc_isf_alloc() … … 321 321 END FUNCTION 322 322 323 SUBROUTINE sbc_isf_bg03(kt) 324 !!==========================================================================325 !! *** SUBROUTINE sbcisf_bg03 ***326 !! add net heat and fresh water flux from ice shelf melting327 !! into the adjacent ocean using the parameterisation by328 !! Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean329 !! interaction for climate models", Ocean Modelling 5(2003) 157-170.330 !! (hereafter BG)331 !!==========================================================================332 !!----------------------------------------------------------------------333 !! sbc_isf_bg03 : routine called from sbcmod334 !!----------------------------------------------------------------------335 !!336 !! ** Purpose : Add heat and fresh water fluxes due to ice shelf melting337 !! ** Reference : Beckmann et Goosse, 2003, Ocean Modelling338 !!339 !! History :340 !! ! 06-02 (C. Wang) Original code341 !!----------------------------------------------------------------------342 343 INTEGER, INTENT ( in ) :: kt344 323 324 SUBROUTINE sbc_isf_bg03(kt) 325 !!========================================================================== 326 !! *** SUBROUTINE sbcisf_bg03 *** 327 !! add net heat and fresh water flux from ice shelf melting 328 !! into the adjacent ocean using the parameterisation by 329 !! Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean 330 !! interaction for climate models", Ocean Modelling 5(2003) 157-170. 331 !! (hereafter BG) 332 !!========================================================================== 333 !!---------------------------------------------------------------------- 334 !! sbc_isf_bg03 : routine called from sbcmod 335 !!---------------------------------------------------------------------- 336 !! 337 !! ** Purpose : Add heat and fresh water fluxes due to ice shelf melting 338 !! ** Reference : Beckmann et Goosse, 2003, Ocean Modelling 339 !! 340 !! History : 341 !! ! 06-02 (C. Wang) Original code 342 !!---------------------------------------------------------------------- 343 INTEGER, INTENT ( in ) :: kt 344 ! 345 345 INTEGER :: ji, jj, jk, jish !temporary integer 346 346 INTEGER :: ijkmin … … 386 386 qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp 387 387 END IF 388 END DO389 END DO388 END DO 389 END DO 390 390 ! 391 391 IF( nn_timing == 1 ) CALL timing_stop('sbc_isf_bg03') 392 ! 392 393 END SUBROUTINE sbc_isf_bg03 394 393 395 394 396 SUBROUTINE sbc_isf_cav( kt ) … … 439 441 ! 440 442 ! 441 !CDIR COLLAPSE442 443 DO jj = 1, jpj 443 444 DO ji = 1, jpi … … 472 473 473 474 nit = nit + 1 474 IF (nit .GE. 100) THEN 475 !WRITE(numout,*) "sbcisf : too many iteration ... ", zhtflx, zhtflx_b,zgammat, rn_gammat0, rn_tfri2, nn_gammablk, ji,jj 476 !WRITE(numout,*) "sbcisf : too many iteration ... ", (zhtflx - zhtflx_b)/zhtflx 477 CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 478 END IF 475 IF (nit .GE. 100) CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 476 479 477 ! save gammat and compute zhtflx_b 480 478 zgammat2d(ji,jj)=zgammat … … 496 494 497 495 ! More complicated 3 equation thermodynamics as in MITgcm 498 !CDIR COLLAPSE499 496 DO jj = 2, jpj 500 497 DO ji = 2, jpi … … 565 562 ! 566 563 IF( nn_timing == 1 ) CALL timing_stop('sbc_isf_cav') 567 564 ! 568 565 END SUBROUTINE sbc_isf_cav 566 569 567 570 568 SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf, ji, jj, lit ) … … 693 691 END IF 694 692 END IF 695 693 ! 696 694 END SUBROUTINE 695 697 696 698 697 SUBROUTINE sbc_isf_tbl( varin, varout, cptin ) … … 756 755 IF (cptin == 'T') CALL lbc_lnk(varout,'T',1.) 757 756 IF (cptin == 'U' .OR. cptin == 'V') CALL lbc_lnk(varout,'T',-1.) 758 757 ! 759 758 END SUBROUTINE sbc_isf_tbl 760 759 … … 794 793 ! test on tmask useless ????? 795 794 DO jk = ikt, mbkt(ji,jj) 796 !IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk795 IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 797 796 END DO 798 797 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. … … 823 822 ! 824 823 END SUBROUTINE sbc_isf_div 825 824 825 826 826 FUNCTION tinsitu( ptem, psal, ppress ) RESULT( pti ) 827 827 !!---------------------------------------------------------------------- … … 874 874 ! 875 875 END FUNCTION tinsitu 876 ! 876 877 877 878 FUNCTION fsatg( pfps, pfpt, pfphp ) 878 879 !!---------------------------------------------------------------------- -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
r5620 r5901 28 28 USE sbcdcy ! surface boundary condition: diurnal cycle 29 29 USE sbcssm ! surface boundary condition: sea-surface mean variables 30 USE sbcapr ! surface boundary condition: atmospheric pressure31 30 USE sbcana ! surface boundary condition: analytical formulation 32 31 USE sbcflx ! surface boundary condition: flux formulation … … 133 132 WRITE(numout,*) ' forced-coupled mixed formulation ln_mixcpl = ', ln_mixcpl 134 133 WRITE(numout,*) ' OASIS coupling (with atm or sas) lk_oasis = ', lk_oasis 135 WRITE(numout,*) ' components of your executable 134 WRITE(numout,*) ' components of your executable nn_components = ', nn_components 136 135 WRITE(numout,*) ' Multicategory heat flux formulation (LIM3) nn_limflx = ', nn_limflx 137 136 WRITE(numout,*) ' Misc. options of sbc : ' … … 176 175 177 176 ! ! allocate sbc arrays 178 IF( sbc_oce_alloc() /= 0 ) CALL ctl_stop( ' STOP', 'sbc_init : unable to allocate sbc_oce arrays' )177 IF( sbc_oce_alloc() /= 0 ) CALL ctl_stop( 'sbc_init : unable to allocate sbc_oce arrays' ) 179 178 180 179 ! ! Checks: 181 180 IF( nn_isf .EQ. 0 ) THEN ! variable initialisation if no ice shelf 182 IF( sbc_isf_alloc() /= 0 ) CALL ctl_stop( ' STOP', 'sbc_init : unable to allocate sbc_isf arrays' )183 fwfisf (:,:) = 0.0_wp184 fwfisf_b(:,:) = 0.0_wp181 IF( sbc_isf_alloc() /= 0 ) CALL ctl_stop( 'sbc_init : unable to allocate sbc_isf arrays' ) 182 fwfisf (:,:) = 0.0_wp ; fwfisf_b (:,:) = 0.0_wp 183 risf_tsc(:,:,:) = 0.0_wp ; risf_tsc_b(:,:,:) = 0.0_wp 185 184 rdivisf = 0.0_wp 186 185 END IF … … 224 223 ENDIF 225 224 ELSE 226 IF ( ln_cdgw .OR. ln_sdw ) &227 & CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but&228 & asked couplingwith drag coefficient (ln_cdgw =T) or Stokes drift (ln_sdw=T) ')225 IF ( ln_cdgw .OR. ln_sdw ) & 226 & CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ', & 227 & 'with drag coefficient (ln_cdgw =T) or Stokes drift (ln_sdw=T) ') 229 228 ENDIF 230 229 ! ! Choice of the Surface Boudary Condition (set nsbc) … … 241 240 IF( nn_components == jp_iam_opa ) & 242 241 & THEN ; nsbc = jp_none ; icpt = icpt + 1 ; ENDIF ! opa coupling via SAS module 243 IF( lk_esopa ) nsbc = jp_esopa ! esopa test, ALL formulations 244 ! 245 IF( icpt /= 1 .AND. .NOT.lk_esopa ) THEN 246 WRITE(numout,*) 247 WRITE(numout,*) ' E R R O R in setting the sbc, one and only one namelist/CPP key option ' 248 WRITE(numout,*) ' must be choosen. You choose ', icpt, ' option(s)' 249 WRITE(numout,*) ' We stop' 250 nstop = nstop + 1 251 ENDIF 242 ! 243 IF( icpt /= 1 ) CALL ctl_stop( 'sbc_init: choose ONE and only ONE sbc option' ) 244 ! 252 245 IF(lwp) THEN 253 246 WRITE(numout,*) 254 IF( nsbc == jp_esopa ) WRITE(numout,*) ' ESOPA test All surface boundary conditions'255 247 IF( nsbc == jp_gyre ) WRITE(numout,*) ' GYRE analytical formulation' 256 248 IF( nsbc == jp_ana ) WRITE(numout,*) ' analytical formulation' … … 267 259 ! 268 260 IF( lk_oasis ) CALL sbc_cpl_init (nn_ice) ! OASIS initialisation. must be done before: (1) first time step 269 ! 261 ! ! (2) the use of nn_fsbc 270 262 271 263 ! nn_fsbc initialization if OPA-SAS coupling via OASIS 272 264 ! sas model time step has to be declared in OASIS (mandatory) -> nn_fsbc has to be modified accordingly 273 265 IF ( nn_components /= jp_iam_nemo ) THEN 274 275 266 IF ( nn_components == jp_iam_opa ) nn_fsbc = cpl_freq('O_SFLX') / NINT(rdt) 276 267 IF ( nn_components == jp_iam_sas ) nn_fsbc = cpl_freq('I_SFLX') / NINT(rdt) … … 344 335 ! ! forcing field computation ! 345 336 ! ! ---------------------------------------- ! 346 !347 IF ( .NOT. lk_bdy ) then348 IF( ln_apr_dyn ) CALL sbc_apr( kt ) ! atmospheric pressure provided at kt+0.5*nn_fsbc349 ENDIF350 ! (caution called before sbc_ssm)351 337 ! 352 338 IF( nn_components /= jp_iam_sas ) CALL sbc_ssm( kt ) ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m) … … 373 359 IF( nn_components == jp_iam_opa ) & 374 360 CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice ) ! OPA-SAS coupling: OPA receiving fields from SAS 375 CASE( jp_esopa )376 CALL sbc_ana ( kt ) ! ESOPA, test ALL the formulations377 CALL sbc_gyre ( kt ) !378 CALL sbc_flx ( kt ) !379 CALL sbc_blk_clio( kt ) !380 CALL sbc_blk_core( kt ) !381 CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice ) !382 361 END SELECT 383 362 -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
r5620 r5901 31 31 PRIVATE 32 32 33 PUBLIC sbc_rnf ! routine call in sbcmod module 34 PUBLIC sbc_rnf_div ! routine called in divcurl module 35 PUBLIC sbc_rnf_alloc ! routine call in sbcmod module 36 PUBLIC sbc_rnf_init ! (PUBLIC for TAM) 33 PUBLIC sbc_rnf ! routine called in sbcmod module 34 PUBLIC sbc_rnf_div ! routine called in divhor module 35 PUBLIC sbc_rnf_alloc ! routine called in sbcmod module 36 PUBLIC sbc_rnf_init ! routine called in sbcmod module 37 37 38 ! !!* namsbc_rnf namelist * 38 39 CHARACTER(len=100) :: cn_dir !: Root directory for location of rnf files -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90
r5038 r5901 107 107 IF( nn_sssr == 1 ) THEN !* Salinity damping term (salt flux only (sfx)) 108 108 zsrp = rn_deds / rday ! from [mm/day] to [kg/m2/s] 109 !CDIR COLLAPSE110 109 DO jj = 1, jpj 111 110 DO ji = 1, jpi … … 121 120 zsrp = rn_deds / rday ! from [mm/day] to [kg/m2/s] 122 121 zerp_bnd = rn_sssr_bnd / rday ! - - 123 !CDIR COLLAPSE124 122 DO jj = 1, jpj 125 123 DO ji = 1, jpi -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
r5620 r5901 4 4 !! Wave module 5 5 !!====================================================================== 6 !! History : 3.3 .1! 2011-09 (Adani M) Original code: Drag Coefficient7 !! : 3.4 6 !! History : 3.3 ! 2011-09 (Adani M) Original code: Drag Coefficient 7 !! : 3.4 ! 2012-10 (Adani M) Stokes Drift 8 8 !!---------------------------------------------------------------------- 9 USE iom ! I/O manager library 10 USE in_out_manager ! I/O manager 11 USE lib_mpp ! distribued memory computing library 9 10 !!---------------------------------------------------------------------- 11 !! sbc_wave : read drag coefficient from wave model in netcdf files 12 !!---------------------------------------------------------------------- 13 USE oce ! 14 USE sbc_oce ! Surface boundary condition: ocean fields 15 USE bdy_oce ! 16 USE domvvl ! 17 ! 18 USE iom ! I/O manager library 19 USE in_out_manager ! I/O manager 20 USE lib_mpp ! distribued memory computing library 12 21 USE fldread ! read input fields 13 USE oce 14 USE sbc_oce ! Surface boundary condition: ocean fields 15 USE domvvl 16 17 18 !!---------------------------------------------------------------------- 19 !! sbc_wave : read drag coefficient from wave model in netcdf files 20 !!---------------------------------------------------------------------- 22 USE wrk_nemo ! 21 23 22 24 IMPLICIT NONE … … 25 27 PUBLIC sbc_wave ! routine called in sbc_blk_core or sbc_blk_mfs 26 28 27 INTEGER , PARAMETER :: jpfld = 3 ! maximum number of files to read for srokes drift 28 INTEGER , PARAMETER :: jp_usd = 1 ! index of stokes drift (i-component) (m/s) at T-point 29 INTEGER , PARAMETER :: jp_vsd = 2 ! index of stokes drift (j-component) (m/s) at T-point 30 INTEGER , PARAMETER :: jp_wn = 3 ! index of wave number (1/m) at T-point 29 INTEGER , PARAMETER :: jpfld = 3 ! maximum number of files to read for srokes drift 30 INTEGER , PARAMETER :: jp_usd = 1 ! index of stokes drift (i-component) (m/s) at T-point 31 INTEGER , PARAMETER :: jp_vsd = 2 ! index of stokes drift (j-component) (m/s) at T-point 32 INTEGER , PARAMETER :: jp_wn = 3 ! index of wave number (1/m) at T-point 33 31 34 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_cd ! structure of input fields (file informations, fields read) Drag Coefficient 32 35 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_sd ! structure of input fields (file informations, fields read) Stokes Drift 33 REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:) :: cdn_wave 34 REAL(wp),ALLOCATABLE,DIMENSION (:,:) :: usd2d,vsd2d,uwavenum,vwavenum 35 REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:) :: usd3d,vsd3d,wsd3d 36 37 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION (:,:) :: cdn_wave 38 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION (:,:,:) :: usd3d, vsd3d, wsd3d 39 REAL(wp), ALLOCATABLE, DIMENSION (:,:) :: usd2d, vsd2d, uwavenum, vwavenum 36 40 37 41 !! * Substitutions 38 42 # include "domzgr_substitute.h90" 43 # include "vectopt_loop_substitute.h90" 39 44 !!---------------------------------------------------------------------- 40 !! NEMO/OPA 4.0 , NEMO Consortium (2011)45 !! NEMO/OPA 3.7 , NEMO Consortium (2014) 41 46 !! $Id$ 42 47 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 56 61 !! - Compute 3d stokes drift using monochromatic 57 62 !! ** action : 58 !!59 63 !!--------------------------------------------------------------------- 60 USE oce, ONLY : un,vn,hdivn,rotn 61 USE divcur 62 USE wrk_nemo 63 #if defined key_bdy 64 USE bdy_oce, ONLY : bdytmask 65 #endif 66 INTEGER, INTENT( in ) :: kt ! ocean time step 67 INTEGER :: ierror ! return error code 68 INTEGER :: ifpr, jj,ji,jk 64 INTEGER, INTENT( in ) :: kt ! ocean time step 65 ! 66 INTEGER :: ierror ! return error code 67 INTEGER :: ifpr, jj,ji,jk 69 68 INTEGER :: ios ! Local integer output status for namelist read 70 REAL(wp),DIMENSION(:,:,:),POINTER :: udummy,vdummy,hdivdummy,rotdummy71 REAL :: z2dt,z1_2dt72 69 TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read 73 70 CHARACTER(len=100) :: cn_dir ! Root directory for location of drag coefficient files 74 71 TYPE(FLD_N) :: sn_cdg, sn_usd, sn_vsd, sn_wn ! informations about the fields to be read 75 !!--------------------------------------------------------------------- 72 REAL(wp), DIMENSION(:,:,:), POINTER :: zusd_t, zvsd_t, ze3hdiv ! 3D workspace 73 !! 76 74 NAMELIST/namsbc_wave/ sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn 77 75 !!--------------------------------------------------------------------- 78 79 !!----------------------------------------------------------------------80 !81 76 ! 82 77 ! ! -------------------- ! … … 92 87 IF(lwm) WRITE ( numond, namsbc_wave ) 93 88 ! 94 95 89 IF ( ln_cdgw ) THEN 96 90 ALLOCATE( sf_cd(1), STAT=ierror ) !* allocate and fill sf_wave with sn_cdg … … 102 96 ALLOCATE( cdn_wave(jpi,jpj) ) 103 97 cdn_wave(:,:) = 0.0 104 ENDIF98 ENDIF 105 99 IF ( ln_sdw ) THEN 106 100 slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd; slf_i(jp_wn) = sn_wn … … 113 107 END DO 114 108 CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 115 ALLOCATE( usd2d(jpi,jpj) ,vsd2d(jpi,jpj),uwavenum(jpi,jpj),vwavenum(jpi,jpj) )109 ALLOCATE( usd2d(jpi,jpj) , vsd2d(jpi,jpj) , uwavenum(jpi,jpj) , vwavenum(jpi,jpj) ) 116 110 ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) ) 117 usd2d(:,:) = 0.0 ; vsd2d(:,:) = 0.0 ; uwavenum(:,:) = 0.0 ; vwavenum(:,:) = 0.0 118 usd3d(:,:,:) = 0.0 ;vsd3d(:,:,:) = 0.0 ; wsd3d(:,:,:) = 0.0 111 usd3d(:,:,:) = 0._wp ; usd2d(:,:) = 0._wp ; uwavenum(:,:) = 0._wp 112 vsd3d(:,:,:) = 0._wp ; vsd2d(:,:) = 0._wp ; vwavenum(:,:) = 0._wp 113 wsd3d(:,:,:) = 0._wp 119 114 ENDIF 120 115 ENDIF 116 ! 117 IF ( ln_cdgw ) THEN !== Neutral drag coefficient ==! 118 CALL fld_read( kt, nn_fsbc, sf_cd ) ! read from external forcing 119 cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 120 ENDIF 121 ! 122 IF ( ln_sdw ) THEN !== Computation of the 3d Stokes Drift ==! 123 ! 124 CALL fld_read( kt, nn_fsbc, sf_sd ) !* read drag coefficient from external forcing 121 125 ! 122 126 ! 123 IF ( ln_cdgw ) THEN 124 CALL fld_read( kt, nn_fsbc, sf_cd ) !* read drag coefficient from external forcing 125 cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 126 ENDIF 127 IF ( ln_sdw ) THEN 128 CALL fld_read( kt, nn_fsbc, sf_sd ) !* read drag coefficient from external forcing 129 130 ! Interpolate wavenumber, stokes drift into the grid_V and grid_V 131 !------------------------------------------------- 132 133 DO jj = 1, jpjm1 134 DO ji = 1, jpim1 135 uwavenum(ji,jj)=0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) & 136 & + sf_sd(3)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) ) 137 138 vwavenum(ji,jj)=0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) & 139 & + sf_sd(3)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) ) 140 141 usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(1)%fnow(ji,jj,1) * tmask(ji,jj,1) & 142 & + sf_sd(1)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) ) 143 144 vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(2)%fnow(ji,jj,1) * tmask(ji,jj,1) & 145 & + sf_sd(2)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) ) 127 CALL wrk_alloc( jpi,jpj,jpk, zusd_t, zvsd_t, ze3hdiv ) 128 ! !* distribute it on the vertical 129 DO jk = 1, jpkm1 130 zusd_t(:,:,jk) = sf_sd(jp_usd)%fnow(:,:,1) * EXP( -2._wp * sf_sd(jp_wn)%fnow(:,:,1) * fsdept_n(:,:,jk) ) 131 zvsd_t(:,:,jk) = sf_sd(jp_vsd)%fnow(:,:,1) * EXP( -2._wp * sf_sd(jp_wn)%fnow(:,:,1) * fsdept_n(:,:,jk) ) 132 END DO 133 ! !* interpolate the stokes drift from t-point to u- and v-points 134 DO jk = 1, jpkm1 135 DO jj = 1, jpjm1 136 DO ji = 1, jpim1 137 usd3d(ji,jj,jk) = 0.5_wp * ( zusd_t(ji ,jj,jk) + zusd_t(ji+1,jj,jk) ) * umask(ji,jj,jk) 138 vsd3d(ji,jj,jk) = 0.5_wp * ( zvsd_t(ji ,jj,jk) + zvsd_t(ji,jj+1,jk) ) * vmask(ji,jj,jk) 139 END DO 146 140 END DO 147 141 END DO 148 149 !Computation of the 3d Stokes Drift 150 DO jk = 1, jpk 151 DO jj = 1, jpj-1 152 DO ji = 1, jpi-1 153 usd3d(ji,jj,jk) = usd2d(ji,jj)*exp(2.0*uwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj ,jk)))) 154 vsd3d(ji,jj,jk) = vsd2d(ji,jj)*exp(2.0*vwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji ,jj+1,jk)))) 155 END DO 156 END DO 157 usd3d(jpi,:,jk) = usd2d(jpi,:)*exp( 2.0*uwavenum(jpi,:)*(-gdept_0(jpi,:,jk)) ) 158 vsd3d(:,jpj,jk) = vsd2d(:,jpj)*exp( 2.0*vwavenum(:,jpj)*(-gdept_0(:,jpj,jk)) ) 159 END DO 160 161 CALL wrk_alloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy) 162 163 udummy(:,:,:)=un(:,:,:) 164 vdummy(:,:,:)=vn(:,:,:) 165 hdivdummy(:,:,:)=hdivn(:,:,:) 166 rotdummy(:,:,:)=rotn(:,:,:) 167 un(:,:,:)=usd3d(:,:,:) 168 vn(:,:,:)=vsd3d(:,:,:) 169 CALL div_cur(kt) 170 ! !------------------------------! 171 ! ! Now Vertical Velocity ! 172 ! !------------------------------! 173 z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog) 174 175 z1_2dt = 1.e0 / z2dt 176 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 177 ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 178 wsd3d(:,:,jk) = wsd3d(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & 179 & - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) & 180 & * tmask(:,:,jk) * z1_2dt 142 CALL lbc_lnk( usd3d(:,:,:), 'U', -1. ) 143 CALL lbc_lnk( vsd3d(:,:,:), 'V', -1. ) 144 ! 145 DO jk = 1, jpkm1 !* e3t * Horizontal divergence ==! 146 DO jj = 2, jpjm1 147 DO ji = fs_2, fs_jpim1 ! vector opt. 148 ze3hdiv(ji,jj,jk) = ( e2u(ji ,jj) * fse3u_n(ji ,jj,jk) * usd3d(ji ,jj,jk) & 149 & - e2u(ji-1,jj) * fse3u_n(ji-1,jj,jk) * usd3d(ji-1,jj,jk) & 150 & + e1v(ji,jj ) * fse3v_n(ji,jj ,jk) * vsd3d(ji,jj ,jk) & 151 & - e1v(ji,jj-1) * fse3v_n(ji,jj-1,jk) * vsd3d(ji,jj-1,jk) ) * r1_e1e2t(ji,jj) 152 END DO 153 END DO 154 IF( .NOT. AGRIF_Root() ) THEN 155 IF( nbondi == 1 .OR. nbondi == 2 ) ze3hdiv(nlci-1, : ,jk) = 0._wp ! east 156 IF( nbondi == -1 .OR. nbondi == 2 ) ze3hdiv( 2 , : ,jk) = 0._wp ! west 157 IF( nbondj == 1 .OR. nbondj == 2 ) ze3hdiv( : ,nlcj-1,jk) = 0._wp ! north 158 IF( nbondj == -1 .OR. nbondj == 2 ) ze3hdiv( : , 2 ,jk) = 0._wp ! south 159 ENDIF 160 END DO 161 CALL lbc_lnk( ze3hdiv, 'T', 1. ) 162 ! 163 DO jk = jpkm1, 1, -1 !* integrate from the bottom the e3t * hor. divergence 164 wsd3d(:,:,jk) = wsd3d(:,:,jk+1) - ze3hdiv(:,:,jk) 165 END DO 181 166 #if defined key_bdy 182 wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:) 167 IF( lk_bdy ) THEN 168 DO jk = 1, jpkm1 169 wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:) 170 END DO 171 ENDIF 183 172 #endif 184 END DO 185 hdivn(:,:,:)=hdivdummy(:,:,:) 186 rotn(:,:,:)=rotdummy(:,:,:) 187 vn(:,:,:)=vdummy(:,:,:) 188 un(:,:,:)=udummy(:,:,:) 189 CALL wrk_dealloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy) 173 CALL wrk_dealloc( jpi,jpj,jpk, zusd_t, zvsd_t, ze3hdiv ) 174 ! 190 175 ENDIF 176 ! 191 177 END SUBROUTINE sbc_wave 192 178
Note: See TracChangeset
for help on using the changeset viewer.