Changeset 7168


Ignore:
Timestamp:
2016-11-02T15:24:08+01:00 (4 years ago)
Author:
jcastill
Message:

First trial to changes needed for wave coupling

Location:
branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r7166 r7168  
    2323   USE sbcapr 
    2424   USE sbcdcy          ! surface boundary condition: diurnal cycle 
     25   USE sbcwave         ! surface boundary condition: waves 
    2526   USE phycst          ! physical constants 
    2627#if defined key_lim3 
     
    105106   INTEGER, PARAMETER ::   jpr_e3t1st = 41            ! first T level thickness  
    106107   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   
    108118 
    109119   INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction sent to the atmosphere 
     
    135145   INTEGER, PARAMETER ::   jps_e3t1st = 27            ! first level depth (vvl) 
    136146   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  
    138152 
    139153   !                                                         !!** namelist namsbc_cpl ** 
     
    149163   ! Received from the atmosphere                     ! 
    150164   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 
    152170   ! Other namelist parameters                        ! 
    153171   INTEGER     ::   nn_cplmodel            ! Maximum number of models to/from which NEMO is potentialy sending/receiving data 
     
    161179 
    162180   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)  
    163184 
    164185   INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo           ! OASIS info argument 
     
    179200      !!             ***  FUNCTION sbc_cpl_alloc  *** 
    180201      !!---------------------------------------------------------------------- 
    181       INTEGER :: ierr(3) 
     202      INTEGER :: ierr(4) 
    182203      !!---------------------------------------------------------------------- 
    183204      ierr(:) = 0 
     
    190211      ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) ) 
    191212      ! 
     213      IF( .NOT. ln_apr_dyn ) ALLOCATE( ssh_ib(jpi,jpj), ssh_ibb(jpi,jpj), apr(jpi, jpj), STAT=ierr(4) )  
     214 
    192215      sbc_cpl_alloc = MAXVAL( ierr ) 
    193216      IF( lk_mpp            )   CALL mpp_sum ( sbc_cpl_alloc ) 
     
    216239      REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos 
    217240      !! 
    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  
    222247      !!--------------------------------------------------------------------- 
    223248      ! 
     
    260285         WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')' 
    261286         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 ), ')'  
    262295         WRITE(numout,*)'  sent fields (multiple ice categories)' 
    263296         WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')' 
    264297         WRITE(numout,*)'      albedo                          = ', TRIM(sn_snd_alb%cldes   ), ' (', TRIM(sn_snd_alb%clcat   ), ')' 
    265298         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 ), ')'  
    266300         WRITE(numout,*)'      surface current                 = ', TRIM(sn_snd_crt%cldes   ), ' (', TRIM(sn_snd_crt%clcat   ), ')' 
    267301         WRITE(numout,*)'                      - referential   = ', sn_snd_crt%clvref  
     
    269303         WRITE(numout,*)'                      - mesh          = ', sn_snd_crt%clvgrd 
    270304         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  
    271311         WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel 
    272312         WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask 
     
    470510      !                                                      ! ------------------------- ! 
    471511      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 
    472518      !                                                      ! ------------------------- ! 
    473519      !                                                      !   topmelt and botmelt     !    
     
    483529         srcv(jpr_topm:jpr_botm)%laction = .TRUE. 
    484530      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      !  
    485575      !                                                      ! ------------------------------- ! 
    486576      !                                                      !   OPA-SAS coupling - rcv by opa !    
     
    637727      !                                                      ! ------------------------- ! 
    638728      ssnd(jps_fice)%clname = 'OIceFrc' 
     729      ssnd(jps_ficet)%clname = 'OIceFrcT'  
    639730      ssnd(jps_hice)%clname = 'OIceTck' 
    640731      ssnd(jps_hsnw)%clname = 'OSnwTck' 
     
    645736      ENDIF 
    646737       
     738      IF (TRIM( sn_snd_ifrac%cldes )  == 'coupled') ssnd(jps_ficet)%laction = .TRUE.  
     739 
    647740      SELECT CASE ( TRIM( sn_snd_thick%cldes ) ) 
    648741      CASE( 'none'         )       ! nothing to do 
     
    665758      ssnd(jps_ocy1)%clname = 'O_OCury1'   ;   ssnd(jps_ivy1)%clname = 'O_IVely1' 
    666759      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'  
    667762      ! 
    668763      ssnd(jps_ocx1:jps_ivz1)%nsgn = -1.   ! vectors: change of the sign at the north fold 
     
    685780      END SELECT 
    686781 
     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 
    687798      !                                                      ! ------------------------- ! 
    688799      !                                                      !          CO2 flux         ! 
    689800      !                                                      ! ------------------------- ! 
    690801      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.  
    691807 
    692808      !                                                      ! ------------------------------- ! 
     
    837953      !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case) 
    838954      !!---------------------------------------------------------------------- 
     955      USE zdf_oce,  ONLY : ln_zdfqiao 
     956 
     957      IMPLICIT NONE 
     958 
    839959      INTEGER, INTENT(in)           ::   kt          ! ocean model time step index 
    840960      INTEGER, INTENT(in)           ::   k_fsbc      ! frequency of sbc (-> ice model) computation  
     
    9921112      IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1) 
    9931113#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) 
    9941172 
    9951173      !  Fields received by SAS when OASIS coupling 
     
    19842162      ENDIF 
    19852163      ! 
     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  
    19862308      ! 
    19872309      !  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  
    2727   PRIVATE 
    2828 
     29   PUBLIC   sbc_stokes, sbc_qiao  ! routines called in sbccpl 
    2930   PUBLIC   sbc_wave    ! routine called in sbcmod 
    3031    
    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 
    3647 
    3748   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient 
     
    5768CONTAINS 
    5869 
     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 
    59179   SUBROUTINE sbc_wave( kt ) 
    60180      !!--------------------------------------------------------------------- 
     
    71191      !! ** action   
    72192      !!--------------------------------------------------------------------- 
    73       USE zdf_oce,  ONLY : ln_zdfqiao 
     193      USE zdf_oce    ,  ONLY : ln_zdfqiao 
     194 
     195      IMPLICIT NONE 
    74196 
    75197      INTEGER, INTENT( in  ) ::   kt       ! ocean time step 
    76198      ! 
    77199      INTEGER                ::   ierror   ! return error code 
    78       INTEGER                ::   ifpr, jj,ji,jk  
     200      INTEGER                ::   ifpr 
    79201      INTEGER                ::   ios      ! Local integer output status for namelist read 
    80202      ! 
    81203      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 
    86205      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  & 
    87206                             &   sn_swh, sn_wmp, sn_wnum, sn_tauoc      ! informations about the fields to be read 
     
    103222         ! 
    104223         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 
    111232            ALLOCATE( cdn_wave(jpi,jpj) ) 
    112233            cdn_wave(:,:) = 0.0 
     
    114235 
    115236         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 
    122245            ALLOCATE( tauoc_wave(jpi,jpj) ) 
    123246            tauoc_wave(:,:) = 0.0 
    124         ENDIF 
     247         ENDIF 
    125248 
    126249         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 
    138287            ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj) ) 
    139288            ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) ) 
     
    145294            swh  (:,:)   = 0._wp   ;    wmp (:,:) = 0._wp ; 
    146295            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 
    152303               ALLOCATE( wnum(jpi,jpj),tsd2d(jpi,jpj) ) 
    153304               wnum(:,:) = 0._wp ; tsd2d(:,:) = 0._wp 
     
    156307      ENDIF 
    157308      ! 
    158       IF ( ln_cdgw ) THEN              !==  Neutral drag coefficient  ==! 
     309      IF ( ln_cdgw .AND. .NOT. cpl_wdrag ) THEN              !==  Neutral drag coefficient  ==! 
    159310         CALL fld_read( kt, nn_fsbc, sf_cd )      ! read from external forcing 
    160311         cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 
    161312      ENDIF 
    162313 
    163       IF ( ln_tauoc ) THEN             !==  Wave induced stress  ==! 
     314      IF ( ln_tauoc .AND. .NOT. cpl_wstrf ) THEN             !==  Wave induced stress  ==! 
    164315         CALL fld_read( kt, nn_fsbc, sf_tauoc )      !* read wave norm stress from external forcing 
    165316         tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1) 
     
    168319      IF ( ln_sdw )  THEN                         !==  Computation of the 3d Stokes Drift  ==!  
    169320         ! 
    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 
    175329         ! 
     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            
    176336         !==  Computation of the 3d Stokes Drift according to Breivik et al.,2014 
    177337         !(DOI: 10.1175/JPO-D-14-0020.1)==!  
    178338         ! 
    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 
    262346      ENDIF 
    263347      ! 
  • branches/UKMO/r5936_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/step.F90

    r7167 r7168  
    231231                         CALL dyn_adv       ( kstp )  ! advection (vector or flux form) 
    232232                         CALL dyn_vor       ( kstp )  ! vorticity term including Coriolis 
     233      IF( ln_stcor    )  CALL dyn_stcor     ( kstp )  ! Stokes-Coriolis forcing 
    233234                         CALL dyn_ldf       ( kstp )  ! lateral mixing 
    234       IF( ln_stcor    )  CALL dyn_stcor     ( kstp )  ! Stokes-Coriolis forcing 
    235235                         CALL dyn_hpg       ( kstp )  ! horizontal gradient of Hydrostatic pressure 
    236236                         CALL dyn_spg       ( kstp )  ! surface pressure gradient 
Note: See TracChangeset for help on using the changeset viewer.