Changeset 2799


Ignore:
Timestamp:
2011-07-12T19:12:16+02:00 (9 years ago)
Author:
vichi
Message:

Updated sbcblk_core.F90 so it can be used for both ORCA2 and ORCA025

Location:
branches/2011/dev_r2784_CMCC1_topbfm/NEMOGCM
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/2011/dev_r2784_CMCC1_topbfm/NEMOGCM/ARCH/arch-PW6_CALYPSO.fcm

    r2798 r2799  
    1717%NCDF_INC            -I/usr/local/netcdf-3.6.3/include 
    1818%NCDF_LIB            -L /usr/local/netcdf-3.6.3/lib -lnetcdf 
    19 %FC                  xlf90_r 
     19%FC                  mpxlf90_r 
    2020%FCFLAGS             -qfree=f90 -O3 -qrealsize=8 -qextname -qsource -q64 -qlargepage -qmaxmem=-1 
    2121%FFLAGS              -qfixed -O3 -qrealsize=8 -qextname -qsource -q64 -qlargepage -qmaxmem=-1 
    22 %LD                  xlf90_r 
     22%LD                  mpxlf90_r 
    2323%LDFLAGS 
    2424%FPPFLAGS            -P -C 
  • branches/2011/dev_r2784_CMCC1_topbfm/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r2777 r2799  
    4444   PUBLIC   blk_ice_core         ! routine called in sbc_ice_lim module 
    4545 
    46    INTEGER , PARAMETER ::   jpfld   = 9           ! maximum number of files to read  
    4746   INTEGER , PARAMETER ::   jp_wndi = 1           ! index of 10m wind velocity (i-component) (m/s)    at T-point 
    4847   INTEGER , PARAMETER ::   jp_wndj = 2           ! index of 10m wind velocity (j-component) (m/s)    at T-point 
     
    5453   INTEGER , PARAMETER ::   jp_snow = 8           ! index of snow (solid prcipitation)       (kg/m2/s) 
    5554   INTEGER , PARAMETER ::   jp_tdif = 9           ! index of tau diff associated to HF tau   (N/m2)   at T-point 
     55#if defined key_orca_r025 
     56   INTEGER , PARAMETER ::   jp_swc  = 10          ! index of GEWEX correction for SW radiation  at T-point 
     57   INTEGER , PARAMETER ::   jp_lwc  = 11          ! index of GEWEX correction for LW radiation  at T-point 
     58   INTEGER , PARAMETER ::   jp_prc  = 12          ! index of PMWC correction forat T-point 
     59   INTEGER , PARAMETER ::   jpfld   = 12          ! maximum number of files to read 
     60#else 
     61   INTEGER , PARAMETER ::   jpfld   = 9           ! maximum number of files to read 
     62#endif 
    5663    
    5764   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read) 
     
    7178   REAL(wp) ::   rn_pfac   = 1.        ! multiplication factor for precipitation 
    7279 
     80#if defined key_orca_r025 
     81   LOGICAL  ::   ln_printdia= .TRUE.     ! logical flag for height of air temp. and hum 
     82   LOGICAL  ::   ln_netsw   = .TRUE.     ! logical flag for height of air temp. and hum 
     83   LOGICAL  ::   ln_core_graceopt=.FALSE., ln_core_spinup=.FALSE. 
     84   LOGICAL  ::   ln_gwxc = .TRUE. 
     85   LOGICAL  ::   ln_corad_antar =.FALSE., ln_corad_arc =.FALSE. , ln_cotair_arc = .FALSE. 
     86   LOGICAL  ::   ln_coprecip =.FALSE. 
     87   REAL(wp) ::   rn_qns_bias = 0._wp     ! heat flux bias 
     88 
     89   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  area  
     90   REAL(wp) :: araux 
     91 
     92   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  zqlw, zqsb  ! long wave and sensible heat fluxes 
     93   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  zqla, zevap ! latent heat fluxes and evaporation 
     94 
     95   REAL(wp), PARAMETER :: zalph = 2.408724e-06_wp, & 
     96   & zbet = -0.006936579_wp, zgam = 449.9094_wp ! GRACE regression coefficients 
     97#endif 
     98 
    7399   !! * Substitutions 
    74100#  include "domzgr_substitute.h90" 
     
    80106   !!---------------------------------------------------------------------- 
    81107CONTAINS 
    82  
     108    
    83109   SUBROUTINE sbc_blk_core( kt ) 
    84110      !!--------------------------------------------------------------------- 
     
    112138      !!              - emp, emps   evaporation minus precipitation 
    113139      !!---------------------------------------------------------------------- 
     140#if defined key_orca_r025 
     141      USE ice_2 
     142#endif 
    114143      INTEGER, INTENT(in) ::   kt   ! ocean time step 
    115144      !! 
     
    117146      INTEGER  ::   ifpr     ! dummy loop indice 
    118147      INTEGER  ::   jfld     ! dummy loop arguments 
     148      INTEGER  ::   ji, jj 
    119149      !! 
    120150      CHARACTER(len=100) ::  cn_dir   !   Root directory for location of core files 
     
    123153      TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !   "                                 " 
    124154      TYPE(FLD_N) ::   sn_tdif                                 !   "                                 " 
     155      TYPE(FLD_N) ::   sn_swc, sn_lwc                          !   "                                 " 
     156      TYPE(FLD_N) ::   sn_prc 
     157      INTEGER  ::   iter_shapiro = 250 
     158      REAL :: zzlat, zzlat1, zzlat2, zfm, zfrld, ztmp 
     159      REAL(wp), DIMENSION(jpi,jpj):: xyt,z_qsr,z_qlw,z_qsr1,z_qlw1, z_hum, z_tair 
     160      REAL(wp), DIMENSION(jpi,jpj):: zqsr_lr, zqsr_hr, zqlw_lr, zqlw_hr, zprec_hr, zprec_lr 
     161      CHARACTER(len=20)  ::  c_kind='ORCA_GLOB' 
     162#if defined key_orca_r025 
     163      NAMELIST/namsbc_core/ cn_dir , ln_2m  , ln_taudif, rn_pfac,           & 
     164         &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           & 
     165         &                  sn_qlw , sn_tair, sn_prec  , sn_snow, sn_tdif   & 
     166         &                  sn_swc , sn_lwc , sn_prc   , ln_gwxc  ,         & 
     167         &                  ln_corad_antar, ln_corad_arc, ln_cotair_arc, ln_coprecip ,  & 
     168         &                  rn_qns_bias, ln_printdia, ln_netsw, ln_core_graceopt,ln_core_spinup 
     169      !!--------------------------------------------------------------------- 
     170#else 
    125171      NAMELIST/namsbc_core/ cn_dir , ln_2m  , ln_taudif, rn_pfac,           & 
    126172         &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           & 
    127173         &                  sn_qlw , sn_tair, sn_prec  , sn_snow, sn_tdif 
    128174      !!--------------------------------------------------------------------- 
     175#endif 
    129176 
    130177      !                                         ! ====================== ! 
    131178      IF( kt == nit000 ) THEN                   !  First call kt=nit000  ! 
    132179         !                                      ! ====================== ! 
     180#if defined key_orca_r025 
     181         !                                      !==  allocate sbc arrays 
     182            IF( sbc_blk_core_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_core_alloc : unable to allocate arrays' ) 
     183#endif 
     184 
    133185         ! set file information (default values) 
    134186         cn_dir = './'       ! directory in which the model is executed 
     
    146198         sn_snow = FLD_N( 'snow'   ,    -1     , 'snow'   ,  .true.    , .false. ,   'yearly'  , ''       , ''       ) 
    147199         sn_tdif = FLD_N( 'taudif' ,    24     , 'taudif' ,  .true.    , .false. ,   'yearly'  , ''       , ''       ) 
     200#if defined key_orca_r025 
     201         sn_swc  = FLD_N( 'swc'    ,    24     ,  'swc'   ,  .true.    , .false. ,   'yearly'  , ''       , ''       ) 
     202         sn_lwc  = FLD_N( 'lwc'    ,    24     ,  'lwc'   ,  .true.    , .false. ,   'yearly'  , ''       , ''       ) 
     203         sn_prc  = FLD_N( 'prc'    ,    24     ,  'prc'   ,  .true.    , .false. ,   'yearly'  , ''       , ''       ) 
     204#endif 
    148205         ! 
    149206         REWIND( numnam )                          ! read in namlist namsbc_core 
     
    163220         slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow 
    164221         slf_i(jp_tdif) = sn_tdif 
     222#if defined key_orca_r025 
     223            slf_i(jp_swc ) = sn_swc 
     224            slf_i(jp_lwc ) = sn_lwc 
     225            slf_i(jp_prc ) = sn_prc 
     226#endif 
    165227         !                  
    166228         lhftau = ln_taudif                        ! do we use HF tau information? 
    167229         jfld = jpfld - COUNT( (/.NOT. lhftau/) ) 
     230#if defined key_orca_r025 
     231         IF( .NOT. ln_gwxc )     jfld = jfld - 2 
     232         IF( .NOT. ln_coprecip ) jfld = jfld - 1 
     233#endif 
    168234         ! 
    169235         ALLOCATE( sf(jfld), STAT=ierror )         ! set sf structure 
     
    176242         CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_core', 'flux formulation for ocean surface boundary condition', 'namsbc_core' ) 
    177243         ! 
     244#if defined key_orca_r025 
     245         IF( ln_printdia .OR. ln_core_graceopt ) THEN 
     246           area = (e1t * e2t) 
     247           araux = sum ( area * tmask(:,:,1) ) 
     248           IF(lk_mpp) CALL mpp_sum ( araux ) 
     249         ENDIF 
     250#endif 
    178251      ENDIF 
    179252 
    180253      CALL fld_read( kt, nn_fsbc, sf )        ! input fields provided at the current time-step 
    181254 
    182       !                                                        ! surface ocean fluxes computed with CLIO bulk formulea 
    183       IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce_core( sf, sst_m, ssu_m, ssv_m ) 
    184       ! 
     255      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 
     256 
     257#if defined key_orca_r025 
     258      ! Introduce ERA-Interim filtering and correction 
     259 
     260         IF( ln_gwxc ) THEN 
     261 
     262           call Shapiro_1D(sf(jp_qsr)%fnow(:,:,1),iter_shapiro, c_kind, zqsr_lr) 
     263           zqsr_hr(:,:)=sf(jp_qsr)%fnow(:,:,1)-zqsr_lr(:,:)          ! We get large scale and small scale 
     264 
     265           call Shapiro_1D(sf(jp_qlw)%fnow(:,:,1),iter_shapiro, c_kind, zqlw_lr) 
     266           zqlw_hr(:,:)=sf(jp_qlw)%fnow(:,:,1)-zqlw_lr(:,:)          ! We get large scale and small scale 
     267 
     268           z_qsr1(:,:)=zqsr_lr(:,:)*sf(jp_swc)%fnow(:,:,1) + zqsr_hr(:,:) 
     269           z_qlw1(:,:)=zqlw_lr(:,:)*sf(jp_lwc)%fnow(:,:,1) + zqlw_hr(:,:) 
     270 
     271           DO jj=1,jpj 
     272             DO ji=1,jpi 
     273               z_qsr1(ji,jj)=max(z_qsr1(ji,jj),0.0) 
     274               z_qlw1(ji,jj)=max(z_qlw1(ji,jj),0.0) 
     275             END DO 
     276           END DO 
     277 
     278         ENDIF 
     279 
     280         IF( ln_coprecip ) THEN 
     281 
     282           call Shapiro_1D(sf(jp_prec)%fnow(:,:,1),iter_shapiro,c_kind,zprec_lr) 
     283           zprec_hr(:,:)=sf(jp_prec)%fnow(:,:,1)-zprec_lr(:,:)       ! We get large scale and small scale 
     284 
     285           DO jj=1,jpj 
     286             DO ji=1,jpi 
     287               IF( zprec_lr(ji,jj) .GT. 0._wp ) THEN 
     288                  ztmp = LOG( ( 1000._wp + sf(jp_prc)%fnow(ji,jj,1) ) * EXP( zprec_lr(ji,jj) ) / 1000._wp ) 
     289                  sf(jp_prec)%fnow(ji,jj,1) = max(ztmp+zprec_hr(ji,jj),0.0) 
     290               ENDIF 
     291             END DO 
     292           END DO 
     293 
     294         ENDIF 
     295 
     296         IF ( ln_corad_antar ) THEN           ! correction of SW and LW in the Southern Ocean 
     297  
     298           z_qsr(:,:)=0.8*z_qsr1(:,:) 
     299           z_qlw(:,:)=1.1*z_qlw1(:,:) 
     300           xyt(:,:) = 0.e0 
     301           zzlat1 = -65. 
     302           zzlat2 = -60. 
     303           DO jj = 1, jpj 
     304             DO ji = 1, jpi 
     305               zzlat = gphit(ji,jj) 
     306               IF ( zzlat >= zzlat1 .AND. zzlat <= zzlat2 ) THEN 
     307                  xyt(ji,jj) = (zzlat2-zzlat)/(zzlat2-zzlat1) 
     308               ELSE IF ( zzlat < zzlat1 ) THEN 
     309                  xyt(ji,jj) = 1 
     310               ENDIF 
     311             END DO 
     312           END DO 
     313           z_qsr1(:,:)=z_qsr(:,:)*xyt(:,:)+(1.0-xyt(:,:))*z_qsr1(:,:) 
     314           z_qlw1(:,:)=z_qlw(:,:)*xyt(:,:)+(1.0-xyt(:,:))*z_qlw1(:,:) 
     315 
     316         ENDIF 
     317 
     318         IF ( ln_corad_arc ) THEN         ! correction of SW in the Arctic Ocean 
     319 
     320           z_qsr(:,:)=0.7*z_qsr1(:,:) 
     321           xyt(:,:) = 0.e0 
     322           zzlat1 = 78. 
     323           zzlat2 = 82. 
     324           DO jj = 1, jpj 
     325             DO ji = 1, jpi 
     326               zzlat = gphit(ji,jj) 
     327               IF ( zzlat >= zzlat1 .AND. zzlat <= zzlat2 ) THEN 
     328                  xyt(ji,jj) = (zzlat-zzlat1)/(zzlat2-zzlat1) 
     329               ELSE IF ( zzlat > zzlat2 ) THEN 
     330                  xyt(ji,jj) = 1 
     331               ENDIF 
     332             END DO 
     333           END DO 
     334           z_qsr1(:,:)=z_qsr(:,:)*xyt(:,:)+(1.0-xyt(:,:))*z_qsr1(:,:) 
     335 
     336         ENDIF 
     337 
     338         sf(jp_qsr)%fnow(:,:,1)=z_qsr1(:,:) 
     339         sf(jp_qlw)%fnow(:,:,1)=z_qlw1(:,:) 
     340 
     341         IF ( ln_cotair_arc ) THEN           ! correction of Air Temperature in the Arctic Ocean 
     342 
     343           z_tair(:,:)=sf(jp_tair)%fnow(:,:,1) - 2.0 
     344           xyt(:,:) = 0.e0 ; zzlat1 = 78. ; zzlat2 = 82. 
     345           DO jj = 1, jpj 
     346             DO ji = 1, jpi 
     347               zzlat = gphit(ji,jj) ; zfrld=frld(ji,jj) 
     348               IF ( zzlat >= zzlat1 .AND. zzlat <= zzlat2 .AND. zfrld < 0.85 ) THEN 
     349                  xyt(ji,jj) = (zzlat-zzlat1)/(zzlat2-zzlat1) 
     350               ELSE IF ( zzlat > zzlat2 .AND. zfrld < 0.85 ) THEN 
     351                  xyt(ji,jj) = 1 
     352               ENDIF 
     353             END DO 
     354           END DO 
     355           sf(jp_tair)%fnow(:,:,1)=z_tair(:,:)*xyt(:,:)+(1.0-xyt(:,:))*sf(jp_tair)%fnow(:,:,1) 
     356 
     357         ENDIF 
     358 
     359#endif 
     360 
     361         CALL blk_oce_core( sf, sst_m, ssu_m, ssv_m )   ! compute the surface ocean fluxes using CLIO bulk formulea 
     362 
     363      ENDIF 
     364      !                                                  ! using CORE bulk formulea 
    185365   END SUBROUTINE sbc_blk_core 
    186366    
     
    315495      IF( lhftau ) THEN  
    316496!CDIR COLLAPSE 
     497#if defined key_orca_r025 
     498         ! Changed!!! Multiply by QSCAT correction 
     499           zwnd_i(:,:) = zwnd_i(:,:) * sf(jp_tdif)%fnow(:,:,1) 
     500           zwnd_j(:,:) = zwnd_j(:,:) * sf(jp_tdif)%fnow(:,:,1) 
     501#endif 
    317502         taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1) 
    318503      ENDIF 
     
    9381123      ! 
    9391124    END FUNCTION psi_h 
     1125 
     1126#if defined key_orca_r025 
     1127    INTEGER FUNCTION sbc_blk_core_alloc() 
     1128      !!---------------------------------------------------------------------- 
     1129      !!                ***  ROUTINE sbc_blk_core_alloc  *** 
     1130      !!---------------------------------------------------------------------- 
     1131      ALLOCATE( area(jpi,jpj)         , zqlw(jpi,jpj)      ,     & 
     1132         &      zqsb(jpi,jpj)         , zqla(jpi,jpj)      ,     & 
     1133         &      zevap(jpi,jpj)        , STAT=sbc_blk_core_alloc ) 
     1134      zqlw=0._wp 
     1135      zqsb=0._wp 
     1136   ! 
     1137      IF( lk_mpp            )   CALL mpp_sum ( sbc_blk_core_alloc ) 
     1138      IF( sbc_blk_core_alloc > 0 )   CALL ctl_warn('sbc_blk_core_alloc: allocation of arrays failed') 
     1139    END FUNCTION sbc_blk_core_alloc 
     1140#endif 
     1141 
     1142    SUBROUTINE Shapiro_1D(rla_varin,id_np, cd_overlap, rlpa_varout) !GIG 
     1143      !!===================================================================== 
     1144      !! 
     1145      !! Description: This function applies a 1D Shapiro filter 
     1146      !!              (3 points filter) horizontally to a 2D field 
     1147      !!              in regular grid 
     1148      !! Arguments : 
     1149      !!            rla_varin   : Input variable to filter 
     1150      !!            zla_mask    : Input mask variable 
     1151      !!            id_np       : Number of Shapiro filter iterations 
     1152      !!            cd_overlap  : Logical argument for periodical condition 
     1153      !!                          (global ocean case) 
     1154      !!            rlpa_varout : Output filtered variable 
     1155      !! 
     1156      !! History : 08/2009  S. CAILLEAU : from 1st version of N. FERRY 
     1157      !!           09/2009  C. REGNIER  : Corrections 
     1158      !! 
     1159      !!===================================================================== 
     1160      IMPLICIT NONE 
     1161      INTEGER, INTENT(IN)                       :: id_np 
     1162      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN)  :: rla_varin !GIG 
     1163      CHARACTER(len=20), INTENT(IN)             :: cd_overlap !GIG 
     1164      REAL(wp), DIMENSION(jpi,jpj), INTENT(OUT) :: rlpa_varout !GIG 
     1165 
     1166      REAL(wp), DIMENSION(jpi,jpj)              :: rlpa_varout_tmp 
     1167      REAL, PARAMETER                           :: rl_alpha = 1./2.    ! fixed stability coefficient (isotrope case) 
     1168      REAL, parameter                           :: rap_aniso_diff_XY=2.25 ! anisotrope case 
     1169      REAL                                      :: alphax,alphay, znum, zden,test 
     1170      INTEGER                                   :: ji, jj, jn, nn 
     1171! 
     1172!! rap_aniso_diff_XY=2.25 : valeur trouvée empiriquement pour 140 itération pour le filtre de Shapiro et 
     1173!! pour un rapport d'anisotopie de 1.5 : on filtre de plus rapidement en x qu'eny. 
     1174!------------------------------------------------------------------------------ 
     1175! 
     1176! Loop on several filter iterations 
     1177 
     1178!     Global ocean case 
     1179      IF (( cd_overlap == 'MERCA_GLOB' )   .OR.   & 
     1180          ( cd_overlap == 'REGULAR_GLOB' ) .OR.   & 
     1181          ( cd_overlap == 'ORCA_GLOB' )) THEN 
     1182             rlpa_varout(:,:) = rla_varin(:,:) 
     1183             rlpa_varout_tmp(:,:) = rlpa_varout(:,:) 
     1184! 
     1185 
     1186       alphax=1./2. 
     1187       alphay=1./2. 
     1188!  Dx/Dy=rap_aniso_diff_XY  , D_ = vitesse de diffusion 
     1189!  140 passes du fitre, Lx/Ly=1.5, le rap_aniso_diff_XY correspondant est: 
     1190       IF ( rap_aniso_diff_XY .GE. 1. ) alphay=alphay/rap_aniso_diff_XY 
     1191       IF ( rap_aniso_diff_XY .LT. 1. ) alphax=alphax*rap_aniso_diff_XY 
     1192 
     1193        DO jn = 1,id_np   ! number of passes of the filter 
     1194            DO ji = 2,jpim1 
     1195               DO jj = 2,jpjm1 
     1196                  ! We crop on the coast 
     1197                   znum = rlpa_varout_tmp(ji,jj)   & 
     1198                          + 0.25*alphax*(rlpa_varout_tmp(ji-1,jj  )-rlpa_varout_tmp(ji,jj))*tmask(ji-1,jj  ,1)  & 
     1199                          + 0.25*alphax*(rlpa_varout_tmp(ji+1,jj  )-rlpa_varout_tmp(ji,jj))*tmask(ji+1,jj  ,1)  & 
     1200                          + 0.25*alphay*(rlpa_varout_tmp(ji  ,jj-1)-rlpa_varout_tmp(ji,jj))*tmask(ji  ,jj-1,1)  & 
     1201                          + 0.25*alphay*(rlpa_varout_tmp(ji  ,jj+1)-rlpa_varout_tmp(ji,jj))*tmask(ji  ,jj+1,1) 
     1202                   rlpa_varout(ji,jj)=znum*tmask(ji,jj,1)+rla_varin(ji,jj)*(1.-tmask(ji,jj,1)) 
     1203                ENDDO  ! end loop ji 
     1204            ENDDO  ! end loop jj 
     1205! 
     1206! 
     1207!           Periodical condition in case of cd_overlap (global ocean) 
     1208!           - on a mercator projection grid we consider that singular point at poles 
     1209!             are a mean of the values at points of the previous latitude 
     1210!           - on ORCA and regular grid we copy the values at points of the previous latitude 
     1211            IF ( cd_overlap == 'MERCAT_GLOB' ) THEN 
     1212!GIG case unchecked 
     1213               rlpa_varout(1,1) = SUM(rlpa_varout(:,2)) / jpi 
     1214               rlpa_varout(jpi,jpj) = SUM(rlpa_varout(:,jpj-1)) / jpi 
     1215            ELSE 
     1216               call lbc_lnk(rlpa_varout, 'T', 1.) ! Boundary condition 
     1217            ENDIF 
     1218            rlpa_varout_tmp(:,:) = rlpa_varout(:,:) 
     1219         ENDDO  ! end loop jn 
     1220      ENDIF 
     1221 
     1222! 
     1223    END SUBROUTINE Shapiro_1D 
     1224 
    9401225   
    9411226   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.