Changeset 6140 for trunk/NEMOGCM/NEMO/OPA_SRC/SBC
- Timestamp:
- 2015-12-21T12:35:23+01:00 (8 years ago)
- Location:
- trunk/NEMOGCM/NEMO/OPA_SRC/SBC
- Files:
-
- 23 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
r5836 r6140 4 4 !! Ocean forcing: read input field for surface boundary condition 5 5 !!===================================================================== 6 !! History : 2.0 ! 06-2006 (S. Masson, G. Madec) Original code 7 !! ! 05-2008 (S. Alderson) Modified for Interpolation in memory 8 !! ! from input grid to model grid 9 !! ! 10-2013 (D. Delrosso, P. Oddo) implement suppression of 10 !! ! land point prior to interpolation 6 !! History : 2.0 ! 06-2006 (S. Masson, G. Madec) Original code 7 !! ! 05-2008 (S. Alderson) Modified for Interpolation in memory from input grid to model grid 8 !! ! 10-2013 (D. Delrosso, P. Oddo) suppression of land point prior to interpolation 11 9 !!---------------------------------------------------------------------- 12 10 … … 15 13 !! surface boundary condition 16 14 !!---------------------------------------------------------------------- 17 USE oce ! ocean dynamics and tracers 18 USE dom_oce ! ocean space and time domain 19 USE phycst ! ??? 20 USE in_out_manager ! I/O manager 21 USE iom ! I/O manager library 22 USE geo2ocean ! for vector rotation on to model grid 23 USE lib_mpp ! MPP library 24 USE wrk_nemo ! work arrays 25 USE lbclnk ! ocean lateral boundary conditions (C1D case) 26 USE ioipsl, ONLY : ymds2ju, ju2ymds ! for calendar 27 USE sbc_oce 15 USE oce ! ocean dynamics and tracers 16 USE dom_oce ! ocean space and time domain 17 USE phycst ! physical constant 18 USE sbc_oce ! surface boundary conditions : fields 19 USE geo2ocean ! for vector rotation on to model grid 20 ! 21 USE in_out_manager ! I/O manager 22 USE iom ! I/O manager library 23 USE ioipsl , ONLY : ymds2ju, ju2ymds ! for calendar 24 USE lib_mpp ! MPP library 25 USE wrk_nemo ! work arrays 26 USE lbclnk ! ocean lateral boundary conditions (C1D case) 28 27 29 28 IMPLICIT NONE … … 60 59 INTEGER , DIMENSION(2) :: nrec_b ! before record (1: index, 2: second since Jan. 1st 00h of nit000 year) 61 60 INTEGER , DIMENSION(2) :: nrec_a ! after record (1: index, 2: second since Jan. 1st 00h of nit000 year) 62 REAL(wp) , ALLOCATABLE, DIMENSION(:,:,: ) :: fnow 63 REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:,:) :: fdta 61 REAL(wp) , ALLOCATABLE, DIMENSION(:,:,: ) :: fnow ! input fields interpolated to now time step 62 REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:,:) :: fdta ! 2 consecutive record of input fields 64 63 CHARACTER(len = 256) :: wgtname ! current name of the NetCDF weight file acting as a key 65 64 ! ! into the WGTLIST structure … … 133 132 INTEGER , INTENT(in ), OPTIONAL :: kit ! subcycle timestep for timesplitting option 134 133 INTEGER , INTENT(in ), OPTIONAL :: kt_offset ! provide fields at time other than "now" 135 ! kt_offset = -1 => fields at "before" time level 136 ! kt_offset = +1 => fields at "after" time level 137 ! etc. 138 !! 139 INTEGER :: itmp ! temporary variable 140 INTEGER :: imf ! size of the structure sd 141 INTEGER :: jf ! dummy indices 142 INTEGER :: isecend ! number of second since Jan. 1st 00h of nit000 year at nitend 143 INTEGER :: isecsbc ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step 144 INTEGER :: it_offset ! local time offset variable 145 LOGICAL :: llnxtyr ! open next year file? 146 LOGICAL :: llnxtmth ! open next month file? 147 LOGICAL :: llstop ! stop is the file does not exist 134 ! ! kt_offset = -1 => fields at "before" time level 135 ! ! kt_offset = +1 => fields at "after" time level 136 ! ! etc. 137 INTEGER :: itmp ! local variable 138 INTEGER :: imf ! size of the structure sd 139 INTEGER :: jf ! dummy indices 140 INTEGER :: isecend ! number of second since Jan. 1st 00h of nit000 year at nitend 141 INTEGER :: isecsbc ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step 142 INTEGER :: it_offset ! local time offset variable 143 LOGICAL :: llnxtyr ! open next year file? 144 LOGICAL :: llnxtmth ! open next month file? 145 LOGICAL :: llstop ! stop is the file does not exist 148 146 LOGICAL :: ll_firstcall ! true if this is the first call to fld_read for this set of fields 149 REAL(wp) :: ztinta ! ratio applied to after records when doing time interpolation150 REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation151 CHARACTER(LEN=1000) :: clfmt 152 TYPE(MAP_POINTER) :: imap ! global-to-local mapping indices147 REAL(wp) :: ztinta ! ratio applied to after records when doing time interpolation 148 REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation 149 CHARACTER(LEN=1000) :: clfmt ! write format 150 TYPE(MAP_POINTER) :: imap ! global-to-local mapping indices 153 151 !!--------------------------------------------------------------------- 154 152 ll_firstcall = kt == nit000 … … 166 164 isecsbc = nsec_year + nsec1jan000 + (kit+it_offset)*NINT( rdt/REAL(nn_baro,wp) ) 167 165 ELSE ! middle of sbc time step 168 isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdt tra(1)) + it_offset * NINT(rdttra(1))166 isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdt) + it_offset * NINT(rdt) 169 167 ENDIF 170 168 imf = SIZE( sd ) … … 193 191 CALL fld_rec( kn_fsbc, sd(jf), kt_offset = it_offset, kit = kit ) ! update after record informations 194 192 195 ! if kn_fsbc*rdt trais larger than nfreqh (which is kind of odd),193 ! if kn_fsbc*rdt is larger than nfreqh (which is kind of odd), 196 194 ! it is possible that the before value is no more the good one... we have to re-read it 197 195 ! if before is not the last record of the file currently opened and after is the first record to be read … … 214 212 IF( sd(jf)%ln_tint ) THEN 215 213 216 ! if kn_fsbc*rdt trais larger than nfreqh (which is kind of odd),214 ! if kn_fsbc*rdt is larger than nfreqh (which is kind of odd), 217 215 ! it is possible that the before value is no more the good one... we have to re-read it 218 216 ! if before record is not just just before the after record... … … 245 243 ! year/month/week/day file to be not present. If the run continue further than the current 246 244 ! year/month/week/day, next year/month/week/day file must exist 247 isecend = nsec_year + nsec1jan000 + (nitend - kt) * NINT(rdt tra(1)) ! second at the end of the run245 isecend = nsec_year + nsec1jan000 + (nitend - kt) * NINT(rdt) ! second at the end of the run 248 246 llstop = isecend > sd(jf)%nrec_a(2) ! read more than 1 record of next year 249 247 ! we suppose that the date of next file is next day (should be ok even for weekly files...) … … 299 297 END DO ! --- end loop over field --- ! 300 298 ! 301 ! ! ====================================== ! 302 ENDIF ! update field at each kn_fsbc time-step ! 303 ! ! ====================================== ! 299 ENDIF 304 300 ! 305 301 END SUBROUTINE fld_read … … 333 329 llprevday = .FALSE. 334 330 isec_week = 0 335 331 ! 336 332 ! define record informations 337 333 CALL fld_rec( kn_fsbc, sdjf, ldbefore = .TRUE. ) ! return before values in sdjf%nrec_a (as we will swap it later) 338 334 ! 339 335 ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar 340 336 ! 341 337 IF( sdjf%ln_tint ) THEN ! we need to read the previous record and we will put it in the current record structure 342 338 ! 343 339 IF( sdjf%nrec_a(1) == 0 ) THEN ! we redefine record sdjf%nrec_a(1) with the last record of previous year file 344 340 IF ( sdjf%nfreqh == -12 ) THEN ! yearly mean … … 391 387 ! 392 388 CALL fld_clopn( sdjf, iyear, imonth, iday, .NOT. llprev ) 393 389 ! 394 390 ! if previous year/month/day file does not exist, we switch to the current year/month/day 395 391 IF( llprev .AND. sdjf%num <= 0 ) THEN … … 401 397 CALL fld_clopn( sdjf ) 402 398 ENDIF 403 399 ! 404 400 IF( llprev ) THEN ! check if the record sdjf%nrec_a(1) exists in the file 405 401 idvar = iom_varid( sdjf%num, sdjf%clvar ) ! id of the variable sdjf%clvar … … 408 404 sdjf%nrec_a(1) = MIN( sdjf%nrec_a(1), inrec ) ! make sure we select an existing record 409 405 ENDIF 410 411 ! read before data in after arrays(as we will swap it later) 412 CALL fld_get( sdjf, map ) 413 406 ! 407 CALL fld_get( sdjf, map ) ! read before data in after arrays(as we will swap it later) 408 ! 414 409 clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')" 415 410 IF(lwp) WRITE(numout, clfmt) TRIM(sdjf%clvar), sdjf%nrec_a(1), REAL(sdjf%nrec_a(2),wp)/rday 416 411 ! 417 412 ENDIF 418 413 ! … … 435 430 LOGICAL , INTENT(in ), OPTIONAL :: ldbefore ! sent back before record values (default = .FALSE.) 436 431 INTEGER , INTENT(in ), OPTIONAL :: kit ! index of barotropic subcycle 437 432 ! ! used only if sdjf%ln_tint = .TRUE. 438 433 INTEGER , INTENT(in ), OPTIONAL :: kt_offset ! Offset of required time level compared to "now" 439 440 ! !434 ! ! time level in units of time steps. 435 ! 441 436 LOGICAL :: llbefore ! local definition of ldbefore 442 437 INTEGER :: iendrec ! end of this record (in seconds) … … 459 454 IF( PRESENT(kt_offset) ) it_offset = kt_offset 460 455 IF( PRESENT(kit) ) THEN ; it_offset = ( kit + it_offset ) * NINT( rdt/REAL(nn_baro,wp) ) 461 ELSE ; it_offset = it_offset * NINT( rdt tra(1))456 ELSE ; it_offset = it_offset * NINT( rdt ) 462 457 ENDIF 463 458 ! … … 536 531 ELSE ; ztmp = REAL(nsec_year ,wp) ! since 00h on Jan 1 of the current year 537 532 ENDIF 538 ztmp = ztmp + 0.5 * REAL(kn_fsbc - 1, wp) * rdt tra(1) + REAL( it_offset, wp )! centrered in the middle of sbc time step539 ztmp = ztmp + 0.01 * rdt tra(1)! avoid truncation error533 ztmp = ztmp + 0.5 * REAL(kn_fsbc - 1, wp) * rdt + REAL( it_offset, wp ) ! centrered in the middle of sbc time step 534 ztmp = ztmp + 0.01 * rdt ! avoid truncation error 540 535 IF( sdjf%ln_tint ) THEN ! time interpolation, shift by 1/2 record 541 536 ! … … 592 587 !! ** Purpose : read the data 593 588 !!---------------------------------------------------------------------- 594 TYPE(FLD) , INTENT(inout) :: sdjf ! input field related variables595 TYPE(MAP_POINTER), INTENT(in) :: map! global-to-local mapping indices596 ! !597 INTEGER :: ipk! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk )598 INTEGER :: iw! index into wgts array599 INTEGER :: ipdom! index of the domain600 INTEGER :: idvar! variable ID601 INTEGER :: idmspc! number of spatial dimensions602 LOGICAL :: lmoor! C1D case: point data589 TYPE(FLD) , INTENT(inout) :: sdjf ! input field related variables 590 TYPE(MAP_POINTER), INTENT(in ) :: map ! global-to-local mapping indices 591 ! 592 INTEGER :: ipk ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk ) 593 INTEGER :: iw ! index into wgts array 594 INTEGER :: ipdom ! index of the domain 595 INTEGER :: idvar ! variable ID 596 INTEGER :: idmspc ! number of spatial dimensions 597 LOGICAL :: lmoor ! C1D case: point data 603 598 !!--------------------------------------------------------------------- 604 599 ! … … 611 606 ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN 612 607 CALL wgt_list( sdjf, iw ) 613 IF( sdjf%ln_tint ) THEN ; CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk , sdjf%fdta(:,:,:,2),&614 &sdjf%nrec_a(1), sdjf%lsmname )615 ELSE ; CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk , sdjf%fnow(:,:,: ),&616 &sdjf%nrec_a(1), sdjf%lsmname )608 IF( sdjf%ln_tint ) THEN ; CALL fld_interp( sdjf%num, sdjf%clvar, iw, ipk, sdjf%fdta(:,:,:,2), & 609 & sdjf%nrec_a(1), sdjf%lsmname ) 610 ELSE ; CALL fld_interp( sdjf%num, sdjf%clvar, iw, ipk, sdjf%fnow(:,:,: ), & 611 & sdjf%nrec_a(1), sdjf%lsmname ) 617 612 ENDIF 618 613 ELSE 619 IF( SIZE(sdjf%fnow, 1) == jpi ) THEN ;ipdom = jpdom_data620 ELSE ;ipdom = jpdom_unknown614 IF( SIZE(sdjf%fnow, 1) == jpi ) THEN ; ipdom = jpdom_data 615 ELSE ; ipdom = jpdom_unknown 621 616 ENDIF 622 617 ! C1D case: If product of spatial dimensions == ipk, then x,y are of 623 618 ! size 1 (point/mooring data): this must be read onto the central grid point 624 619 idvar = iom_varid( sdjf%num, sdjf%clvar ) 625 idmspc = iom_file ( sdjf%num )%ndims( idvar )620 idmspc = iom_file ( sdjf%num )%ndims( idvar ) 626 621 IF( iom_file( sdjf%num )%luld( idvar ) ) idmspc = idmspc - 1 627 lmoor = ( idmspc == 0 .OR. PRODUCT( iom_file( sdjf%num )%dimsz( 1:MAX(idmspc,1) ,idvar ) ) == ipk)622 lmoor = ( idmspc == 0 .OR. PRODUCT( iom_file( sdjf%num )%dimsz( 1:MAX(idmspc,1) ,idvar ) ) == ipk ) 628 623 ! 629 624 SELECT CASE( ipk ) … … 660 655 ! 661 656 sdjf%rotn(2) = .false. ! vector not yet rotated 662 657 ! 663 658 END SUBROUTINE fld_get 659 664 660 665 661 SUBROUTINE fld_map( num, clvar, dta, nrec, map ) … … 688 684 REAL(wp), POINTER, DIMENSION(:,:,:) :: dta_read ! work space for global data 689 685 !!--------------------------------------------------------------------- 690 686 ! 691 687 ipi = SIZE( dta, 1 ) 692 688 ipj = 1 693 689 ipk = SIZE( dta, 3 ) 694 690 ! 695 691 idvar = iom_varid( num, clvar ) 696 692 ilendta = iom_file(num)%dimsz(1,idvar) … … 698 694 #if defined key_bdy 699 695 ipj = iom_file(num)%dimsz(2,idvar) 700 IF ( map%ll_unstruc) THEN! unstructured open boundary data file696 IF( map%ll_unstruc) THEN ! unstructured open boundary data file 701 697 dta_read => dta_global 702 ELSE ! structured open boundary data file698 ELSE ! structured open boundary data file 703 699 dta_read => dta_global2 704 700 ENDIF 705 701 #endif 706 702 707 IF(lwp) WRITE(numout,*) 'Dim size for ', TRIM(clvar),' is ', ilendta703 IF(lwp) WRITE(numout,*) 'Dim size for ', TRIM(clvar),' is ', ilendta 708 704 IF(lwp) WRITE(numout,*) 'Number of levels for ',TRIM(clvar),' is ', ipk 709 705 ! 710 706 SELECT CASE( ipk ) 711 707 CASE(1) ; CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1 ), nrec ) … … 713 709 END SELECT 714 710 ! 715 IF ( map%ll_unstruc ) THEN! unstructured open boundary data file711 IF( map%ll_unstruc ) THEN ! unstructured open boundary data file 716 712 DO ib = 1, ipi 717 713 DO ik = 1, ipk … … 728 724 END DO 729 725 ENDIF 730 726 ! 731 727 END SUBROUTINE fld_map 732 728 … … 738 734 !! ** Purpose : Vector fields may need to be rotated onto the local grid direction 739 735 !!---------------------------------------------------------------------- 740 INTEGER , INTENT(in ) :: kt ! ocean time step 741 TYPE(FLD), INTENT(inout), DIMENSION(:) :: sd ! input field related variables 742 !! 743 INTEGER :: ju,jv,jk,jn ! loop indices 744 INTEGER :: imf ! size of the structure sd 745 INTEGER :: ill ! character length 746 INTEGER :: iv ! indice of V component 736 INTEGER , INTENT(in ) :: kt ! ocean time step 737 TYPE(FLD), DIMENSION(:), INTENT(inout) :: sd ! input field related variables 738 ! 739 INTEGER :: ju, jv, jk, jn ! loop indices 740 INTEGER :: imf ! size of the structure sd 741 INTEGER :: ill ! character length 742 INTEGER :: iv ! indice of V component 743 CHARACTER (LEN=100) :: clcomp ! dummy weight name 747 744 REAL(wp), POINTER, DIMENSION(:,:) :: utmp, vtmp ! temporary arrays for vector rotation 748 CHARACTER (LEN=100) :: clcomp ! dummy weight name 749 !!--------------------------------------------------------------------- 750 751 CALL wrk_alloc( jpi,jpj, utmp, vtmp ) 752 745 !!--------------------------------------------------------------------- 746 ! 747 CALL wrk_alloc( jpi,jpj, utmp, vtmp ) 748 ! 753 749 !! (sga: following code should be modified so that pairs arent searched for each time 754 750 ! … … 786 782 END DO 787 783 ! 788 CALL wrk_dealloc( jpi,jpj, utmp, vtmp )784 CALL wrk_dealloc( jpi,jpj, utmp, vtmp ) 789 785 ! 790 786 END SUBROUTINE fld_rot … … 802 798 INTEGER, OPTIONAL, INTENT(in ) :: kday ! day value 803 799 LOGICAL, OPTIONAL, INTENT(in ) :: ldstop ! stop if open to read a non-existing file (default = .TRUE.) 804 ! !800 ! 805 801 LOGICAL :: llprevyr ! are we reading previous year file? 806 802 LOGICAL :: llprevmth ! are we reading previous month file? … … 853 849 ! 854 850 IF( TRIM(clname) /= TRIM(sdjf%clname) .OR. sdjf%num == 0 ) THEN ! new file to be open 855 851 ! 856 852 sdjf%clname = TRIM(clname) 857 853 IF( sdjf%num /= 0 ) CALL iom_close( sdjf%num ) ! close file if already open 858 854 CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof = LEN(TRIM(sdjf%wgtname)) > 0 ) 859 855 ! 860 856 ! find the last record to be read -> update sdjf%nreclast 861 857 indexyr = iyear - nyear + 1 … … 866 862 CASE ( 2 ) ; imonth_len = 31 ! next year -> imonth = 1 867 863 END SELECT 868 864 ! 869 865 ! last record to be read in the current file 870 866 IF ( sdjf%nfreqh == -12 ) THEN ; sdjf%nreclast = 1 ! yearly mean … … 880 876 ENDIF 881 877 ENDIF 882 878 ! 883 879 ENDIF 884 880 ! … … 901 897 INTEGER :: jf ! dummy indices 902 898 !!--------------------------------------------------------------------- 903 899 ! 904 900 DO jf = 1, SIZE(sdf) 905 901 sdf(jf)%clrootname = TRIM( cdir )//TRIM( sdf_n(jf)%clname ) … … 923 919 sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn 924 920 END DO 925 921 ! 926 922 IF(lwp) THEN ! control print 927 923 WRITE(numout,*) … … 943 939 END DO 944 940 ENDIF 945 941 ! 946 942 END SUBROUTINE fld_fill 947 943 … … 958 954 TYPE( FLD ), INTENT(in ) :: sd ! field with name of weights file 959 955 INTEGER , INTENT(inout) :: kwgt ! index of weights 960 ! !956 ! 961 957 INTEGER :: kw, nestid ! local integer 962 958 LOGICAL :: found ! local logical … … 966 962 !! weights filename is either present or we hit the end of the list 967 963 found = .FALSE. 968 964 ! 969 965 !! because agrif nest part of filenames are now added in iom_open 970 966 !! to distinguish between weights files on the different grids, need to track … … 1028 1024 TYPE( FLD ), INTENT(in) :: sd ! field with name of weights file 1029 1025 !! 1030 INTEGER :: jn ! dummy loop indices 1031 INTEGER :: inum ! temporary logical unit 1032 INTEGER :: id ! temporary variable id 1033 INTEGER :: ipk ! temporary vertical dimension 1034 CHARACTER (len=5) :: aname 1026 INTEGER :: jn ! dummy loop indices 1027 INTEGER :: inum ! local logical unit 1028 INTEGER :: id ! local variable id 1029 INTEGER :: ipk ! local vertical dimension 1030 INTEGER :: zwrap ! local integer 1031 LOGICAL :: cyclical ! 1032 CHARACTER (len=5) :: aname ! 1035 1033 INTEGER , DIMENSION(:), ALLOCATABLE :: ddims 1036 1034 INTEGER , POINTER, DIMENSION(:,:) :: data_src 1037 1035 REAL(wp), POINTER, DIMENSION(:,:) :: data_tmp 1038 LOGICAL :: cyclical 1039 INTEGER :: zwrap ! local integer 1040 !!---------------------------------------------------------------------- 1041 ! 1042 CALL wrk_alloc( jpi,jpj, data_src ) ! integer 1043 CALL wrk_alloc( jpi,jpj, data_tmp ) 1036 !!---------------------------------------------------------------------- 1037 ! 1038 CALL wrk_alloc( jpi,jpj, data_src ) ! integer 1039 CALL wrk_alloc( jpi,jpj, data_tmp ) 1044 1040 ! 1045 1041 IF( nxt_wgt > tot_wgts ) THEN … … 1151 1147 ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) ) 1152 1148 IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col(1,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) ) 1153 1149 ! 1154 1150 nxt_wgt = nxt_wgt + 1 1155 1151 ! 1156 1152 ELSE 1157 1153 CALL ctl_stop( ' fld_weight : unable to read the file ' ) … … 1166 1162 1167 1163 1168 SUBROUTINE apply_seaoverland( clmaskfile,zfieldo,jpi1_lsm,jpi2_lsm,jpj1_lsm,&1169 & jpj2_lsm, itmpi,itmpj,itmpz,rec1_lsm,recn_lsm)1164 SUBROUTINE apply_seaoverland( clmaskfile, zfieldo, jpi1_lsm, jpi2_lsm, jpj1_lsm, & 1165 & jpj2_lsm, itmpi, itmpj, itmpz, rec1_lsm, recn_lsm ) 1170 1166 !!--------------------------------------------------------------------- 1171 1167 !! *** ROUTINE apply_seaoverland *** … … 1176 1172 !! D. Delrosso INGV 1177 1173 !!---------------------------------------------------------------------- 1178 INTEGER :: inum,jni,jnj,jnz,jc ! temporary indices1179 INTEGER, INTENT(in ) :: itmpi,itmpj,itmpz ! lengths1180 INTEGER, INTENT(in) :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices1181 INTEGER, DIMENSION(3), INTENT(in) :: rec1_lsm,recn_lsm ! temporary arrays for start and length1182 REAL(wp),DIMENSION (:,:,:),INTENT(inout) :: zfieldo ! input/output array for seaoverland application1183 REAL(wp),DIMENSION (:,:,:),ALLOCATABLE :: zslmec1 ! temporary array for land point detection1184 REAL(wp),DIMENSION (:,:), ALLOCATABLE :: zfieldn ! array of forcing field with undeff for land points1185 REAL(wp),DIMENSION (:,: ), ALLOCATABLE :: zfield ! array of forcing field1186 CHARACTER (len=100), INTENT(in) :: clmaskfile ! land/sea mask file name1187 !!---------------------------------------------------------------------1188 ALLOCATE ( zslmec1(itmpi,itmpj,itmpz) )1189 ALLOCATE ( zfieldn(itmpi,itmpj) )1190 ALLOCATE ( z field(itmpi,itmpj) )1191 1174 INTEGER, INTENT(in ) :: itmpi,itmpj,itmpz ! lengths 1175 INTEGER, INTENT(in ) :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices 1176 INTEGER, DIMENSION(3), INTENT(in ) :: rec1_lsm,recn_lsm ! temporary arrays for start and length 1177 REAL(wp),DIMENSION (:,:,:),INTENT(inout) :: zfieldo ! input/output array for seaoverland application 1178 CHARACTER (len=100), INTENT(in ) :: clmaskfile ! land/sea mask file name 1179 ! 1180 INTEGER :: inum,jni,jnj,jnz,jc ! local indices 1181 REAL(wp),DIMENSION (:,:,:),ALLOCATABLE :: zslmec1 ! local array for land point detection 1182 REAL(wp),DIMENSION (:,:), ALLOCATABLE :: zfieldn ! array of forcing field with undeff for land points 1183 REAL(wp),DIMENSION (:,:), ALLOCATABLE :: zfield ! array of forcing field 1184 !!--------------------------------------------------------------------- 1185 ! 1186 ALLOCATE ( zslmec1(itmpi,itmpj,itmpz), zfieldn(itmpi,itmpj), zfield(itmpi,itmpj) ) 1187 ! 1192 1188 ! Retrieve the land sea mask data 1193 1189 CALL iom_open( clmaskfile, inum ) 1194 1190 SELECT CASE( SIZE(zfieldo(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) ) 1195 1191 CASE(1) 1196 1192 CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), 1, rec1_lsm, recn_lsm) 1197 1193 CASE DEFAULT 1198 1194 CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), 1, rec1_lsm, recn_lsm) 1199 1195 END SELECT 1200 1196 CALL iom_close( inum ) 1201 1202 DO jnz=1,rec1_lsm(3) !! Loop over k dimension 1203 1204 DO jni=1,itmpi !! copy the original field into a tmp array 1205 DO jnj=1,itmpj !! substituting undeff over land points 1206 zfieldn(jni,jnj) = zfieldo(jni,jnj,jnz) 1207 IF ( zslmec1(jni,jnj,jnz) == 1. ) THEN 1208 zfieldn(jni,jnj) = undeff_lsm 1209 ENDIF 1197 ! 1198 DO jnz=1,rec1_lsm(3) !! Loop over k dimension 1199 ! 1200 DO jni = 1, itmpi !! copy the original field into a tmp array 1201 DO jnj = 1, itmpj !! substituting undeff over land points 1202 zfieldn(jni,jnj) = zfieldo(jni,jnj,jnz) 1203 IF( zslmec1(jni,jnj,jnz) == 1. ) zfieldn(jni,jnj) = undeff_lsm 1210 1204 END DO 1211 1205 END DO 1212 1213 CALL seaoverland(zfieldn,itmpi,itmpj,zfield) 1214 DO jc=1,nn_lsm 1215 CALL seaoverland(zfield,itmpi,itmpj,zfield) 1216 END DO 1217 1218 ! Check for Undeff and substitute original values 1219 IF(ANY(zfield==undeff_lsm)) THEN 1220 DO jni=1,itmpi 1221 DO jnj=1,itmpj 1222 IF (zfield(jni,jnj)==undeff_lsm) THEN 1223 zfield(jni,jnj) = zfieldo(jni,jnj,jnz) 1224 ENDIF 1225 ENDDO 1226 ENDDO 1227 ENDIF 1228 1229 zfieldo(:,:,jnz)=zfield(:,:) 1230 1231 END DO !! End Loop over k dimension 1232 1233 DEALLOCATE ( zslmec1 ) 1234 DEALLOCATE ( zfieldn ) 1235 DEALLOCATE ( zfield ) 1236 1206 ! 1207 CALL seaoverland( zfieldn, itmpi, itmpj, zfield ) 1208 DO jc = 1, nn_lsm 1209 CALL seaoverland( zfield, itmpi, itmpj, zfield ) 1210 END DO 1211 ! 1212 ! Check for Undeff and substitute original values 1213 IF( ANY(zfield==undeff_lsm) ) THEN 1214 DO jni = 1, itmpi 1215 DO jnj = 1, itmpj 1216 IF( zfield(jni,jnj)==undeff_lsm ) zfield(jni,jnj) = zfieldo(jni,jnj,jnz) 1217 END DO 1218 END DO 1219 ENDIF 1220 ! 1221 zfieldo(:,:,jnz) = zfield(:,:) 1222 ! 1223 END DO !! End Loop over k dimension 1224 ! 1225 DEALLOCATE ( zslmec1, zfieldn, zfield ) 1226 ! 1237 1227 END SUBROUTINE apply_seaoverland 1238 1228 1239 1229 1240 SUBROUTINE seaoverland( zfieldn,ileni,ilenj,zfield)1230 SUBROUTINE seaoverland( zfieldn, ileni, ilenj, zfield ) 1241 1231 !!--------------------------------------------------------------------- 1242 1232 !! *** ROUTINE seaoverland *** … … 1245 1235 !! D. Delrosso INGV 1246 1236 !!---------------------------------------------------------------------- 1247 INTEGER,INTENT(in) :: ileni,ilenj ! lengths 1248 REAL,DIMENSION (ileni,ilenj),INTENT(in) :: zfieldn ! array of forcing field with undeff for land points 1249 REAL,DIMENSION (ileni,ilenj),INTENT(out) :: zfield ! array of forcing field 1250 REAL,DIMENSION (ileni,ilenj) :: zmat1,zmat2,zmat3,zmat4 ! temporary arrays for seaoverland application 1251 REAL,DIMENSION (ileni,ilenj) :: zmat5,zmat6,zmat7,zmat8 ! temporary arrays for seaoverland application 1252 REAL,DIMENSION (ileni,ilenj) :: zlsm2d ! temporary arrays for seaoverland application 1253 REAL,DIMENSION (ileni,ilenj,8) :: zlsm3d ! temporary arrays for seaoverland application 1254 LOGICAL,DIMENSION (ileni,ilenj,8) :: ll_msknan3d ! logical mask for undeff detection 1255 LOGICAL,DIMENSION (ileni,ilenj) :: ll_msknan2d ! logical mask for undeff detection 1237 INTEGER , INTENT(in ) :: ileni,ilenj ! lengths 1238 REAL, DIMENSION (ileni,ilenj), INTENT(in ) :: zfieldn ! array of forcing field with undeff for land points 1239 REAL, DIMENSION (ileni,ilenj), INTENT( out) :: zfield ! array of forcing field 1240 ! 1241 REAL , DIMENSION (ileni,ilenj) :: zmat1, zmat2, zmat3, zmat4 ! local arrays 1242 REAL , DIMENSION (ileni,ilenj) :: zmat5, zmat6, zmat7, zmat8 ! - - 1243 REAL , DIMENSION (ileni,ilenj) :: zlsm2d ! - - 1244 REAL , DIMENSION (ileni,ilenj,8) :: zlsm3d ! - - 1245 LOGICAL, DIMENSION (ileni,ilenj,8) :: ll_msknan3d ! logical mask for undeff detection 1246 LOGICAL, DIMENSION (ileni,ilenj) :: ll_msknan2d ! logical mask for undeff detection 1256 1247 !!---------------------------------------------------------------------- 1257 zmat8 = eoshift( zfieldn , SHIFT=-1, BOUNDARY = (/zfieldn(:,1)/) ,DIM=2)1258 zmat1 = eoshift( zmat8 , SHIFT=-1, BOUNDARY = (/zmat8(1,:)/) ,DIM=1)1259 zmat2 = eoshift( zfieldn , SHIFT=-1, BOUNDARY = (/zfieldn(1,:)/) ,DIM=1)1260 zmat4 = eoshift( zfieldn , SHIFT= 1, BOUNDARY = (/zfieldn(:,ilenj)/),DIM=2)1261 zmat3 = eoshift( zmat4 , SHIFT=-1, BOUNDARY = (/zmat4(1,:)/) ,DIM=1)1262 zmat5 = eoshift( zmat4 , SHIFT= 1, BOUNDARY = (/zmat4(ileni,:)/) ,DIM=1)1263 zmat6 = eoshift( zfieldn , SHIFT= 1, BOUNDARY = (/zfieldn(ileni,:)/),DIM=1)1264 zmat7 = eoshift( zmat8 , SHIFT= 1, BOUNDARY = (/zmat8(ileni,:)/) ,DIM=1)1265 1248 zmat8 = eoshift( zfieldn , SHIFT=-1 , BOUNDARY = (/zfieldn(:,1)/) , DIM=2 ) 1249 zmat1 = eoshift( zmat8 , SHIFT=-1 , BOUNDARY = (/zmat8(1,:)/) , DIM=1 ) 1250 zmat2 = eoshift( zfieldn , SHIFT=-1 , BOUNDARY = (/zfieldn(1,:)/) , DIM=1 ) 1251 zmat4 = eoshift( zfieldn , SHIFT= 1 , BOUNDARY = (/zfieldn(:,ilenj)/) , DIM=2 ) 1252 zmat3 = eoshift( zmat4 , SHIFT=-1 , BOUNDARY = (/zmat4(1,:)/) , DIM=1 ) 1253 zmat5 = eoshift( zmat4 , SHIFT= 1 , BOUNDARY = (/zmat4(ileni,:)/) , DIM=1 ) 1254 zmat6 = eoshift( zfieldn , SHIFT= 1 , BOUNDARY = (/zfieldn(ileni,:)/) , DIM=1 ) 1255 zmat7 = eoshift( zmat8 , SHIFT= 1 , BOUNDARY = (/zmat8(ileni,:)/) , DIM=1 ) 1256 ! 1266 1257 zlsm3d = RESHAPE( (/ zmat1, zmat2, zmat3, zmat4, zmat5, zmat6, zmat7, zmat8 /), (/ ileni, ilenj, 8 /)) 1267 ll_msknan3d = .not.(zlsm3d==undeff_lsm) 1268 ll_msknan2d = .not.(zfieldn==undeff_lsm) ! FALSE where is Undeff (land) 1269 zlsm2d = (SUM ( zlsm3d, 3 , ll_msknan3d ) )/(MAX(1,(COUNT( ll_msknan3d , 3 )) )) 1270 WHERE ((COUNT( ll_msknan3d , 3 )) == 0.0_wp) zlsm2d = undeff_lsm 1271 zfield = MERGE (zfieldn,zlsm2d,ll_msknan2d) 1258 ll_msknan3d = .NOT.( zlsm3d == undeff_lsm ) 1259 ll_msknan2d = .NOT.( zfieldn == undeff_lsm ) ! FALSE where is Undeff (land) 1260 zlsm2d = SUM( zlsm3d, 3 , ll_msknan3d ) / MAX( 1 , COUNT( ll_msknan3d , 3 ) ) 1261 WHERE( COUNT( ll_msknan3d , 3 ) == 0._wp ) zlsm2d = undeff_lsm 1262 zfield = MERGE( zfieldn, zlsm2d, ll_msknan2d ) 1263 ! 1272 1264 END SUBROUTINE seaoverland 1273 1265 … … 1288 1280 INTEGER , INTENT(in ) :: nrec ! record number to read (ie time slice) 1289 1281 CHARACTER(LEN=*) , INTENT(in ) :: lsmfile ! land sea mask file name 1290 !! 1291 REAL(wp),DIMENSION(:,:,:),ALLOCATABLE :: ztmp_fly_dta,zfieldo ! temporary array of values on input grid 1292 INTEGER, DIMENSION(3) :: rec1,recn ! temporary arrays for start and length 1293 INTEGER, DIMENSION(3) :: rec1_lsm,recn_lsm ! temporary arrays for start and length in case of seaoverland 1294 INTEGER :: ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2 ! temporary indices 1295 INTEGER :: jk, jn, jm, jir, jjr ! loop counters 1296 INTEGER :: ni, nj ! lengths 1297 INTEGER :: jpimin,jpiwid ! temporary indices 1298 INTEGER :: jpimin_lsm,jpiwid_lsm ! temporary indices 1299 INTEGER :: jpjmin,jpjwid ! temporary indices 1300 INTEGER :: jpjmin_lsm,jpjwid_lsm ! temporary indices 1301 INTEGER :: jpi1,jpi2,jpj1,jpj2 ! temporary indices 1302 INTEGER :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices 1303 INTEGER :: itmpi,itmpj,itmpz ! lengths 1304 1282 ! 1283 INTEGER, DIMENSION(3) :: rec1, recn ! temporary arrays for start and length 1284 INTEGER, DIMENSION(3) :: rec1_lsm, recn_lsm ! temporary arrays for start and length in case of seaoverland 1285 INTEGER :: ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2 ! temporary indices 1286 INTEGER :: jk, jn, jm, jir, jjr ! loop counters 1287 INTEGER :: ni, nj ! lengths 1288 INTEGER :: jpimin,jpiwid ! temporary indices 1289 INTEGER :: jpimin_lsm,jpiwid_lsm ! temporary indices 1290 INTEGER :: jpjmin,jpjwid ! temporary indices 1291 INTEGER :: jpjmin_lsm,jpjwid_lsm ! temporary indices 1292 INTEGER :: jpi1,jpi2,jpj1,jpj2 ! temporary indices 1293 INTEGER :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices 1294 INTEGER :: itmpi,itmpj,itmpz ! lengths 1295 REAL(wp),DIMENSION(:,:,:), ALLOCATABLE :: ztmp_fly_dta ! local array of values on input grid 1305 1296 !!---------------------------------------------------------------------- 1306 1297 ! … … 1355 1346 1356 1347 1357 itmpi= SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),1)1358 itmpj= SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),2)1348 itmpi=jpi2_lsm-jpi1_lsm+1 1349 itmpj=jpj2_lsm-jpj1_lsm+1 1359 1350 itmpz=kk 1360 1351 ALLOCATE(ztmp_fly_dta(itmpi,itmpj,itmpz)) -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90
r5836 r6140 5 5 !!====================================================================== 6 6 !! History : OPA ! 07-1996 (O. Marti) Original code 7 !! NEMO 1.0 ! 02-2008 (G. Madec) F90: Free form 8 !! 3.0 ! 7 !! NEMO 1.0 ! 06-2006 (G. Madec ) Free form, F90 + opt. 8 !! ! 04-2007 (S. Masson) angle: Add T, F points and bugfix in cos lateral boundary 9 !! 3.0 ! 07-2008 (G. Madec) geo2oce suppress lon/lat agruments 10 !! 3.7 ! 11-2015 (G. Madec) remove the unused repere and repcmo routines 9 11 !!---------------------------------------------------------------------- 10 12 11 13 !!---------------------------------------------------------------------- 12 !! r epcmo :13 !! angle :14 !! geo2oce :15 !! repere : old routine suppress it ???14 !! rot_rep : Rotate the Repere: geographic grid <==> stretched coordinates grid 15 !! angle : 16 !! geo2oce : 17 !! oce2geo : 16 18 !!---------------------------------------------------------------------- 17 USE dom_oce ! mesh and scale factors 18 USE phycst ! physical constants 19 USE in_out_manager ! I/O manager 20 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 21 USE lib_mpp ! MPP library 19 USE dom_oce ! mesh and scale factors 20 USE phycst ! physical constants 21 ! 22 USE in_out_manager ! I/O manager 23 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 24 USE lib_mpp ! MPP library 22 25 23 26 IMPLICIT NONE 24 27 PRIVATE 25 28 26 PUBLIC rot_rep, repcmo, repere, geo2oce, oce2geo ! only rot_rep should be used 27 ! repcmo and repere are keep only for compatibility. 28 ! they are only a useless overlay of rot_rep 29 30 PUBLIC obs_rot 31 32 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: & 33 gsint, gcost, & ! cos/sin between model grid lines and NP direction at T point 34 gsinu, gcosu, & ! cos/sin between model grid lines and NP direction at U point 35 gsinv, gcosv, & ! cos/sin between model grid lines and NP direction at V point 36 gsinf, gcosf ! cos/sin between model grid lines and NP direction at F point 29 PUBLIC rot_rep ! called in sbccpl, fldread, and cyclone 30 PUBLIC geo2oce ! called in sbccpl 31 PUBLIC oce2geo ! called in sbccpl 32 PUBLIC obs_rot ! called in obs_rot_vel and obs_write 33 34 ! ! cos/sin between model grid lines and NP direction 35 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: gsint, gcost ! at T point 36 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: gsinu, gcosu ! at U point 37 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: gsinv, gcosv ! at V point 38 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: gsinf, gcosf ! at F point 37 39 38 40 LOGICAL , SAVE, DIMENSION(4) :: linit = .FALSE. 39 41 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gsinlon, gcoslon, gsinlat, gcoslat 40 42 41 LOGICAL :: lmust_init = .TRUE. !: used to initialize the cos/sin variables (se above)43 LOGICAL :: lmust_init = .TRUE. !: used to initialize the cos/sin variables (see above) 42 44 43 45 !! * Substitutions … … 50 52 CONTAINS 51 53 52 SUBROUTINE repcmo ( pxu1, pyu1, pxv1, pyv1, &53 px2 , py2 )54 !!----------------------------------------------------------------------55 !! *** ROUTINE repcmo ***56 !!57 !! ** Purpose : Change vector componantes from a geographic grid to a58 !! stretched coordinates grid.59 !!60 !! ** Method : Initialization of arrays at the first call.61 !!62 !! ** Action : - px2 : first componante (defined at u point)63 !! - py2 : second componante (defined at v point)64 !!----------------------------------------------------------------------65 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pxu1, pyu1 ! geographic vector componantes at u-point66 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pxv1, pyv1 ! geographic vector componantes at v-point67 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: px2 ! i-componante (defined at u-point)68 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: py2 ! j-componante (defined at v-point)69 !!----------------------------------------------------------------------70 71 ! Change from geographic to stretched coordinate72 ! ----------------------------------------------73 CALL rot_rep( pxu1, pyu1, 'U', 'en->i',px2 )74 CALL rot_rep( pxv1, pyv1, 'V', 'en->j',py2 )75 76 END SUBROUTINE repcmo77 78 79 54 SUBROUTINE rot_rep ( pxin, pyin, cd_type, cdtodo, prot ) 80 55 !!---------------------------------------------------------------------- … … 83 58 !! ** Purpose : Rotate the Repere: Change vector componantes between 84 59 !! geographic grid <--> stretched coordinates grid. 85 !! 86 !! History : 87 !! 9.2 ! 07-04 (S. Masson) 88 !! (O. Marti ) Original code (repere and repcmo) 89 !!---------------------------------------------------------------------- 90 REAL(wp), DIMENSION(jpi,jpj), INTENT( IN ) :: pxin, pyin ! vector componantes 91 CHARACTER(len=1), INTENT( IN ) :: cd_type ! define the nature of pt2d array grid-points 92 CHARACTER(len=5), INTENT( IN ) :: cdtodo ! specify the work to do: 93 !! ! 'en->i' east-north componantes to model i componante 94 !! ! 'en->j' east-north componantes to model j componante 95 !! ! 'ij->e' model i-j componantes to east componante 96 !! ! 'ij->n' model i-j componantes to east componante 97 REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: prot 98 !!---------------------------------------------------------------------- 99 100 ! Initialization of gsin* and gcos* at first call 101 ! ----------------------------------------------- 102 103 IF( lmust_init ) THEN 60 !!---------------------------------------------------------------------- 61 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pxin, pyin ! vector componantes 62 CHARACTER(len=1), INTENT(in ) :: cd_type ! define the nature of pt2d array grid-points 63 CHARACTER(len=5), INTENT(in ) :: cdtodo ! type of transpormation: 64 ! ! 'en->i' = east-north to i-component 65 ! ! 'en->j' = east-north to j-component 66 ! ! 'ij->e' = (i,j) components to east 67 ! ! 'ij->n' = (i,j) components to north 68 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: prot 69 !!---------------------------------------------------------------------- 70 ! 71 IF( lmust_init ) THEN ! at 1st call only: set gsin. & gcos. 104 72 IF(lwp) WRITE(numout,*) 105 IF(lwp) WRITE(numout,*) ' rot_rep : geographic <--> stretched'106 IF(lwp) WRITE(numout,*) ' ~~~~~ coordinate transformation'73 IF(lwp) WRITE(numout,*) ' rot_rep: coordinate transformation : geographic <==> model (i,j)-components' 74 IF(lwp) WRITE(numout,*) ' ~~~~~~~~ ' 107 75 ! 108 76 CALL angle ! initialization of the transformation 109 77 lmust_init = .FALSE. 110 78 ENDIF 111 112 SELECT CASE (cdtodo) 113 CASE ('en->i') ! 'en->i' est-north componantes to model i componante 79 ! 80 SELECT CASE( cdtodo ) ! type of rotation 81 ! 82 CASE( 'en->i' ) ! east-north to i-component 114 83 SELECT CASE (cd_type) 115 84 CASE ('T') ; prot(:,:) = pxin(:,:) * gcost(:,:) + pyin(:,:) * gsint(:,:) … … 119 88 CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 120 89 END SELECT 121 CASE ('en->j') ! 'en->j' est-north componantes to model j componante90 CASE ('en->j') ! east-north to j-component 122 91 SELECT CASE (cd_type) 123 92 CASE ('T') ; prot(:,:) = pyin(:,:) * gcost(:,:) - pxin(:,:) * gsint(:,:) … … 127 96 CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 128 97 END SELECT 129 CASE ('ij->e') ! 'ij->e' model i-j componantes to est componante98 CASE ('ij->e') ! (i,j)-components to east 130 99 SELECT CASE (cd_type) 131 100 CASE ('T') ; prot(:,:) = pxin(:,:) * gcost(:,:) - pyin(:,:) * gsint(:,:) … … 135 104 CASE DEFAULT ; CALL ctl_stop( 'Only T, U, V and F grid points are coded' ) 136 105 END SELECT 137 CASE ('ij->n') ! 'ij->n' model i-j componantes to est componante106 CASE ('ij->n') ! (i,j)-components to north 138 107 SELECT CASE (cd_type) 139 108 CASE ('T') ; prot(:,:) = pyin(:,:) * gcost(:,:) + pxin(:,:) * gsint(:,:) … … 144 113 END SELECT 145 114 CASE DEFAULT ; CALL ctl_stop( 'rot_rep: Syntax Error in the definition of cdtodo' ) 115 ! 146 116 END SELECT 147 117 ! 148 118 END SUBROUTINE rot_rep 149 119 … … 155 125 !! ** Purpose : Compute angles between model grid lines and the North direction 156 126 !! 157 !! ** Method : 158 !! 159 !! ** Action : Compute (gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf) arrays: 160 !! sinus and cosinus of the angle between the north-south axe and the 161 !! j-direction at t, u, v and f-points 162 !! 163 !! History : 164 !! 7.0 ! 96-07 (O. Marti ) Original code 165 !! 8.0 ! 98-06 (G. Madec ) 166 !! 8.5 ! 98-06 (G. Madec ) Free form, F90 + opt. 167 !! 9.2 ! 07-04 (S. Masson) Add T, F points and bugfix in cos lateral boundary 168 !!---------------------------------------------------------------------- 169 INTEGER :: ji, jj ! dummy loop indices 170 INTEGER :: ierr ! local integer 171 REAL(wp) :: & 172 zlam, zphi, & ! temporary scalars 173 zlan, zphh, & ! " " 174 zxnpt, zynpt, znnpt, & ! x,y components and norm of the vector: T point to North Pole 175 zxnpu, zynpu, znnpu, & ! x,y components and norm of the vector: U point to North Pole 176 zxnpv, zynpv, znnpv, & ! x,y components and norm of the vector: V point to North Pole 177 zxnpf, zynpf, znnpf, & ! x,y components and norm of the vector: F point to North Pole 178 zxvvt, zyvvt, znvvt, & ! x,y components and norm of the vector: between V points below and above a T point 179 zxffu, zyffu, znffu, & ! x,y components and norm of the vector: between F points below and above a U point 180 zxffv, zyffv, znffv, & ! x,y components and norm of the vector: between F points left and right a V point 181 zxuuf, zyuuf, znuuf ! x,y components and norm of the vector: between U points below and above a F point 182 !!---------------------------------------------------------------------- 183 127 !! ** Method : sinus and cosinus of the angle between the north-south axe 128 !! and the j-direction at t, u, v and f-points 129 !! dot and cross products are used to obtain cos and sin, resp. 130 !! 131 !! ** Action : - gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf 132 !!---------------------------------------------------------------------- 133 INTEGER :: ji, jj ! dummy loop indices 134 INTEGER :: ierr ! local integer 135 REAL(wp) :: zlam, zphi ! local scalars 136 REAL(wp) :: zlan, zphh ! - - 137 REAL(wp) :: zxnpt, zynpt, znnpt ! x,y components and norm of the vector: T point to North Pole 138 REAL(wp) :: zxnpu, zynpu, znnpu ! x,y components and norm of the vector: U point to North Pole 139 REAL(wp) :: zxnpv, zynpv, znnpv ! x,y components and norm of the vector: V point to North Pole 140 REAL(wp) :: zxnpf, zynpf, znnpf ! x,y components and norm of the vector: F point to North Pole 141 REAL(wp) :: zxvvt, zyvvt, znvvt ! x,y components and norm of the vector: between V points below and above a T point 142 REAL(wp) :: zxffu, zyffu, znffu ! x,y components and norm of the vector: between F points below and above a U point 143 REAL(wp) :: zxffv, zyffv, znffv ! x,y components and norm of the vector: between F points left and right a V point 144 REAL(wp) :: zxuuf, zyuuf, znuuf ! x,y components and norm of the vector: between U points below and above a F point 145 !!---------------------------------------------------------------------- 146 ! 184 147 ALLOCATE( gsint(jpi,jpj), gcost(jpi,jpj), & 185 148 & gsinu(jpi,jpj), gcosu(jpi,jpj), & … … 187 150 & gsinf(jpi,jpj), gcosf(jpi,jpj), STAT=ierr ) 188 151 IF(lk_mpp) CALL mpp_sum( ierr ) 189 IF( ierr /= 0 ) CALL ctl_stop( 'angle: unable to allocate arrays' )190 152 IF( ierr /= 0 ) CALL ctl_stop( 'angle: unable to allocate arrays' ) 153 ! 191 154 ! ============================= ! 192 155 ! Compute the cosinus and sinus ! 193 156 ! ============================= ! 194 157 ! (computation done on the north stereographic polar plane) 195 158 ! 196 159 DO jj = 2, jpjm1 197 160 DO ji = fs_2, jpi ! vector opt. 198 199 ! north pole direction & modulous (at t-point) 200 zlam = glamt(ji,jj) 161 ! 162 zlam = glamt(ji,jj) ! north pole direction & modulous (at t-point) 201 163 zphi = gphit(ji,jj) 202 164 zxnpt = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 203 165 zynpt = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 204 166 znnpt = zxnpt*zxnpt + zynpt*zynpt 205 206 ! north pole direction & modulous (at u-point) 207 zlam = glamu(ji,jj) 167 ! 168 zlam = glamu(ji,jj) ! north pole direction & modulous (at u-point) 208 169 zphi = gphiu(ji,jj) 209 170 zxnpu = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 210 171 zynpu = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 211 172 znnpu = zxnpu*zxnpu + zynpu*zynpu 212 213 ! north pole direction & modulous (at v-point) 214 zlam = glamv(ji,jj) 173 ! 174 zlam = glamv(ji,jj) ! north pole direction & modulous (at v-point) 215 175 zphi = gphiv(ji,jj) 216 176 zxnpv = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 217 177 zynpv = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 218 178 znnpv = zxnpv*zxnpv + zynpv*zynpv 219 220 ! north pole direction & modulous (at f-point) 221 zlam = glamf(ji,jj) 179 ! 180 zlam = glamf(ji,jj) ! north pole direction & modulous (at f-point) 222 181 zphi = gphif(ji,jj) 223 182 zxnpf = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 224 183 zynpf = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 225 184 znnpf = zxnpf*zxnpf + zynpf*zynpf 226 227 ! j-direction: v-point segment direction (around t-point) 228 zlam = glamv(ji,jj ) 185 ! 186 zlam = glamv(ji,jj ) ! j-direction: v-point segment direction (around t-point) 229 187 zphi = gphiv(ji,jj ) 230 188 zlan = glamv(ji,jj-1) … … 236 194 znvvt = SQRT( znnpt * ( zxvvt*zxvvt + zyvvt*zyvvt ) ) 237 195 znvvt = MAX( znvvt, 1.e-14 ) 238 239 ! j-direction: f-point segment direction (around u-point) 240 zlam = glamf(ji,jj ) 196 ! 197 zlam = glamf(ji,jj ) ! j-direction: f-point segment direction (around u-point) 241 198 zphi = gphif(ji,jj ) 242 199 zlan = glamf(ji,jj-1) … … 248 205 znffu = SQRT( znnpu * ( zxffu*zxffu + zyffu*zyffu ) ) 249 206 znffu = MAX( znffu, 1.e-14 ) 250 251 ! i-direction: f-point segment direction (around v-point) 252 zlam = glamf(ji ,jj) 207 ! 208 zlam = glamf(ji ,jj) ! i-direction: f-point segment direction (around v-point) 253 209 zphi = gphif(ji ,jj) 254 210 zlan = glamf(ji-1,jj) … … 260 216 znffv = SQRT( znnpv * ( zxffv*zxffv + zyffv*zyffv ) ) 261 217 znffv = MAX( znffv, 1.e-14 ) 262 263 ! j-direction: u-point segment direction (around f-point) 264 zlam = glamu(ji,jj+1) 218 ! 219 zlam = glamu(ji,jj+1) ! j-direction: u-point segment direction (around f-point) 265 220 zphi = gphiu(ji,jj+1) 266 221 zlan = glamu(ji,jj ) … … 272 227 znuuf = SQRT( znnpf * ( zxuuf*zxuuf + zyuuf*zyuuf ) ) 273 228 znuuf = MAX( znuuf, 1.e-14 ) 274 275 ! cosinus and sinus using scalar and vectorialproducts229 ! 230 ! ! cosinus and sinus using dot and cross products 276 231 gsint(ji,jj) = ( zxnpt*zyvvt - zynpt*zxvvt ) / znvvt 277 232 gcost(ji,jj) = ( zxnpt*zxvvt + zynpt*zyvvt ) / znvvt 278 233 ! 279 234 gsinu(ji,jj) = ( zxnpu*zyffu - zynpu*zxffu ) / znffu 280 235 gcosu(ji,jj) = ( zxnpu*zxffu + zynpu*zyffu ) / znffu 281 236 ! 282 237 gsinf(ji,jj) = ( zxnpf*zyuuf - zynpf*zxuuf ) / znuuf 283 238 gcosf(ji,jj) = ( zxnpf*zxuuf + zynpf*zyuuf ) / znuuf 284 285 ! (caution, rotation of 90 degres) 239 ! 286 240 gsinv(ji,jj) = ( zxnpv*zxffv + zynpv*zyffv ) / znffv 287 gcosv(ji,jj) =-( zxnpv*zyffv - zynpv*zxffv ) / znffv 288 241 gcosv(ji,jj) =-( zxnpv*zyffv - zynpv*zxffv ) / znffv ! (caution, rotation of 90 degres) 242 ! 289 243 END DO 290 244 END DO … … 318 272 ! Lateral boundary conditions ! 319 273 ! =========================== ! 320 321 ! lateral boundary cond.: T-, U-, V-, F-pts, sgn 274 ! ! lateral boundary cond.: T-, U-, V-, F-pts, sgn 322 275 CALL lbc_lnk( gcost, 'T', -1. ) ; CALL lbc_lnk( gsint, 'T', -1. ) 323 276 CALL lbc_lnk( gcosu, 'U', -1. ) ; CALL lbc_lnk( gsinu, 'U', -1. ) 324 277 CALL lbc_lnk( gcosv, 'V', -1. ) ; CALL lbc_lnk( gsinv, 'V', -1. ) 325 278 CALL lbc_lnk( gcosf, 'F', -1. ) ; CALL lbc_lnk( gsinf, 'F', -1. ) 326 279 ! 327 280 END SUBROUTINE angle 328 281 329 282 330 SUBROUTINE geo2oce ( pxx, pyy, pzz, cgrid, & 331 pte, ptn ) 283 SUBROUTINE geo2oce ( pxx, pyy, pzz, cgrid, pte, ptn ) 332 284 !!---------------------------------------------------------------------- 333 285 !! *** ROUTINE geo2oce *** … … 335 287 !! ** Purpose : 336 288 !! 337 !! ** Method : Change wind stress from geocentric to east/north 338 !! 339 !! History : 340 !! ! (O. Marti) Original code 341 !! ! 91-03 (G. Madec) 342 !! ! 92-07 (M. Imbard) 343 !! ! 99-11 (M. Imbard) NetCDF format with IOIPSL 344 !! ! 00-08 (D. Ludicone) Reduced section at Bab el Mandeb 345 !! 8.5 ! 02-06 (G. Madec) F90: Free form 346 !! 3.0 ! 07-08 (G. Madec) geo2oce suppress lon/lat agruments 289 !! ** Method : Change a vector from geocentric to east/north 290 !! 347 291 !!---------------------------------------------------------------------- 348 292 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pxx, pyy, pzz 349 293 CHARACTER(len=1) , INTENT(in ) :: cgrid 350 294 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pte, ptn 351 ! !295 ! 352 296 REAL(wp), PARAMETER :: rpi = 3.141592653e0 353 297 REAL(wp), PARAMETER :: rad = rpi / 180.e0 … … 355 299 INTEGER :: ierr ! local integer 356 300 !!---------------------------------------------------------------------- 357 301 ! 358 302 IF( .NOT. ALLOCATED( gsinlon ) ) THEN 359 303 ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) , & … … 361 305 IF( lk_mpp ) CALL mpp_sum( ierr ) 362 306 IF( ierr /= 0 ) CALL ctl_stop('geo2oce: unable to allocate arrays' ) 307 ENDIF 308 ! 309 SELECT CASE( cgrid) 310 CASE ( 'T' ) 311 ig = 1 312 IF( .NOT. linit(ig) ) THEN 313 gsinlon(:,:,ig) = SIN( rad * glamt(:,:) ) 314 gcoslon(:,:,ig) = COS( rad * glamt(:,:) ) 315 gsinlat(:,:,ig) = SIN( rad * gphit(:,:) ) 316 gcoslat(:,:,ig) = COS( rad * gphit(:,:) ) 317 linit(ig) = .TRUE. 318 ENDIF 319 CASE ( 'U' ) 320 ig = 2 321 IF( .NOT. linit(ig) ) THEN 322 gsinlon(:,:,ig) = SIN( rad * glamu(:,:) ) 323 gcoslon(:,:,ig) = COS( rad * glamu(:,:) ) 324 gsinlat(:,:,ig) = SIN( rad * gphiu(:,:) ) 325 gcoslat(:,:,ig) = COS( rad * gphiu(:,:) ) 326 linit(ig) = .TRUE. 327 ENDIF 328 CASE ( 'V' ) 329 ig = 3 330 IF( .NOT. linit(ig) ) THEN 331 gsinlon(:,:,ig) = SIN( rad * glamv(:,:) ) 332 gcoslon(:,:,ig) = COS( rad * glamv(:,:) ) 333 gsinlat(:,:,ig) = SIN( rad * gphiv(:,:) ) 334 gcoslat(:,:,ig) = COS( rad * gphiv(:,:) ) 335 linit(ig) = .TRUE. 336 ENDIF 337 CASE ( 'F' ) 338 ig = 4 339 IF( .NOT. linit(ig) ) THEN 340 gsinlon(:,:,ig) = SIN( rad * glamf(:,:) ) 341 gcoslon(:,:,ig) = COS( rad * glamf(:,:) ) 342 gsinlat(:,:,ig) = SIN( rad * gphif(:,:) ) 343 gcoslat(:,:,ig) = COS( rad * gphif(:,:) ) 344 linit(ig) = .TRUE. 345 ENDIF 346 CASE default 347 WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid 348 CALL ctl_stop( ctmp1 ) 349 END SELECT 350 ! 351 pte = - gsinlon(:,:,ig) * pxx + gcoslon(:,:,ig) * pyy 352 ptn = - gcoslon(:,:,ig) * gsinlat(:,:,ig) * pxx & 353 & - gsinlon(:,:,ig) * gsinlat(:,:,ig) * pyy & 354 & + gcoslat(:,:,ig) * pzz 355 ! 356 END SUBROUTINE geo2oce 357 358 359 SUBROUTINE oce2geo ( pte, ptn, cgrid, pxx , pyy , pzz ) 360 !!---------------------------------------------------------------------- 361 !! *** ROUTINE oce2geo *** 362 !! 363 !! ** Purpose : 364 !! 365 !! ** Method : Change vector from east/north to geocentric 366 !! 367 !! History : ! (A. Caubel) oce2geo - Original code 368 !!---------------------------------------------------------------------- 369 REAL(wp), DIMENSION(jpi,jpj), INTENT( IN ) :: pte, ptn 370 CHARACTER(len=1) , INTENT( IN ) :: cgrid 371 REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT ) :: pxx , pyy , pzz 372 !! 373 REAL(wp), PARAMETER :: rpi = 3.141592653E0 374 REAL(wp), PARAMETER :: rad = rpi / 180.e0 375 INTEGER :: ig ! 376 INTEGER :: ierr ! local integer 377 !!---------------------------------------------------------------------- 378 379 IF( .NOT. ALLOCATED( gsinlon ) ) THEN 380 ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) , & 381 & gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr ) 382 IF( lk_mpp ) CALL mpp_sum( ierr ) 383 IF( ierr /= 0 ) CALL ctl_stop('oce2geo: unable to allocate arrays' ) 363 384 ENDIF 364 385 … … 404 425 CALL ctl_stop( ctmp1 ) 405 426 END SELECT 406 407 pte = - gsinlon(:,:,ig) * pxx + gcoslon(:,:,ig) * pyy 408 ptn = - gcoslon(:,:,ig) * gsinlat(:,:,ig) * pxx & 409 - gsinlon(:,:,ig) * gsinlat(:,:,ig) * pyy & 410 + gcoslat(:,:,ig) * pzz 411 !!$ ptv = gcoslon(:,:,ig) * gcoslat(:,:,ig) * pxx & 412 !!$ + gsinlon(:,:,ig) * gcoslat(:,:,ig) * pyy & 413 !!$ + gsinlat(:,:,ig) * pzz 414 ! 415 END SUBROUTINE geo2oce 416 417 SUBROUTINE oce2geo ( pte, ptn, cgrid, & 418 pxx , pyy , pzz ) 419 !!---------------------------------------------------------------------- 420 !! *** ROUTINE oce2geo *** 421 !! 422 !! ** Purpose : 423 !! 424 !! ** Method : Change vector from east/north to geocentric 425 !! 426 !! History : 427 !! ! (A. Caubel) oce2geo - Original code 428 !!---------------------------------------------------------------------- 429 REAL(wp), DIMENSION(jpi,jpj), INTENT( IN ) :: pte, ptn 430 CHARACTER(len=1) , INTENT( IN ) :: cgrid 431 REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT ) :: pxx , pyy , pzz 432 !! 433 REAL(wp), PARAMETER :: rpi = 3.141592653E0 434 REAL(wp), PARAMETER :: rad = rpi / 180.e0 435 INTEGER :: ig ! 436 INTEGER :: ierr ! local integer 437 !!---------------------------------------------------------------------- 438 439 IF( .NOT. ALLOCATED( gsinlon ) ) THEN 440 ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) , & 441 & gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr ) 442 IF( lk_mpp ) CALL mpp_sum( ierr ) 443 IF( ierr /= 0 ) CALL ctl_stop('oce2geo: unable to allocate arrays' ) 444 ENDIF 445 446 SELECT CASE( cgrid) 447 CASE ( 'T' ) 448 ig = 1 449 IF( .NOT. linit(ig) ) THEN 450 gsinlon(:,:,ig) = SIN( rad * glamt(:,:) ) 451 gcoslon(:,:,ig) = COS( rad * glamt(:,:) ) 452 gsinlat(:,:,ig) = SIN( rad * gphit(:,:) ) 453 gcoslat(:,:,ig) = COS( rad * gphit(:,:) ) 454 linit(ig) = .TRUE. 455 ENDIF 456 CASE ( 'U' ) 457 ig = 2 458 IF( .NOT. linit(ig) ) THEN 459 gsinlon(:,:,ig) = SIN( rad * glamu(:,:) ) 460 gcoslon(:,:,ig) = COS( rad * glamu(:,:) ) 461 gsinlat(:,:,ig) = SIN( rad * gphiu(:,:) ) 462 gcoslat(:,:,ig) = COS( rad * gphiu(:,:) ) 463 linit(ig) = .TRUE. 464 ENDIF 465 CASE ( 'V' ) 466 ig = 3 467 IF( .NOT. linit(ig) ) THEN 468 gsinlon(:,:,ig) = SIN( rad * glamv(:,:) ) 469 gcoslon(:,:,ig) = COS( rad * glamv(:,:) ) 470 gsinlat(:,:,ig) = SIN( rad * gphiv(:,:) ) 471 gcoslat(:,:,ig) = COS( rad * gphiv(:,:) ) 472 linit(ig) = .TRUE. 473 ENDIF 474 CASE ( 'F' ) 475 ig = 4 476 IF( .NOT. linit(ig) ) THEN 477 gsinlon(:,:,ig) = SIN( rad * glamf(:,:) ) 478 gcoslon(:,:,ig) = COS( rad * glamf(:,:) ) 479 gsinlat(:,:,ig) = SIN( rad * gphif(:,:) ) 480 gcoslat(:,:,ig) = COS( rad * gphif(:,:) ) 481 linit(ig) = .TRUE. 482 ENDIF 483 CASE default 484 WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid 485 CALL ctl_stop( ctmp1 ) 486 END SELECT 487 488 pxx = - gsinlon(:,:,ig) * pte - gcoslon(:,:,ig) * gsinlat(:,:,ig) * ptn 489 pyy = gcoslon(:,:,ig) * pte - gsinlon(:,:,ig) * gsinlat(:,:,ig) * ptn 490 pzz = gcoslat(:,:,ig) * ptn 491 492 427 ! 428 pxx = - gsinlon(:,:,ig) * pte - gcoslon(:,:,ig) * gsinlat(:,:,ig) * ptn 429 pyy = gcoslon(:,:,ig) * pte - gsinlon(:,:,ig) * gsinlat(:,:,ig) * ptn 430 pzz = gcoslat(:,:,ig) * ptn 431 ! 493 432 END SUBROUTINE oce2geo 494 433 495 434 496 SUBROUTINE repere ( px1, py1, px2, py2, kchoix, cd_type ) 497 !!---------------------------------------------------------------------- 498 !! *** ROUTINE repere *** 499 !! 500 !! ** Purpose : Change vector componantes between a geopgraphic grid 501 !! and a stretched coordinates grid. 502 !! 503 !! ** Method : 504 !! 505 !! ** Action : 506 !! 507 !! History : 508 !! ! 89-03 (O. Marti) original code 509 !! ! 92-02 (M. Imbard) 510 !! ! 93-03 (M. Guyon) symetrical conditions 511 !! ! 98-05 (B. Blanke) 512 !! 8.5 ! 02-08 (G. Madec) F90: Free form 513 !!---------------------------------------------------------------------- 514 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: px1, py1 ! two horizontal components to be rotated 515 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: px2, py2 ! the two horizontal components in the model repere 516 INTEGER , INTENT(in ) :: kchoix ! type of transformation 517 ! ! = 1 change from geographic to model grid. 518 ! ! =-1 change from model to geographic grid 519 CHARACTER(len=1), INTENT(in ), OPTIONAL :: cd_type ! define the nature of pt2d array grid-points 520 ! 521 CHARACTER(len=1) :: cl_type ! define the nature of pt2d array grid-points (T point by default) 522 !!---------------------------------------------------------------------- 523 524 cl_type = 'T' 525 IF( PRESENT(cd_type) ) cl_type = cd_type 526 ! 527 SELECT CASE (kchoix) 528 CASE ( 1) ! change from geographic to model grid. 529 CALL rot_rep( px1, py1, cl_type, 'en->i', px2 ) 530 CALL rot_rep( px1, py1, cl_type, 'en->j', py2 ) 531 CASE (-1) ! change from model to geographic grid 532 CALL rot_rep( px1, py1, cl_type, 'ij->e', px2 ) 533 CALL rot_rep( px1, py1, cl_type, 'ij->n', py2 ) 534 CASE DEFAULT ; CALL ctl_stop( 'repere: Syntax Error in the definition of kchoix (1 OR -1' ) 535 END SELECT 536 537 END SUBROUTINE repere 538 539 540 SUBROUTINE obs_rot ( psinu, pcosu, psinv, pcosv ) 435 SUBROUTINE obs_rot( psinu, pcosu, psinv, pcosv ) 541 436 !!---------------------------------------------------------------------- 542 437 !! *** ROUTINE obs_rot *** … … 546 441 !! current at observation points 547 442 !! 548 !! History : 549 !! 9.2 ! 09-02 (K. Mogensen) 443 !! History : 9.2 ! 09-02 (K. Mogensen) 550 444 !!---------------------------------------------------------------------- 551 445 REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT ):: psinu, pcosu, psinv, pcosv ! copy of data 552 446 !!---------------------------------------------------------------------- 553 447 ! 554 448 ! Initialization of gsin* and gcos* at first call 555 449 ! ----------------------------------------------- 556 557 450 IF( lmust_init ) THEN 558 451 IF(lwp) WRITE(numout,*) 559 452 IF(lwp) WRITE(numout,*) ' obs_rot : geographic <--> stretched' 560 453 IF(lwp) WRITE(numout,*) ' ~~~~~~~ coordinate transformation' 561 562 454 CALL angle ! initialization of the transformation 563 455 lmust_init = .FALSE. 564 565 456 ENDIF 566 457 ! 567 458 psinu(:,:) = gsinu(:,:) 568 459 pcosu(:,:) = gcosu(:,:) 569 460 psinv(:,:) = gsinv(:,:) 570 461 pcosv(:,:) = gcosv(:,:) 571 462 ! 572 463 END SUBROUTINE obs_rot 573 464 -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
r5836 r6140 44 44 LOGICAL , PUBLIC :: ln_dm2dc !: Daily mean to Diurnal Cycle short wave (qsr) 45 45 LOGICAL , PUBLIC :: ln_rnf !: runoffs / runoff mouths 46 LOGICAL , PUBLIC :: ln_isf !: ice shelf melting 46 47 LOGICAL , PUBLIC :: ln_ssr !: Sea Surface restoring on SST and/or SSS 47 48 LOGICAL , PUBLIC :: ln_apr_dyn !: Atmospheric pressure forcing used on dynamics (ocean & ice) 48 49 INTEGER , PUBLIC :: nn_ice !: flag for ice in the surface boundary condition (=0/1/2/3) 49 INTEGER , PUBLIC :: nn_isf !: flag for isf in the surface boundary condition (=0/1/2/3/4)50 50 INTEGER , PUBLIC :: nn_ice_embd !: flag for levitating/embedding sea-ice in the ocean 51 51 ! !: =0 levitating ice (no mass exchange, concentration/dilution effect) … … 172 172 & ssv_m (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) ) 173 173 ! 174 #if defined key_vvl175 174 ALLOCATE( e3t_m(jpi,jpj) , STAT=ierr(5) ) 176 #endif177 175 ! 178 176 sbc_oce_alloc = MAXVAL( ierr ) -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90
r5836 r6140 37 37 38 38 !! * Substitutions 39 # include "domzgr_substitute.h90"40 39 # include "vectopt_loop_substitute.h90" 41 40 !!---------------------------------------------------------------------- -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90
r5930 r6140 13 13 USE sbc_oce ! surface boundary condition 14 14 USE phycst ! physical constants 15 ! 15 16 USE fldread ! read input fields 16 17 USE in_out_manager ! I/O manager … … 38 39 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_apr ! structure of input fields (file informations, fields read) 39 40 40 !! * Substitutions41 # include "domzgr_substitute.h90"42 41 !!---------------------------------------------------------------------- 43 42 !! NEMO/OPA 4.0 , NEMO Consortium (2011) … … 112 111 ENDIF 113 112 !jc: stop below should rather be a warning 114 IF( ( ln_apr_obc ) .AND. .NOT.ln_apr_dyn ) &113 IF( ln_apr_obc .AND. .NOT.ln_apr_dyn ) & 115 114 CALL ctl_stop( 'sbc_apr: use inverse barometer ssh at open boundary ONLY requires ln_apr_dyn=T' ) 116 115 ENDIF -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
r5583 r6140 19 19 20 20 !!---------------------------------------------------------------------- 21 !! sbc_blk_core 22 !! blk_oce_core 23 !! blk_ice_core 24 !! turb_core_2z 25 !! cd_neutral_10m 26 !! psi_m 27 !! psi_h 21 !! sbc_blk_core : bulk formulation as ocean surface boundary condition (forced mode, CORE bulk formulea) 22 !! blk_oce_core : computes momentum, heat and freshwater fluxes over ocean 23 !! blk_ice_core : computes momentum, heat and freshwater fluxes over ice 24 !! turb_core_2z : Computes turbulent transfert coefficients 25 !! cd_neutral_10m: Estimate of the neutral drag coefficient at 10m 26 !! psi_m : universal profile stability function for momentum 27 !! psi_h : universal profile stability function for temperature and humidity 28 28 !!---------------------------------------------------------------------- 29 USE oce ! ocean dynamics and tracers 30 USE dom_oce ! ocean space and time domain 31 USE phycst ! physical constants 32 USE fldread ! read input fields 33 USE sbc_oce ! Surface boundary condition: ocean fields 34 USE cyclone ! Cyclone 10m wind form trac of cyclone centres 35 USE sbcdcy ! surface boundary condition: diurnal cycle 36 USE iom ! I/O manager library 37 USE in_out_manager ! I/O manager 38 USE lib_mpp ! distribued memory computing library 39 USE wrk_nemo ! work arrays 40 USE timing ! Timing 41 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 42 USE prtctl ! Print control 43 USE sbcwave, ONLY : cdn_wave ! wave module 44 USE sbc_ice ! Surface boundary condition: ice fields 45 USE lib_fortran ! to use key_nosignedzero 29 USE oce ! ocean dynamics and tracers 30 USE dom_oce ! ocean space and time domain 31 USE phycst ! physical constants 32 USE fldread ! read input fields 33 USE sbc_oce ! Surface boundary condition: ocean fields 34 USE cyclone ! Cyclone 10m wind form trac of cyclone centres 35 USE sbcdcy ! surface boundary condition: diurnal cycle 36 USE sbcwave , ONLY : cdn_wave ! wave module 37 USE sbc_ice ! Surface boundary condition: ice fields 38 USE lib_fortran ! to use key_nosignedzero 46 39 #if defined key_lim3 47 USE ice , ONLY :u_ice, v_ice, jpl, pfrld, a_i_b48 USE limthd_dh 40 USE ice , ONLY : u_ice, v_ice, jpl, pfrld, a_i_b 41 USE limthd_dh ! for CALL lim_thd_snwblow 49 42 #elif defined key_lim2 50 USE ice_2 , ONLY :u_ice, v_ice51 USE par_ice_2 43 USE ice_2 , ONLY : u_ice, v_ice 44 USE par_ice_2 ! LIM-2 parameters 52 45 #endif 46 ! 47 USE iom ! I/O manager library 48 USE in_out_manager ! I/O manager 49 USE lib_mpp ! distribued memory computing library 50 USE wrk_nemo ! work arrays 51 USE timing ! Timing 52 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 53 USE prtctl ! Print control 53 54 54 55 IMPLICIT NONE … … 84 85 REAL(wp), PARAMETER :: albo = 0.066 ! ocean albedo assumed to be constant 85 86 86 ! 87 ! !!* Namelist namsbc_core : CORE bulk parameters 87 88 LOGICAL :: ln_taudif ! logical flag to use the "mean of stress module - module of mean stress" data 88 89 REAL(wp) :: rn_pfac ! multiplication factor for precipitation … … 93 94 94 95 !! * Substitutions 95 # include "domzgr_substitute.h90"96 96 # include "vectopt_loop_substitute.h90" 97 97 !!---------------------------------------------------------------------- … … 149 149 TYPE(FLD_N) :: sn_tdif ! " " 150 150 NAMELIST/namsbc_core/ cn_dir , ln_taudif, rn_pfac, rn_efac, rn_vfac, & 151 & sn_wndi, sn_wndj , sn_humi, sn_qsr , &152 & sn_qlw , sn_tair , sn_prec, sn_snow, &153 & sn_tdif, rn_zqt , rn_zu151 & sn_wndi, sn_wndj , sn_humi, sn_qsr , & 152 & sn_qlw , sn_tair , sn_prec, sn_snow, & 153 & sn_tdif, rn_zqt , rn_zu 154 154 !!--------------------------------------------------------------------- 155 155 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90
r5836 r6140 8 8 9 9 !!---------------------------------------------------------------------- 10 !! sbc_blk_mfs : bulk formulation as ocean surface boundary condition10 !! sbc_blk_mfs : bulk formulation as ocean surface boundary condition 11 11 !! (forced mode, mfs bulk formulae) 12 !! blk_oce_mfs : ocean: computes momentum, heat and freshwater fluxes12 !! blk_oce_mfs : ocean: computes momentum, heat and freshwater fluxes 13 13 !!---------------------------------------------------------------------- 14 USE oce ! ocean dynamics and tracers 15 USE dom_oce ! ocean space and time domain 16 USE phycst ! physical constants 17 USE fldread ! read input fields 18 USE sbc_oce ! Surface boundary condition: ocean fields 19 USE iom ! I/O manager library 20 USE in_out_manager ! I/O manager 21 USE lib_mpp ! distribued memory computing library 22 USE wrk_nemo ! work arrays 23 USE timing ! Timing 24 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 25 USE prtctl ! Print control 26 USE sbcwave,ONLY : cdn_wave !wave module 14 USE oce ! ocean dynamics and tracers 15 USE dom_oce ! ocean space and time domain 16 USE phycst ! physical constants 17 USE fldread ! read input fields 18 USE sbc_oce ! Surface boundary condition: ocean fields 19 USE sbcwave ,ONLY : cdn_wave !wave module 20 ! 21 USE iom ! I/O manager library 22 USE in_out_manager ! I/O manager 23 USE lib_mpp ! distribued memory computing library 24 USE wrk_nemo ! work arrays 25 USE timing ! Timing 26 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 27 USE prtctl ! Print control 27 28 28 29 IMPLICIT NONE … … 42 43 43 44 !! * Substitutions 44 # include "domzgr_substitute.h90"45 45 # include "vectopt_loop_substitute.h90" 46 46 !!---------------------------------------------------------------------- … … 49 49 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 50 50 !!---------------------------------------------------------------------- 51 52 51 CONTAINS 53 54 52 55 53 SUBROUTINE sbc_blk_mfs( kt ) -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
r6121 r6140 18 18 !! sbc_cpl_snd : send fields to the atmosphere 19 19 !!---------------------------------------------------------------------- 20 USE dom_oce 21 USE sbc_oce 22 USE sbc_ice 23 USE sbcapr 24 USE sbcdcy 25 USE phycst 20 USE dom_oce ! ocean space and time domain 21 USE sbc_oce ! Surface boundary condition: ocean fields 22 USE sbc_ice ! Surface boundary condition: ice fields 23 USE sbcapr ! Stochastic param. : ??? 24 USE sbcdcy ! surface boundary condition: diurnal cycle 25 USE phycst ! physical constants 26 26 #if defined key_lim3 27 USE ice 27 USE ice ! ice variables 28 28 #endif 29 29 #if defined key_lim2 30 USE par_ice_2 31 USE ice_2 30 USE par_ice_2 ! ice parameters 31 USE ice_2 ! ice variables 32 32 #endif 33 USE cpl_oasis3 34 USE geo2ocean 33 USE cpl_oasis3 ! OASIS3 coupling 34 USE geo2ocean ! 35 35 USE oce , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev 36 USE albedo ! 37 USE in_out_manager ! I/O manager 38 USE iom ! NetCDF library 39 USE lib_mpp ! distribued memory computing library 40 USE wrk_nemo ! work arrays 41 USE timing ! Timing 42 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 43 USE eosbn2 44 USE sbcrnf , ONLY : l_rnfcpl 36 USE albedo ! 37 USE eosbn2 ! 38 USE sbcrnf , ONLY : l_rnfcpl 45 39 #if defined key_cpl_carbon_cycle 46 40 USE p4zflx, ONLY : oce_co2 … … 50 44 #endif 51 45 #if defined key_lim3 52 USE limthd_dh 46 USE limthd_dh ! for CALL lim_thd_snwblow 53 47 #endif 48 ! 49 USE in_out_manager ! I/O manager 50 USE iom ! NetCDF library 51 USE lib_mpp ! distribued memory computing library 52 USE wrk_nemo ! work arrays 53 USE timing ! Timing 54 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 54 55 55 56 IMPLICIT NONE 56 57 PRIVATE 57 58 58 PUBLIC sbc_cpl_init 59 PUBLIC sbc_cpl_rcv 60 PUBLIC sbc_cpl_snd 61 PUBLIC sbc_cpl_ice_tau 62 PUBLIC sbc_cpl_ice_flx 63 PUBLIC sbc_cpl_alloc 64 65 INTEGER, PARAMETER :: jpr_otx1 = 1 66 INTEGER, PARAMETER :: jpr_oty1 = 2 67 INTEGER, PARAMETER :: jpr_otz1 = 3 68 INTEGER, PARAMETER :: jpr_otx2 = 4 69 INTEGER, PARAMETER :: jpr_oty2 = 5 70 INTEGER, PARAMETER :: jpr_otz2 = 6 71 INTEGER, PARAMETER :: jpr_itx1 = 7 72 INTEGER, PARAMETER :: jpr_ity1 = 8 73 INTEGER, PARAMETER :: jpr_itz1 = 9 74 INTEGER, PARAMETER :: jpr_itx2 = 10 75 INTEGER, PARAMETER :: jpr_ity2 = 11 76 INTEGER, PARAMETER :: jpr_itz2 = 12 77 INTEGER, PARAMETER :: jpr_qsroce = 13 78 INTEGER, PARAMETER :: jpr_qsrice = 14 59 PUBLIC sbc_cpl_init ! routine called by sbcmod.F90 60 PUBLIC sbc_cpl_rcv ! routine called by sbc_ice_lim(_2).F90 61 PUBLIC sbc_cpl_snd ! routine called by step.F90 62 PUBLIC sbc_cpl_ice_tau ! routine called by sbc_ice_lim(_2).F90 63 PUBLIC sbc_cpl_ice_flx ! routine called by sbc_ice_lim(_2).F90 64 PUBLIC sbc_cpl_alloc ! routine called in sbcice_cice.F90 65 66 INTEGER, PARAMETER :: jpr_otx1 = 1 ! 3 atmosphere-ocean stress components on grid 1 67 INTEGER, PARAMETER :: jpr_oty1 = 2 ! 68 INTEGER, PARAMETER :: jpr_otz1 = 3 ! 69 INTEGER, PARAMETER :: jpr_otx2 = 4 ! 3 atmosphere-ocean stress components on grid 2 70 INTEGER, PARAMETER :: jpr_oty2 = 5 ! 71 INTEGER, PARAMETER :: jpr_otz2 = 6 ! 72 INTEGER, PARAMETER :: jpr_itx1 = 7 ! 3 atmosphere-ice stress components on grid 1 73 INTEGER, PARAMETER :: jpr_ity1 = 8 ! 74 INTEGER, PARAMETER :: jpr_itz1 = 9 ! 75 INTEGER, PARAMETER :: jpr_itx2 = 10 ! 3 atmosphere-ice stress components on grid 2 76 INTEGER, PARAMETER :: jpr_ity2 = 11 ! 77 INTEGER, PARAMETER :: jpr_itz2 = 12 ! 78 INTEGER, PARAMETER :: jpr_qsroce = 13 ! Qsr above the ocean 79 INTEGER, PARAMETER :: jpr_qsrice = 14 ! Qsr above the ice 79 80 INTEGER, PARAMETER :: jpr_qsrmix = 15 80 INTEGER, PARAMETER :: jpr_qnsoce = 16 81 INTEGER, PARAMETER :: jpr_qnsice = 17 81 INTEGER, PARAMETER :: jpr_qnsoce = 16 ! Qns above the ocean 82 INTEGER, PARAMETER :: jpr_qnsice = 17 ! Qns above the ice 82 83 INTEGER, PARAMETER :: jpr_qnsmix = 18 83 INTEGER, PARAMETER :: jpr_rain = 19 84 INTEGER, PARAMETER :: jpr_snow = 20 85 INTEGER, PARAMETER :: jpr_tevp = 21 86 INTEGER, PARAMETER :: jpr_ievp = 22 87 INTEGER, PARAMETER :: jpr_sbpr = 23 88 INTEGER, PARAMETER :: jpr_semp = 24 89 INTEGER, PARAMETER :: jpr_oemp = 25 90 INTEGER, PARAMETER :: jpr_w10m = 26 91 INTEGER, PARAMETER :: jpr_dqnsdt = 27 92 INTEGER, PARAMETER :: jpr_rnf = 28 93 INTEGER, PARAMETER :: jpr_cal = 29 94 INTEGER, PARAMETER :: jpr_taum = 30 84 INTEGER, PARAMETER :: jpr_rain = 19 ! total liquid precipitation (rain) 85 INTEGER, PARAMETER :: jpr_snow = 20 ! solid precipitation over the ocean (snow) 86 INTEGER, PARAMETER :: jpr_tevp = 21 ! total evaporation 87 INTEGER, PARAMETER :: jpr_ievp = 22 ! solid evaporation (sublimation) 88 INTEGER, PARAMETER :: jpr_sbpr = 23 ! sublimation - liquid precipitation - solid precipitation 89 INTEGER, PARAMETER :: jpr_semp = 24 ! solid freshwater budget (sublimation - snow) 90 INTEGER, PARAMETER :: jpr_oemp = 25 ! ocean freshwater budget (evap - precip) 91 INTEGER, PARAMETER :: jpr_w10m = 26 ! 10m wind 92 INTEGER, PARAMETER :: jpr_dqnsdt = 27 ! d(Q non solar)/d(temperature) 93 INTEGER, PARAMETER :: jpr_rnf = 28 ! runoffs 94 INTEGER, PARAMETER :: jpr_cal = 29 ! calving 95 INTEGER, PARAMETER :: jpr_taum = 30 ! wind stress module 95 96 INTEGER, PARAMETER :: jpr_co2 = 31 96 INTEGER, PARAMETER :: jpr_topm = 32 97 INTEGER, PARAMETER :: jpr_botm = 33 98 INTEGER, PARAMETER :: jpr_sflx = 34 99 INTEGER, PARAMETER :: jpr_toce = 35 100 INTEGER, PARAMETER :: jpr_soce = 36 101 INTEGER, PARAMETER :: jpr_ocx1 = 37 102 INTEGER, PARAMETER :: jpr_ocy1 = 38 103 INTEGER, PARAMETER :: jpr_ssh = 39 104 INTEGER, PARAMETER :: jpr_fice = 40 105 INTEGER, PARAMETER :: jpr_e3t1st = 41 106 INTEGER, PARAMETER :: jpr_fraqsr = 42 107 INTEGER, PARAMETER :: jprcv = 42 108 109 INTEGER, PARAMETER :: jps_fice = 1 110 INTEGER, PARAMETER :: jps_toce = 2 111 INTEGER, PARAMETER :: jps_tice = 3 112 INTEGER, PARAMETER :: jps_tmix = 4 113 INTEGER, PARAMETER :: jps_albice = 5 114 INTEGER, PARAMETER :: jps_albmix = 6 115 INTEGER, PARAMETER :: jps_hice = 7 116 INTEGER, PARAMETER :: jps_hsnw = 8 117 INTEGER, PARAMETER :: jps_ocx1 = 9 118 INTEGER, PARAMETER :: jps_ocy1 = 10 119 INTEGER, PARAMETER :: jps_ocz1 = 11 120 INTEGER, PARAMETER :: jps_ivx1 = 12 121 INTEGER, PARAMETER :: jps_ivy1 = 13 122 INTEGER, PARAMETER :: jps_ivz1 = 14 97 INTEGER, PARAMETER :: jpr_topm = 32 ! topmeltn 98 INTEGER, PARAMETER :: jpr_botm = 33 ! botmeltn 99 INTEGER, PARAMETER :: jpr_sflx = 34 ! salt flux 100 INTEGER, PARAMETER :: jpr_toce = 35 ! ocean temperature 101 INTEGER, PARAMETER :: jpr_soce = 36 ! ocean salinity 102 INTEGER, PARAMETER :: jpr_ocx1 = 37 ! ocean current on grid 1 103 INTEGER, PARAMETER :: jpr_ocy1 = 38 ! 104 INTEGER, PARAMETER :: jpr_ssh = 39 ! sea surface height 105 INTEGER, PARAMETER :: jpr_fice = 40 ! ice fraction 106 INTEGER, PARAMETER :: jpr_e3t1st = 41 ! first T level thickness 107 INTEGER, PARAMETER :: jpr_fraqsr = 42 ! fraction of solar net radiation absorbed in the first ocean level 108 INTEGER, PARAMETER :: jprcv = 42 ! total number of fields received 109 110 INTEGER, PARAMETER :: jps_fice = 1 ! ice fraction sent to the atmosphere 111 INTEGER, PARAMETER :: jps_toce = 2 ! ocean temperature 112 INTEGER, PARAMETER :: jps_tice = 3 ! ice temperature 113 INTEGER, PARAMETER :: jps_tmix = 4 ! mixed temperature (ocean+ice) 114 INTEGER, PARAMETER :: jps_albice = 5 ! ice albedo 115 INTEGER, PARAMETER :: jps_albmix = 6 ! mixed albedo 116 INTEGER, PARAMETER :: jps_hice = 7 ! ice thickness 117 INTEGER, PARAMETER :: jps_hsnw = 8 ! snow thickness 118 INTEGER, PARAMETER :: jps_ocx1 = 9 ! ocean current on grid 1 119 INTEGER, PARAMETER :: jps_ocy1 = 10 ! 120 INTEGER, PARAMETER :: jps_ocz1 = 11 ! 121 INTEGER, PARAMETER :: jps_ivx1 = 12 ! ice current on grid 1 122 INTEGER, PARAMETER :: jps_ivy1 = 13 ! 123 INTEGER, PARAMETER :: jps_ivz1 = 14 ! 123 124 INTEGER, PARAMETER :: jps_co2 = 15 124 INTEGER, PARAMETER :: jps_soce = 16 125 INTEGER, PARAMETER :: jps_ssh = 17 126 INTEGER, PARAMETER :: jps_qsroce = 18 127 INTEGER, PARAMETER :: jps_qnsoce = 19 128 INTEGER, PARAMETER :: jps_oemp = 20 129 INTEGER, PARAMETER :: jps_sflx = 21 130 INTEGER, PARAMETER :: jps_otx1 = 22 131 INTEGER, PARAMETER :: jps_oty1 = 23 132 INTEGER, PARAMETER :: jps_rnf = 24 133 INTEGER, PARAMETER :: jps_taum = 25 134 INTEGER, PARAMETER :: jps_fice2 = 26 135 INTEGER, PARAMETER :: jps_e3t1st = 27 136 INTEGER, PARAMETER :: jps_fraqsr = 28 137 INTEGER, PARAMETER :: jpsnd = 28 138 139 ! 140 TYPE :: FLD_C 141 CHARACTER(len = 32) :: cldes 142 CHARACTER(len = 32) :: clcat 143 CHARACTER(len = 32) :: clvref 144 CHARACTER(len = 32) :: clvor 145 CHARACTER(len = 32) :: clvgrd 125 INTEGER, PARAMETER :: jps_soce = 16 ! ocean salinity 126 INTEGER, PARAMETER :: jps_ssh = 17 ! sea surface height 127 INTEGER, PARAMETER :: jps_qsroce = 18 ! Qsr above the ocean 128 INTEGER, PARAMETER :: jps_qnsoce = 19 ! Qns above the ocean 129 INTEGER, PARAMETER :: jps_oemp = 20 ! ocean freshwater budget (evap - precip) 130 INTEGER, PARAMETER :: jps_sflx = 21 ! salt flux 131 INTEGER, PARAMETER :: jps_otx1 = 22 ! 2 atmosphere-ocean stress components on grid 1 132 INTEGER, PARAMETER :: jps_oty1 = 23 ! 133 INTEGER, PARAMETER :: jps_rnf = 24 ! runoffs 134 INTEGER, PARAMETER :: jps_taum = 25 ! wind stress module 135 INTEGER, PARAMETER :: jps_fice2 = 26 ! ice fraction sent to OPA (by SAS when doing SAS-OPA coupling) 136 INTEGER, PARAMETER :: jps_e3t1st = 27 ! first level depth (vvl) 137 INTEGER, PARAMETER :: jps_fraqsr = 28 ! fraction of solar net radiation absorbed in the first ocean level 138 INTEGER, PARAMETER :: jpsnd = 28 ! total number of fields sended 139 140 ! !!** namelist namsbc_cpl ** 141 TYPE :: FLD_C ! 142 CHARACTER(len = 32) :: cldes ! desciption of the coupling strategy 143 CHARACTER(len = 32) :: clcat ! multiple ice categories strategy 144 CHARACTER(len = 32) :: clvref ! reference of vector ('spherical' or 'cartesian') 145 CHARACTER(len = 32) :: clvor ! orientation of vector fields ('eastward-northward' or 'local grid') 146 CHARACTER(len = 32) :: clvgrd ! grids on which is located the vector fields 146 147 END TYPE FLD_C 147 ! Send to the atmosphere !148 ! ! Send to the atmosphere 148 149 TYPE(FLD_C) :: sn_snd_temp, sn_snd_alb, sn_snd_thick, sn_snd_crt, sn_snd_co2 149 ! Received from the atmosphere !150 ! ! Received from the atmosphere 150 151 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 152 TYPE(FLD_C) :: sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2 152 ! Other namelist parameters !153 INTEGER :: nn_cplmodel 154 LOGICAL :: ln_usecplmask 155 153 ! ! Other namelist parameters 154 INTEGER :: nn_cplmodel ! Maximum number of models to/from which NEMO is potentialy sending/receiving data 155 LOGICAL :: ln_usecplmask ! use a coupling mask file to merge data received from several models 156 ! -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel) 156 157 TYPE :: DYNARR 157 158 REAL(wp), POINTER, DIMENSION(:,:,:) :: z3 158 159 END TYPE DYNARR 159 160 160 TYPE( DYNARR ), SAVE, DIMENSION(jprcv) :: frcv 161 162 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: albedo_oce_mix 163 164 INTEGER , ALLOCATABLE, SAVE, DIMENSION( :) :: nrcvinfo 161 TYPE( DYNARR ), SAVE, DIMENSION(jprcv) :: frcv ! all fields recieved from the atmosphere 162 163 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: albedo_oce_mix ! ocean albedo sent to atmosphere (mix clear/overcast sky) 164 165 INTEGER , ALLOCATABLE, SAVE, DIMENSION( :) :: nrcvinfo ! OASIS info argument 165 166 166 167 !! Substitution 167 # include "domzgr_substitute.h90"168 168 # include "vectopt_loop_substitute.h90" 169 169 !!---------------------------------------------------------------------- 170 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)170 !! NEMO/OPA 3.7 , NEMO Consortium (2015) 171 171 !! $Id$ 172 172 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 173 173 !!---------------------------------------------------------------------- 174 175 174 CONTAINS 176 175 … … 209 208 !! * initialise the OASIS coupler 210 209 !!---------------------------------------------------------------------- 211 INTEGER, INTENT(in) :: k_ice ! ice management in the sbc (=0/1/2/3) 212 !! 213 INTEGER :: jn ! dummy loop index 214 INTEGER :: ios ! Local integer output status for namelist read 215 INTEGER :: inum 210 INTEGER, INTENT(in) :: k_ice ! ice management in the sbc (=0/1/2/3) 211 ! 212 INTEGER :: jn ! dummy loop index 213 INTEGER :: ios, inum ! Local integer 216 214 REAL(wp), POINTER, DIMENSION(:,:) :: zacs, zaos 217 215 !! … … 222 220 !!--------------------------------------------------------------------- 223 221 ! 224 IF( nn_timing == 1 ) CALL timing_start('sbc_cpl_init')225 ! 226 CALL wrk_alloc( jpi,jpj, zacs, zaos )222 IF( nn_timing == 1 ) CALL timing_start('sbc_cpl_init') 223 ! 224 CALL wrk_alloc( jpi,jpj, zacs, zaos ) 227 225 228 226 ! ================================ ! 229 227 ! Namelist informations ! 230 228 ! ================================ ! 231 229 ! 232 230 REWIND( numnam_ref ) ! Namelist namsbc_cpl in reference namelist : Variables for OASIS coupling 233 231 READ ( numnam_ref, namsbc_cpl, IOSTAT = ios, ERR = 901) 234 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cpl in reference namelist', lwp )235 232 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cpl in reference namelist', lwp ) 233 ! 236 234 REWIND( numnam_cfg ) ! Namelist namsbc_cpl in configuration namelist : Variables for OASIS coupling 237 235 READ ( numnam_cfg, namsbc_cpl, IOSTAT = ios, ERR = 902 ) 238 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cpl in configuration namelist', lwp )236 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cpl in configuration namelist', lwp ) 239 237 IF(lwm) WRITE ( numond, namsbc_cpl ) 240 238 ! 241 239 IF(lwp) THEN ! control print 242 240 WRITE(numout,*) … … 374 372 srcv(jpr_ity1)%clgrid = 'V' ! i.e. it is always at U- & V-points for i- & j-comp. resp. 375 373 ENDIF 376 374 ! 377 375 ! ! ------------------------- ! 378 376 ! ! freshwater budget ! E-P … … 396 394 CASE default ; CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_emp%cldes' ) 397 395 END SELECT 398 396 ! 399 397 ! ! ------------------------- ! 400 398 ! ! Runoffs & Calving ! … … 410 408 ! 411 409 srcv(jpr_cal )%clname = 'OCalving' ; IF( TRIM( sn_rcv_cal%cldes ) == 'coupled' ) srcv(jpr_cal)%laction = .TRUE. 412 410 ! 413 411 ! ! ------------------------- ! 414 412 ! ! non solar radiation ! Qns … … 535 533 IF( .NOT. ln_cpl ) srcv(:)%nsgn = 1. ! force default definition in case of opa <-> sas coupling 536 534 srcv( (/jpr_toce, jpr_soce, jpr_ssh, jpr_fraqsr, jpr_ocx1, jpr_ocy1/) )%laction = .TRUE. 537 srcv( jpr_e3t1st )%laction = lk_vvl535 srcv( jpr_e3t1st )%laction = .NOT.ln_linssh 538 536 srcv(jpr_ocx1)%clgrid = 'U' ! oce components given at U-point 539 537 srcv(jpr_ocy1)%clgrid = 'V' ! and V-point … … 701 699 ssnd(:)%laction = .FALSE. ! force default definition in case of opa <-> sas coupling 702 700 ssnd( (/jps_toce, jps_soce, jps_ssh, jps_fraqsr, jps_ocx1, jps_ocy1/) )%laction = .TRUE. 703 ssnd( jps_e3t1st )%laction = lk_vvl701 ssnd( jps_e3t1st )%laction = .NOT.ln_linssh 704 702 ! vector definition: not used but cleaner... 705 703 ssnd(jps_ocx1)%clgrid = 'U' ! oce components given at U-point … … 785 783 ncpl_qsr_freq = 86400 / ncpl_qsr_freq 786 784 787 CALL wrk_dealloc( jpi,jpj, zacs, zaos )788 ! 789 IF( nn_timing == 1 ) CALL timing_stop('sbc_cpl_init')785 CALL wrk_dealloc( jpi,jpj, zacs, zaos ) 786 ! 787 IF( nn_timing == 1 ) CALL timing_stop('sbc_cpl_init') 790 788 ! 791 789 END SUBROUTINE sbc_cpl_init … … 837 835 !! emp upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case) 838 836 !!---------------------------------------------------------------------- 839 INTEGER, INTENT(in) :: kt! ocean model time step index840 INTEGER, INTENT(in) :: k_fsbc! frequency of sbc (-> ice model) computation841 INTEGER, INTENT(in) :: k_ice! ice management in the sbc (=0/1/2/3)837 INTEGER, INTENT(in) :: kt ! ocean model time step index 838 INTEGER, INTENT(in) :: k_fsbc ! frequency of sbc (-> ice model) computation 839 INTEGER, INTENT(in) :: k_ice ! ice management in the sbc (=0/1/2/3) 842 840 843 841 !! 844 842 LOGICAL :: llnewtx, llnewtau ! update wind stress components and module?? 845 843 INTEGER :: ji, jj, jn ! dummy loop indices 846 INTEGER :: isec ! number of seconds since nit000 (assuming rdt tradid not change since nit000)844 INTEGER :: isec ! number of seconds since nit000 (assuming rdt did not change since nit000) 847 845 REAL(wp) :: zcumulneg, zcumulpos ! temporary scalars 848 846 REAL(wp) :: zcoef ! temporary scalar … … 853 851 !!---------------------------------------------------------------------- 854 852 ! 855 IF( nn_timing == 1 ) CALL timing_start('sbc_cpl_rcv')856 ! 857 CALL wrk_alloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr )853 IF( nn_timing == 1 ) CALL timing_start('sbc_cpl_rcv') 854 ! 855 CALL wrk_alloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr ) 858 856 ! 859 857 IF( ln_mixcpl ) zmsk(:,:) = 1. - xcplmask(:,:,0) … … 862 860 ! ! Receive all the atmos. fields (including ice information) 863 861 ! ! ======================================================= ! 864 isec = ( kt - nit000 ) * NINT( rdt tra(1) )! date of exchanges862 isec = ( kt - nit000 ) * NINT( rdt ) ! date of exchanges 865 863 DO jn = 1, jprcv ! received fields sent by the atmosphere 866 864 IF( srcv(jn)%laction ) CALL cpl_rcv( jn, isec, frcv(jn)%z3, xcplmask(:,:,1:nn_cplmodel), nrcvinfo(jn) ) … … 1104 1102 IF( srcv(jpr_fice )%laction ) fr_i(:,:) = frcv(jpr_fice )%z3(:,:,1) 1105 1103 ! 1106 1107 ENDIF 1108 ! 1109 CALL wrk_dealloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr ) 1110 ! 1111 IF( nn_timing == 1 ) CALL timing_stop('sbc_cpl_rcv') 1104 ENDIF 1105 ! 1106 CALL wrk_dealloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr ) 1107 ! 1108 IF( nn_timing == 1 ) CALL timing_stop('sbc_cpl_rcv') 1112 1109 ! 1113 1110 END SUBROUTINE sbc_cpl_rcv … … 1150 1147 REAL(wp), INTENT(out), DIMENSION(:,:) :: p_tauj ! at I-point (B-grid) or U & V-point (C-grid) 1151 1148 !! 1152 INTEGER :: ji, jj 1153 INTEGER :: itx 1149 INTEGER :: ji, jj ! dummy loop indices 1150 INTEGER :: itx ! index of taux over ice 1154 1151 REAL(wp), POINTER, DIMENSION(:,:) :: ztx, zty 1155 1152 !!---------------------------------------------------------------------- 1156 1153 ! 1157 IF( nn_timing == 1 ) CALL timing_start('sbc_cpl_ice_tau')1158 ! 1159 CALL wrk_alloc( jpi,jpj, ztx, zty )1154 IF( nn_timing == 1 ) CALL timing_start('sbc_cpl_ice_tau') 1155 ! 1156 CALL wrk_alloc( jpi,jpj, ztx, zty ) 1160 1157 1161 1158 IF( srcv(jpr_itx1)%laction ) THEN ; itx = jpr_itx1 … … 1165 1162 ! do something only if we just received the stress from atmosphere 1166 1163 IF( nrcvinfo(itx) == OASIS_Rcv ) THEN 1167 1168 1164 ! ! ======================= ! 1169 1165 IF( srcv(jpr_itx1)%laction ) THEN ! ice stress received ! … … 1318 1314 ENDIF 1319 1315 ! 1320 CALL wrk_dealloc( jpi,jpj, ztx, zty )1321 ! 1322 IF( nn_timing == 1 ) CALL timing_stop('sbc_cpl_ice_tau')1316 CALL wrk_dealloc( jpi,jpj, ztx, zty ) 1317 ! 1318 IF( nn_timing == 1 ) CALL timing_stop('sbc_cpl_ice_tau') 1323 1319 ! 1324 1320 END SUBROUTINE sbc_cpl_ice_tau … … 1365 1361 !! sprecip solid precipitation over the ocean 1366 1362 !!---------------------------------------------------------------------- 1367 REAL(wp), INTENT(in ), DIMENSION(:,:) :: p_frld ! lead fraction[0 to 1]1363 REAL(wp), INTENT(in ), DIMENSION(:,:) :: p_frld ! lead fraction [0 to 1] 1368 1364 ! optional arguments, used only in 'mixed oce-ice' case 1369 REAL(wp), INTENT(in ), DIMENSION(:,:,:), OPTIONAL :: palbi 1370 REAL(wp), INTENT(in ), DIMENSION(:,: ), OPTIONAL :: psst ! sea surface temperature[Celsius]1371 REAL(wp), INTENT(in ), DIMENSION(:,:,:), OPTIONAL :: pist ! ice surface temperature[Kelvin]1372 ! 1373 INTEGER :: jl 1365 REAL(wp), INTENT(in ), DIMENSION(:,:,:), OPTIONAL :: palbi ! all skies ice albedo 1366 REAL(wp), INTENT(in ), DIMENSION(:,: ), OPTIONAL :: psst ! sea surface temperature [Celsius] 1367 REAL(wp), INTENT(in ), DIMENSION(:,:,:), OPTIONAL :: pist ! ice surface temperature [Kelvin] 1368 ! 1369 INTEGER :: jl ! dummy loop index 1374 1370 REAL(wp), POINTER, DIMENSION(:,: ) :: zcptn, ztmp, zicefr, zmsk 1375 1371 REAL(wp), POINTER, DIMENSION(:,: ) :: zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot … … 1378 1374 !!---------------------------------------------------------------------- 1379 1375 ! 1380 IF( nn_timing == 1 ) CALL timing_start('sbc_cpl_ice_flx')1381 ! 1382 CALL wrk_alloc( jpi,jpj, zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot )1383 CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice )1376 IF( nn_timing == 1 ) CALL timing_start('sbc_cpl_ice_flx') 1377 ! 1378 CALL wrk_alloc( jpi,jpj, zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot ) 1379 CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice ) 1384 1380 1385 1381 IF( ln_mixcpl ) zmsk(:,:) = 1. - xcplmask(:,:,0) … … 1554 1550 CALL wrk_dealloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce ) 1555 1551 #else 1556 1557 ! clem: this formulation is certainly wrong... but better than it was ...1552 ! 1553 ! clem: this formulation is certainly wrong... but better than it was before... 1558 1554 zqns_tot(:,:) = zqns_tot(:,:) & ! zqns_tot update over free ocean with: 1559 1555 & - ztmp(:,:) & ! remove the latent heat flux of solid precip. melting … … 1571 1567 qns_ice(:,:,:) = zqns_ice(:,:,:) 1572 1568 ENDIF 1573 1569 ! 1574 1570 #endif 1575 1576 1571 ! ! ========================= ! 1577 1572 SELECT CASE( TRIM( sn_rcv_qsr%cldes ) ) ! solar heat fluxes ! (qsr) … … 1682 1677 fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 1683 1678 1684 CALL wrk_dealloc( jpi,jpj, zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot )1685 CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice )1686 ! 1687 IF( nn_timing == 1 ) CALL timing_stop('sbc_cpl_ice_flx')1679 CALL wrk_dealloc( jpi,jpj, zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot ) 1680 CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice ) 1681 ! 1682 IF( nn_timing == 1 ) CALL timing_stop('sbc_cpl_ice_flx') 1688 1683 ! 1689 1684 END SUBROUTINE sbc_cpl_ice_flx … … 1708 1703 !!---------------------------------------------------------------------- 1709 1704 ! 1710 IF( nn_timing == 1 ) CALL timing_start('sbc_cpl_snd')1711 ! 1712 CALL wrk_alloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )1713 CALL wrk_alloc( jpi,jpj,jpl, ztmp3, ztmp4 )1714 1715 isec = ( kt - nit000 ) * NINT( rdttra(1)) ! date of exchanges1705 IF( nn_timing == 1 ) CALL timing_start('sbc_cpl_snd') 1706 ! 1707 CALL wrk_alloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 ) 1708 CALL wrk_alloc( jpi,jpj,jpl, ztmp3, ztmp4 ) 1709 1710 isec = ( kt - nit000 ) * NINT( rdt ) ! date of exchanges 1716 1711 1717 1712 zfr_l(:,:) = 1.- fr_i(:,:) … … 2031 2026 ! ! first T level thickness 2032 2027 IF( ssnd(jps_e3t1st )%laction ) THEN 2033 CALL cpl_snd( jps_e3t1st, isec, RESHAPE ( fse3t_n(:,:,1) , (/jpi,jpj,1/) ), info )2028 CALL cpl_snd( jps_e3t1st, isec, RESHAPE ( e3t_n(:,:,1) , (/jpi,jpj,1/) ), info ) 2034 2029 ENDIF 2035 2030 ! ! Qsr fraction … … 2049 2044 IF( ssnd(jps_taum )%laction ) CALL cpl_snd( jps_taum , isec, RESHAPE ( taum, (/jpi,jpj,1/) ), info ) 2050 2045 2051 CALL wrk_dealloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )2052 CALL wrk_dealloc( jpi,jpj,jpl, ztmp3, ztmp4 )2053 ! 2054 IF( nn_timing == 1 ) CALL timing_stop('sbc_cpl_snd')2046 CALL wrk_dealloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 ) 2047 CALL wrk_dealloc( jpi,jpj,jpl, ztmp3, ztmp4 ) 2048 ! 2049 IF( nn_timing == 1 ) CALL timing_stop('sbc_cpl_snd') 2055 2050 ! 2056 2051 END SUBROUTINE sbc_cpl_snd -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcdcy.F90
r3764 r6140 90 90 91 91 ! When are we during the day (from 0 to 1) 92 zlo = ( REAL(nsec_day, wp) - 0.5_wp * rdt tra(1)) / rday93 zup = zlo + ( REAL(nn_fsbc, wp) * rdt tra(1)) / rday92 zlo = ( REAL(nsec_day, wp) - 0.5_wp * rdt ) / rday 93 zup = zlo + ( REAL(nn_fsbc, wp) * rdt ) / rday 94 94 ! 95 95 IF( nday_qsr == -1 ) THEN ! first time step only … … 189 189 END DO 190 190 ! 191 ztmp = rday / ( rdt tra(1)* REAL(nn_fsbc, wp) )191 ztmp = rday / ( rdt * REAL(nn_fsbc, wp) ) 192 192 rscal(:,:) = rscal(:,:) * ztmp 193 193 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90
r5836 r6140 17 17 USE sbcdcy ! surface boundary condition: diurnal cycle on qsr 18 18 USE phycst ! physical constants 19 ! 19 20 USE fldread ! read input fields 20 21 USE iom ! IOM library … … 37 38 38 39 !! * Substitutions 39 # include "domzgr_substitute.h90"40 40 # include "vectopt_loop_substitute.h90" 41 41 !!---------------------------------------------------------------------- … … 165 165 WRITE(numout,*) 166 166 WRITE(numout,*) ' day: ', ndastp , TRIM(sf(jf)%clvar), ' * ', zfact 167 CALL prihre( sf(jf)%fnow, jpi, jpj, 1, jpi, 20, 1, jpj, 10, zfact, numout )168 167 END DO 169 CALL FLUSH(numout)170 168 ENDIF 171 169 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90
r5643 r6140 12 12 13 13 !!---------------------------------------------------------------------- 14 !! sbc_fwb : freshwater budget for global ocean configurations15 !! in free surface and forced mode16 !!----------------------------------------------------------------------17 USE oce ! ocean dynamics and tracers18 USE dom_oce ! ocean space and time domain19 USE sbc_oce ! surface ocean boundary condition20 USE phycst ! physical constants21 USE sbc rnf ! ocean runoffs22 USE sbc isf ! ice shelf melting contribution23 USE sbcssr ! SS damping terms24 USE in_out_manager 25 USE lib_mpp 26 USE wrk_nemo 27 USE timing 28 USE lbclnk 29 USE lib_fortran 14 !! sbc_fwb : freshwater budget for global ocean configurations (free surface & forced mode) 15 !!---------------------------------------------------------------------- 16 USE oce ! ocean dynamics and tracers 17 USE dom_oce ! ocean space and time domain 18 USE sbc_oce ! surface ocean boundary condition 19 USE phycst ! physical constants 20 USE sbcrnf ! ocean runoffs 21 USE sbcisf ! ice shelf melting contribution 22 USE sbcssr ! Sea-Surface damping terms 23 ! 24 USE in_out_manager ! I/O manager 25 USE lib_mpp ! distribued memory computing library 26 USE wrk_nemo ! work arrays 27 USE timing ! Timing 28 USE lbclnk ! ocean lateral boundary conditions 29 USE lib_fortran ! 30 30 31 31 IMPLICIT NONE … … 40 40 41 41 !! * Substitutions 42 # include "domzgr_substitute.h90"43 42 # include "vectopt_loop_substitute.h90" 44 43 !!---------------------------------------------------------------------- … … 129 128 ENDIF 130 129 ! ! Update fwfold if new year start 131 ikty = 365 * 86400 / rdt tra(1)!!bug use of 365 days leap year or 360d year !!!!!!!130 ikty = 365 * 86400 / rdt !!bug use of 365 days leap year or 360d year !!!!!!! 132 131 IF( MOD( kt, ikty ) == 0 ) THEN 133 132 a_fwb_b = a_fwb ! mean sea level taking into account the ice+snow -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90
r5836 r6140 9 9 !!---------------------------------------------------------------------- 10 10 !! sbc_ice_cice : sea-ice model time-stepping and update ocean sbc over ice-covered area 11 !!12 !!13 11 !!---------------------------------------------------------------------- 14 12 USE oce ! ocean dynamics and tracers … … 92 90 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:), PRIVATE :: png ! local array used in sbc_cice_ice 93 91 94 !! * Substitutions95 # include "domzgr_substitute.h90"96 92 !!---------------------------------------------------------------------- 97 93 !! NEMO/OPA 3.7 , NEMO-consortium (2015) … … 244 240 snwice_mass_b(:,:) = 0.0_wp ! no mass exchanges 245 241 ENDIF 246 IF( .NOT. 242 IF( .NOT.ln_rstart ) THEN 247 243 IF( nn_ice_embd == 2 ) THEN ! full embedment (case 2) deplete the initial ssh below sea-ice area 248 244 sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0 249 245 sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0 250 #if defined key_vvl 251 ! key_vvl necessary? clem: yes for compilation purpose 252 DO jk = 1,jpkm1 ! adjust initial vertical scale factors 253 fse3t_n(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 254 fse3t_b(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 255 ENDDO 256 fse3t_a(:,:,:) = fse3t_b(:,:,:) 257 ! Reconstruction of all vertical scale factors at now and before time 258 ! steps 259 ! ============================================================================= 260 ! Horizontal scale factor interpolations 261 ! -------------------------------------- 262 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 263 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 264 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 265 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 266 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' ) 267 ! Vertical scale factor interpolations 268 ! ------------------------------------ 269 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 270 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 271 CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 272 CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 273 CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 274 ! t- and w- points depth 275 ! ---------------------- 276 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 277 fsdepw_n(:,:,1) = 0.0_wp 278 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 279 DO jk = 2, jpk 280 fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 281 fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 282 fsde3w_n(:,:,jk) = fsdept_n(:,:,jk ) - sshn (:,:) 283 END DO 284 #endif 246 247 !!gm This should be put elsewhere.... (same remark for limsbc) 248 !!gm especially here it is assumed zstar coordinate, but it can be ztilde.... 249 IF( .NOT.ln_linssh ) THEN 250 ! 251 DO jk = 1,jpkm1 ! adjust initial vertical scale factors 252 e3t_n(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 253 e3t_b(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 254 ENDDO 255 e3t_a(:,:,:) = e3t_b(:,:,:) 256 ! Reconstruction of all vertical scale factors at now and before time-steps 257 ! ============================================================================= 258 ! Horizontal scale factor interpolations 259 ! -------------------------------------- 260 CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) 261 CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) 262 CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 263 CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 264 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 265 ! Vertical scale factor interpolations 266 ! ------------------------------------ 267 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) 268 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 269 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 270 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 271 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 272 ! t- and w- points depth 273 ! ---------------------- 274 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 275 gdepw_n(:,:,1) = 0.0_wp 276 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 277 DO jk = 2, jpk 278 gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk) 279 gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1) 280 gde3w_n(:,:,jk) = gdept_n(:,:,jk ) - sshn (:,:) 281 END DO 282 ENDIF 285 283 ENDIF 286 284 ENDIF 287 285 ! 288 286 CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 ) 289 287 ! … … 448 446 ! Freezing/melting potential 449 447 ! Calculated over NEMO leapfrog timestep (hence 2*dt) 450 nfrzmlt(:,:) =rau0*rcp*fse3t_m(:,:)*(Tocnfrz-sst_m(:,:))/(2.0*dt)448 nfrzmlt(:,:) = rau0 * rcp * e3t_m(:,:) * ( Tocnfrz-sst_m(:,:) ) / ( 2.0*dt ) 451 449 452 450 ztmp(:,:) = nfrzmlt(:,:) -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_if.F90
r5541 r6140 34 34 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_ice ! structure of input ice-cover (file informations, fields read) 35 35 36 !! * Substitutions37 # include "domzgr_substitute.h90"38 36 !!---------------------------------------------------------------------- 39 37 !! NEMO/OPA 3.3 , NEMO Consortium (2010) -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
r5541 r6140 25 25 USE thd_ice ! LIM-3: thermodynamical variables 26 26 USE dom_ice ! LIM-3: ice domain 27 27 ! 28 28 USE sbc_oce ! Surface boundary condition: ocean fields 29 29 USE sbc_ice ! Surface boundary condition: ice fields … … 32 32 USE sbccpl ! Surface boundary condition: coupled interface 33 33 USE albedo ! ocean & ice albedo 34 34 ! 35 35 USE phycst ! Define parameters for the routines 36 36 USE eosbn2 ! equation of state … … 47 47 USE limupdate2 ! update of global variables 48 48 USE limvar ! Ice variables switch 49 49 USE limctl ! 50 50 USE limmsh ! LIM mesh 51 51 USE limistate ! LIM initial state 52 52 USE limthd_sal ! LIM ice thermodynamics: salinity 53 53 ! 54 54 USE c1d ! 1D vertical configuration 55 USE in_out_manager ! I/O manager 56 USE iom ! I/O manager library 57 USE prtctl ! Print control 58 USE lib_fortran ! 55 59 USE lbclnk ! lateral boundary condition - MPP link 56 60 USE lib_mpp ! MPP library 57 61 USE wrk_nemo ! work arrays 58 62 USE timing ! Timing 59 USE iom ! I/O manager library60 USE in_out_manager ! I/O manager61 USE prtctl ! Print control62 USE lib_fortran !63 USE limctl64 63 65 64 #if defined key_bdy … … 74 73 75 74 !! * Substitutions 76 # include "domzgr_substitute.h90"77 75 # include "vectopt_loop_substitute.h90" 78 76 !!---------------------------------------------------------------------- … … 82 80 !!---------------------------------------------------------------------- 83 81 CONTAINS 84 85 !!======================================================================86 82 87 83 SUBROUTINE sbc_ice_lim( kt, kblk ) … … 270 266 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ via Louvain la Neuve Ice Model (LIM-3) time stepping' 271 267 ! 272 268 ! ! Open the reference and configuration namelist files and namelist output file 273 269 CALL ctl_opn( numnam_ice_ref, 'namelist_ice_ref', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 274 270 CALL ctl_opn( numnam_ice_cfg, 'namelist_ice_cfg', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 275 271 IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 ) 276 272 ! 277 273 CALL ice_run ! set some ice run parameters 278 274 ! … … 348 344 REWIND( numnam_ice_ref ) ! Namelist namicerun in reference namelist : Parameters for ice 349 345 READ ( numnam_ice_ref, namicerun, IOSTAT = ios, ERR = 901) 350 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp )351 346 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp ) 347 ! 352 348 REWIND( numnam_ice_cfg ) ! Namelist namicerun in configuration namelist : Parameters for ice 353 349 READ ( numnam_ice_cfg, namicerun, IOSTAT = ios, ERR = 902 ) 354 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp ) 355 IF(lwm) WRITE ( numoni, namicerun ) 356 ! 350 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp ) 351 IF(lwm) WRITE( numoni, namicerun ) 357 352 ! 358 353 IF(lwp) THEN ! control print … … 373 368 ! 374 369 ! sea-ice timestep and inverse 375 rdt_ice = nn_fsbc * rdt tra(1)370 rdt_ice = nn_fsbc * rdt 376 371 r1_rdtice = 1._wp / rdt_ice 377 372 … … 405 400 REWIND( numnam_ice_ref ) ! Namelist namiceitd in reference namelist : Parameters for ice 406 401 READ ( numnam_ice_ref, namiceitd, IOSTAT = ios, ERR = 903) 407 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in reference namelist', lwp )408 402 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in reference namelist', lwp ) 403 ! 409 404 REWIND( numnam_ice_cfg ) ! Namelist namiceitd in configuration namelist : Parameters for ice 410 405 READ ( numnam_ice_cfg, namiceitd, IOSTAT = ios, ERR = 904 ) 411 904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in configuration namelist', lwp ) 412 IF(lwm) WRITE ( numoni, namiceitd ) 413 ! 406 904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in configuration namelist', lwp ) 407 IF(lwm) WRITE( numoni, namiceitd ) 414 408 ! 415 409 IF(lwp) THEN ! control print … … 417 411 WRITE(numout,*) 'ice_itd : ice cat distribution' 418 412 WRITE(numout,*) ' ~~~~~~' 419 WRITE(numout,*) ' shape of ice categories distribution 420 WRITE(numout,*) ' mean ice thickness in the domain ( only active if nn_catbnd=2)rn_himean = ', rn_himean413 WRITE(numout,*) ' shape of ice categories distribution nn_catbnd = ', nn_catbnd 414 WRITE(numout,*) ' mean ice thickness in the domain (used if nn_catbnd=2) rn_himean = ', rn_himean 421 415 ENDIF 422 416 ! 423 417 !---------------------------------- 424 418 !- Thickness categories boundaries … … 427 421 IF(lwp) WRITE(numout,*) 'lim_itd_init : Initialization of ice cat distribution ' 428 422 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 429 423 ! 430 424 hi_max(:) = 0._wp 431 432 SELECT CASE ( nn_catbnd ) 433 !---------------------- 434 CASE (1) ! tanh function (CICE) 435 !---------------------- 425 ! 426 SELECT CASE ( nn_catbnd ) ! type of ice categories distribution 427 ! 428 CASE (1) !== tanh function (CICE) ==! 436 429 zc1 = 3._wp / REAL( jpl, wp ) 437 430 zc2 = 10._wp * zc1 438 431 zc3 = 3._wp 439 440 432 DO jl = 1, jpl 441 433 zx1 = REAL( jl-1, wp ) / REAL( jpl, wp ) 442 434 hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1._wp + TANH( zc3 * (zx1 - 1._wp ) ) ) 443 435 END DO 444 445 !---------------------- 446 CASE (2) ! h^(-alpha) function 447 !---------------------- 448 zalpha = 0.05 ! exponent of the transform function 449 450 zhmax = 3.*rn_himean 451 436 ! 437 CASE (2) !== h^(-alpha) function ==! 438 zalpha = 0.05_wp 439 zhmax = 3._wp * rn_himean 452 440 DO jl = 1, jpl 453 441 znum = jpl * ( zhmax+1 )**zalpha 454 zden = ( jpl - jl ) * ( zhmax+1 )**zalpha + jl442 zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp ) 455 443 hi_max(jl) = ( znum / zden )**(1./zalpha) - 1 456 444 END DO 457 445 ! 458 446 END SELECT 459 460 DO jl = 1, jpl 447 ! 448 DO jl = 1, jpl ! mean thickness by category 461 449 hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 462 450 END DO 463 464 ! Set hi_max(jpl) to a big value to ensure that all ice is thinner than hi_max(jpl) 465 hi_max(jpl) = 99._wp 466 451 ! 452 hi_max(jpl) = 99._wp ! set to a big value to ensure that all ice is thinner than hi_max(jpl) 453 ! 467 454 IF(lwp) WRITE(numout,*) ' Thickness category boundaries ' 468 455 IF(lwp) WRITE(numout,*) ' hi_max ', hi_max(0:jpl) … … 471 458 472 459 473 SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_limflx ) 460 SUBROUTINE ice_lim_flx( ptn_ice , palb_ice, pqns_ice , & 461 & pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_limflx ) 474 462 !!--------------------------------------------------------------------- 475 463 !! *** ROUTINE ice_lim_flx *** … … 483 471 !!--------------------------------------------------------------------- 484 472 INTEGER , INTENT(in ) :: k_limflx ! =-1 do nothing; =0 average ; 485 473 ! ! =1 average and redistribute ; =2 redistribute 486 474 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptn_ice ! ice surface temperature 487 475 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: palb_ice ! ice albedo … … 503 491 REAL(wp), POINTER, DIMENSION(:,:) :: z_devap_m ! Mean d(evap)/dT over all categories 504 492 !!---------------------------------------------------------------------- 505 493 ! 506 494 IF( nn_timing == 1 ) CALL timing_start('ice_lim_flx') 507 !508 495 ! 509 496 SELECT CASE( k_limflx ) !== averaged on all ice categories ==! … … 529 516 CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m) 530 517 END SELECT 531 518 ! 532 519 SELECT CASE( k_limflx ) !== redistribution on all ice categories ==! 533 520 CASE( 1 , 2 ) … … 548 535 ! 549 536 END SUBROUTINE ice_lim_flx 537 550 538 551 539 SUBROUTINE sbc_lim_bef … … 564 552 u_ice_b(:,:) = u_ice(:,:) 565 553 v_ice_b(:,:) = v_ice(:,:) 566 554 ! 567 555 END SUBROUTINE sbc_lim_bef 556 568 557 569 558 SUBROUTINE sbc_lim_diag0 … … 580 569 sfx_bom(:,:) = 0._wp ; sfx_sum(:,:) = 0._wp 581 570 sfx_res(:,:) = 0._wp 582 571 ! 583 572 wfx_snw(:,:) = 0._wp ; wfx_ice(:,:) = 0._wp 584 573 wfx_sni(:,:) = 0._wp ; wfx_opw(:,:) = 0._wp … … 587 576 wfx_res(:,:) = 0._wp ; wfx_sub(:,:) = 0._wp 588 577 wfx_spr(:,:) = 0._wp ; 589 578 ! 590 579 hfx_thd(:,:) = 0._wp ; 591 580 hfx_snw(:,:) = 0._wp ; hfx_opw(:,:) = 0._wp … … 596 585 hfx_err(:,:) = 0._wp ; hfx_err_rem(:,:) = 0._wp 597 586 hfx_err_dif(:,:) = 0._wp ; 598 587 ! 599 588 afx_tot(:,:) = 0._wp ; 600 589 afx_dyn(:,:) = 0._wp ; afx_thd(:,:) = 0._wp 601 590 ! 602 591 diag_heat(:,:) = 0._wp ; diag_smvi(:,:) = 0._wp ; 603 592 diag_vice(:,:) = 0._wp ; diag_vsnw(:,:) = 0._wp ; 604 593 ! 605 594 END SUBROUTINE sbc_lim_diag0 606 595 … … 634 623 END FUNCTION fice_ice_ave 635 624 636 637 625 #else 638 626 !!---------------------------------------------------------------------- -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90
r5541 r6140 64 64 65 65 !! * Substitutions 66 # include "domzgr_substitute.h90"67 66 # include "vectopt_loop_substitute.h90" 68 67 !!---------------------------------------------------------------------- -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
r5836 r6140 35 35 ! public in order to be able to output then 36 36 37 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_tsc_b, risf_tsc 38 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qisf !: net heat flux from ice shelf37 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_tsc_b, risf_tsc !: before and now T & S isf contents [K.m/s & PSU.m/s] 38 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qisf !: net heat flux from ice shelf [W/m2] 39 39 REAL(wp), PUBLIC :: rn_hisf_tbl !: thickness of top boundary layer [m] 40 LOGICAL , PUBLIC :: ln_divisf !: flag to correct divergence 41 INTEGER , PUBLIC :: nn_isfblk !: 42 INTEGER , PUBLIC :: nn_gammablk !: 43 LOGICAL , PUBLIC :: ln_conserve !: 44 REAL(wp), PUBLIC :: rn_gammat0 !: temperature exchange coeficient 45 REAL(wp), PUBLIC :: rn_gammas0 !: salinity exchange coeficient 46 REAL(wp), PUBLIC :: rdivisf !: flag to test if fwf apply on divergence 40 INTEGER , PUBLIC :: nn_isf !: flag to choose between explicit/param/specified 41 INTEGER , PUBLIC :: nn_isfblk !: flag to choose the bulk formulation to compute the ice shelf melting 42 INTEGER , PUBLIC :: nn_gammablk !: flag to choose how the exchange coefficient is computed 43 REAL(wp), PUBLIC :: rn_gammat0 !: temperature exchange coeficient [] 44 REAL(wp), PUBLIC :: rn_gammas0 !: salinity exchange coeficient [] 47 45 48 46 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: rzisf_tbl !:depth of calving front (shallowest point) nn_isf ==2/3 49 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: rhisf_tbl, rhisf_tbl_0 !:thickness of tbl 47 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: rhisf_tbl, rhisf_tbl_0 !:thickness of tbl [m] 50 48 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: r1_hisf_tbl !:1/thickness of tbl 51 49 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: ralpha !:proportion of bottom cell influenced by tbl 52 50 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: risfLeff !:effective length (Leff) BG03 nn_isf==2 53 51 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: ttbl, stbl, utbl, vtbl !:top boundary layer variable at T point 54 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: misfkt, misfkb !:Level of ice shelf base55 56 REAL(wp), PUBLIC, SAVE :: rcpi = 2000.0_wp ! phycst ?57 REAL(wp), PUBLIC, SAVE :: kappa = 1.54e-6_wp ! phycst ?58 REAL(wp), PUBLIC, SAVE :: rhoisf = 920.0_wp ! phycst ?59 REAL(wp), PUBLIC, SAVE :: tsurf = -20.0_wp ! phycst ?60 REAL(wp), PUBLIC, SAVE :: lfusisf= 0.334e6_wp ! phycst ?52 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: misfkt, misfkb !:Level of ice shelf base 53 54 REAL(wp), PUBLIC, SAVE :: rcpi = 2000.0_wp ! specific heat of ice shelf [J/kg/K] 55 REAL(wp), PUBLIC, SAVE :: rkappa = 1.54e-6_wp ! heat diffusivity through the ice-shelf [m2/s] 56 REAL(wp), PUBLIC, SAVE :: rhoisf = 920.0_wp ! volumic mass of ice shelf [kg/m3] 57 REAL(wp), PUBLIC, SAVE :: tsurf = -20.0_wp ! air temperature on top of ice shelf [C] 58 REAL(wp), PUBLIC, SAVE :: rlfusisf = 0.334e6_wp ! latent heat of fusion of ice shelf [J/kg] 61 59 62 60 !: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3) 63 CHARACTER(len=100), PUBLIC :: cn_dirisf = './' !: Root directory for location of ssr files 64 TYPE(FLD_N) , PUBLIC :: sn_qisf, sn_fwfisf !: information about the runoff file to be read 65 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_qisf, sf_fwfisf 66 TYPE(FLD_N) , PUBLIC :: sn_rnfisf !: information about the runoff file to be read 67 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf 68 TYPE(FLD_N) , PUBLIC :: sn_depmax_isf, sn_depmin_isf, sn_Leff_isf !: information about the runoff file to be read 61 CHARACTER(len=100), PUBLIC :: cn_dirisf = './' !: Root directory for location of ssr files 62 TYPE(FLD_N) , PUBLIC :: sn_fwfisf !: information about the isf melting file to be read 63 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_fwfisf 64 TYPE(FLD_N) , PUBLIC :: sn_rnfisf !: information about the isf melting param. file to be read 65 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf 66 TYPE(FLD_N) , PUBLIC :: sn_depmax_isf !: information about the grounding line depth file to be read 67 TYPE(FLD_N) , PUBLIC :: sn_depmin_isf !: information about the calving line depth file to be read 68 TYPE(FLD_N) , PUBLIC :: sn_Leff_isf !: information about the effective length file to be read 69 69 70 !! * Substitutions71 # include "domzgr_substitute.h90"72 70 !!---------------------------------------------------------------------- 73 71 !! NEMO/OPA 3.7 , LOCEAN-IPSL (2015) … … 77 75 CONTAINS 78 76 79 77 SUBROUTINE sbc_isf(kt) 80 78 !!--------------------------------------------------------------------- 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 !! 95 NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, ln_divisf, ln_conserve, rn_gammat0, rn_gammas0, nn_gammablk, & 96 & sn_fwfisf, sn_qisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf 79 !! *** ROUTINE sbc_isf *** 80 !! 81 !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf 82 !! melting and freezing 83 !! 84 !! ** Method : 4 parameterizations are available according to nn_isf 85 !! nn_isf = 1 : Realistic ice_shelf formulation 86 !! 2 : Beckmann & Goose parameterization 87 !! 3 : Specified runoff in deptht (Mathiot & al. ) 88 !! 4 : specified fwf and heat flux forcing beneath the ice shelf 89 !!---------------------------------------------------------------------- 90 INTEGER, INTENT( in ) :: kt ! ocean time step 91 ! 92 INTEGER :: ji, jj ! loop index 93 REAL(wp), DIMENSION (:,:), POINTER :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep) 97 94 !!--------------------------------------------------------------------- 98 95 ! … … 100 97 IF( kt == nit000 ) THEN ! First call kt=nit000 ! 101 98 ! ! ====================== ! 102 REWIND( numnam_ref ) ! Namelist namsbc_rnf in reference namelist : Runoffs 103 READ ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901) 104 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp ) 105 106 REWIND( numnam_cfg ) ! Namelist namsbc_rnf in configuration namelist : Runoffs 107 READ ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 ) 108 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp ) 109 IF(lwm) WRITE ( numond, namsbc_isf ) 110 111 IF ( lwp ) WRITE(numout,*) 112 IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf' 113 IF ( lwp ) WRITE(numout,*) '~~~~~~~~~' 114 IF ( lwp ) WRITE(numout,*) 'sbcisf :' 115 IF ( lwp ) WRITE(numout,*) '~~~~~~~~' 116 IF ( lwp ) WRITE(numout,*) ' nn_isf = ', nn_isf 117 IF ( lwp ) WRITE(numout,*) ' nn_isfblk = ', nn_isfblk 118 IF ( lwp ) WRITE(numout,*) ' rn_hisf_tbl = ', rn_hisf_tbl 119 IF ( lwp ) WRITE(numout,*) ' ln_divisf = ', ln_divisf 120 IF ( lwp ) WRITE(numout,*) ' nn_gammablk = ', nn_gammablk 121 IF ( lwp ) WRITE(numout,*) ' rn_tfri2 = ', rn_tfri2 122 IF (ln_divisf) THEN ! keep it in the namelist ??? used true anyway as for runoff ? (PM) 123 rdivisf = 1._wp 124 ELSE 125 rdivisf = 0._wp 126 END IF 127 ! 128 ! Allocate public variable 129 IF ( sbc_isf_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' ) 130 ! 131 ! initialisation 132 qisf(:,:) = 0._wp ; fwfisf(:,:) = 0._wp 133 risf_tsc(:,:,:) = 0._wp 134 ! 135 ! define isf tbl tickness, top and bottom indice 136 IF (nn_isf == 1) THEN 137 rhisf_tbl(:,:) = rn_hisf_tbl 138 misfkt(:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv 139 ELSE IF ((nn_isf == 3) .OR. (nn_isf == 2)) THEN 140 ALLOCATE( sf_rnfisf(1), STAT=ierror ) 141 ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) ) 142 CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 143 144 !: read effective lenght (BG03) 145 IF (nn_isf == 2) THEN 146 ! Read Data and save some integral values 147 CALL iom_open( sn_Leff_isf%clname, inum ) 148 cvarLeff = 'soLeff' !: variable name for Efficient Length scale 149 CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1) 150 CALL iom_close(inum) 151 ! 152 risfLeff = risfLeff*1000 !: convertion in m 153 END IF 154 155 ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth) 156 CALL iom_open( sn_depmax_isf%clname, inum ) 157 cvarhisf = TRIM(sn_depmax_isf%clvar) 158 CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base 159 CALL iom_close(inum) 160 ! 161 CALL iom_open( sn_depmin_isf%clname, inum ) 162 cvarzisf = TRIM(sn_depmin_isf%clvar) 163 CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base 164 CALL iom_close(inum) 165 ! 166 rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:) !: tickness isf boundary layer 167 168 !! compute first level of the top boundary layer 169 DO ji = 1, jpi 170 DO jj = 1, jpj 171 jk = 2 172 DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdepw(ji,jj,jk) < rzisf_tbl(ji,jj) ) ; jk = jk + 1 ; END DO 173 misfkt(ji,jj) = jk-1 174 END DO 175 END DO 176 177 ELSE IF ( nn_isf == 4 ) THEN 178 ! as in nn_isf == 1 179 rhisf_tbl(:,:) = rn_hisf_tbl 180 misfkt(:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv 181 182 ! load variable used in fldread (use for temporal interpolation of isf fwf forcing) 183 ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror ) 184 ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) ) 185 ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) ) 186 CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 187 !CALL fld_fill( sf_qisf , (/ sn_qisf /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data' , 'namsbc_isf' ) 188 END IF 189 190 ! compute bottom level of isf tbl and thickness of tbl below the ice shelf 191 rhisf_tbl_0(:,:) = rhisf_tbl(:,:) 192 DO jj = 1,jpj 193 DO ji = 1,jpi 194 ikt = misfkt(ji,jj) 195 ikb = misfkt(ji,jj) 196 ! thickness of boundary layer at least the top level thickness 197 rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt)) 198 199 ! determine the deepest level influenced by the boundary layer 200 ! test on tmask useless ????? 201 DO jk = ikt, mbkt(ji,jj) 202 IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 203 END DO 204 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 205 misfkb(ji,jj) = ikb ! last wet level of the tbl 206 r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 207 208 zhk = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 209 ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer 210 END DO 211 END DO 99 CALL sbc_isf_init 100 ! ! ---------------------------------------- ! 101 ELSE ! Swap of forcing fields ! 102 ! ! ---------------------------------------- ! 103 fwfisf_b (:,: ) = fwfisf (:,: ) ! Swap the ocean forcing fields except at nit000 104 risf_tsc_b(:,:,:) = risf_tsc(:,:,:) ! where before fields are set at the end of the routine 212 105 ! 213 106 END IF 214 107 215 ! ! ---------------------------------------- !216 IF( kt /= nit000 ) THEN ! Swap of forcing fields !217 ! ! ---------------------------------------- !218 fwfisf_b (:,: ) = fwfisf (:,: ) ! Swap the ocean forcing fields except at nit000219 risf_tsc_b(:,:,:) = risf_tsc(:,:,:) ! where before fields are set at the end of the routine220 !221 ENDIF222 223 108 IF( MOD( kt-1, nn_fsbc) == 0 ) THEN 224 225 226 ! compute salf and heat flux 227 IF (nn_isf == 1) THEN 228 ! realistic ice shelf formulation 109 ! allocation 110 CALL wrk_alloc( jpi,jpj, zt_frz, zdep ) 111 112 ! compute salt and heat flux 113 SELECT CASE ( nn_isf ) 114 CASE ( 1 ) ! realistic ice shelf formulation 229 115 ! compute T/S/U/V for the top boundary layer 230 116 CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T') 231 117 CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T') 232 CALL sbc_isf_tbl(un(:,:,:) ,utbl(:,:),'U')233 CALL sbc_isf_tbl(vn(:,:,:) ,vtbl(:,:),'V')118 CALL sbc_isf_tbl(un(:,:,:) ,utbl(:,:),'U') 119 CALL sbc_isf_tbl(vn(:,:,:) ,vtbl(:,:),'V') 234 120 ! iom print 235 121 CALL iom_put('ttbl',ttbl(:,:)) 236 122 CALL iom_put('stbl',stbl(:,:)) 237 CALL iom_put('utbl',utbl(:,:) )238 CALL iom_put('vtbl',vtbl(:,:) )123 CALL iom_put('utbl',utbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)) 124 CALL iom_put('vtbl',vtbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)) 239 125 ! compute fwf and heat flux 240 126 CALL sbc_isf_cav (kt) 241 127 242 ELSE IF (nn_isf == 2) THEN 243 ! Beckmann and Goosse parametrisation 128 CASE ( 2 ) ! Beckmann and Goosse parametrisation 244 129 stbl(:,:) = soce 245 130 CALL sbc_isf_bg03(kt) 246 131 247 ELSE IF (nn_isf == 3) THEN 248 ! specified runoff in depth (Mathiot et al., XXXX in preparation) 132 CASE ( 3 ) ! specified runoff in depth (Mathiot et al., XXXX in preparation) 249 133 CALL fld_read ( kt, nn_fsbc, sf_rnfisf ) 250 fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1) ! f resh waterflux from the isf (fwfisf <0 mean melting)251 qisf(:,:) = fwfisf(:,:) * lfusisf ! heatflux134 fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1) ! fwf flux from the isf (fwfisf <0 mean melting) 135 qisf(:,:) = fwfisf(:,:) * rlfusisf ! heat flux 252 136 stbl(:,:) = soce 253 137 254 ELSE IF (nn_isf == 4) THEN 255 ! specified fwf and heat flux forcing beneath the ice shelf 138 CASE ( 4 ) ! specified fwf and heat flux forcing beneath the ice shelf 256 139 CALL fld_read ( kt, nn_fsbc, sf_fwfisf ) 257 !CALL fld_read ( kt, nn_fsbc, sf_qisf ) 258 fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1) ! fwf 259 qisf(:,:) = fwfisf(:,:) * lfusisf ! heat flux 260 !qisf(:,:) = sf_qisf(1)%fnow(:,:,1) ! heat flux 140 fwfisf(:,:) = - sf_fwfisf(1)%fnow(:,:,1) ! fwf flux from the isf (fwfisf <0 mean melting) 141 qisf(:,:) = fwfisf(:,:) * rlfusisf ! heat flux 261 142 stbl(:,:) = soce 262 143 263 END IF 144 END SELECT 145 264 146 ! compute tsc due to isf 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 ! 147 ! isf melting implemented as a volume flux and we assume that melt water is at 0 PSU. 148 ! WARNING water add at temp = 0C, need to add a correction term (fwfisf * tfreez / rau0). 149 ! compute freezing point beneath ice shelf (or top cell if nn_isf = 3) 150 DO jj = 1,jpj 151 DO ji = 1,jpi 152 zdep(ji,jj)=gdepw_n(ji,jj,misfkt(ji,jj)) 153 END DO 154 END DO 155 CALL eos_fzp( stbl(:,:), zt_frz(:,:), zdep(:,:) ) 269 156 270 ! salt effect already take into account in vertical advection 271 risf_tsc(:,:,jp_sal) = (1.0_wp-rdivisf) * fwfisf(:,:) * stbl(:,:) * r1_rau0 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 157 risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rau0 ! 158 risf_tsc(:,:,jp_sal) = 0.0_wp 159 280 160 ! lbclnk 281 161 CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.) 282 162 CALL lbc_lnk(risf_tsc(:,:,jp_sal),'T',1.) 283 CALL lbc_lnk(fwfisf(:,:) ,'T',1.)284 CALL lbc_lnk(qisf(:,:) ,'T',1.)285 286 IF( kt == nit000 ) THEN 163 CALL lbc_lnk(fwfisf(:,:) ,'T',1.) 164 CALL lbc_lnk(qisf(:,:) ,'T',1.) 165 166 IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 ! 287 167 IF( ln_rstart .AND. & ! Restart: read in restart file 288 168 & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN … … 295 175 risf_tsc_b(:,:,:)= risf_tsc(:,:,:) 296 176 END IF 297 END IF177 END IF 298 178 ! 179 ! output 180 CALL iom_put('qisf' , qisf) 181 CALL iom_put('fwfisf', fwfisf) 182 183 ! deallocation 184 CALL wrk_dealloc( jpi,jpj, zt_frz, zdep ) 299 185 END IF 300 186 ! … … 315 201 & STAT= sbc_isf_alloc ) 316 202 ! 317 IF( lk_mpp 203 IF( lk_mpp ) CALL mpp_sum ( sbc_isf_alloc ) 318 204 IF( sbc_isf_alloc /= 0 ) CALL ctl_warn('sbc_isf_alloc: failed to allocate arrays.') 319 205 ! 320 END IF206 END IF 321 207 END FUNCTION 322 208 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 !! 209 SUBROUTINE sbc_isf_init 210 !!--------------------------------------------------------------------- 211 !! *** ROUTINE sbc_isf_init *** 212 !! 213 !! ** Purpose : Initialisation of variables for iceshelf fluxes formulation 214 !! 215 !! ** Method : 4 parameterizations are available according to nn_isf 216 !! nn_isf = 1 : Realistic ice_shelf formulation 217 !! 2 : Beckmann & Goose parameterization 218 !! 3 : Specified runoff in deptht (Mathiot & al. ) 219 !! 4 : specified fwf and heat flux forcing beneath the ice shelf 220 !!---------------------------------------------------------------------- 221 INTEGER :: ji, jj, jk ! loop index 222 INTEGER :: ik ! current level index 223 INTEGER :: ikt, ikb ! top and bottom level of the isf boundary layer 224 INTEGER :: inum, ierror 225 INTEGER :: ios ! Local integer output status for namelist read 226 REAL(wp) :: zhk 227 CHARACTER(len=256) :: cvarzisf, cvarhisf ! name for isf file 228 CHARACTER(LEN=32 ) :: cvarLeff ! variable name for efficient Length scale 229 !!---------------------------------------------------------------------- 230 NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf, & 231 & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf 232 !!---------------------------------------------------------------------- 233 234 REWIND( numnam_ref ) ! Namelist namsbc_rnf in reference namelist : Runoffs 235 READ ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901) 236 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp ) 237 238 REWIND( numnam_cfg ) ! Namelist namsbc_rnf in configuration namelist : Runoffs 239 READ ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 ) 240 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp ) 241 IF(lwm) WRITE ( numond, namsbc_isf ) 242 243 IF ( lwp ) WRITE(numout,*) 244 IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf' 245 IF ( lwp ) WRITE(numout,*) '~~~~~~~~~' 246 IF ( lwp ) WRITE(numout,*) 'sbcisf :' 247 IF ( lwp ) WRITE(numout,*) '~~~~~~~~' 248 IF ( lwp ) WRITE(numout,*) ' nn_isf = ', nn_isf 249 IF ( lwp ) WRITE(numout,*) ' nn_isfblk = ', nn_isfblk 250 IF ( lwp ) WRITE(numout,*) ' rn_hisf_tbl = ', rn_hisf_tbl 251 IF ( lwp ) WRITE(numout,*) ' nn_gammablk = ', nn_gammablk 252 IF ( lwp ) WRITE(numout,*) ' rn_gammat0 = ', rn_gammat0 253 IF ( lwp ) WRITE(numout,*) ' rn_gammas0 = ', rn_gammas0 254 IF ( lwp ) WRITE(numout,*) ' rn_tfri2 = ', rn_tfri2 255 ! 256 ! Allocate public variable 257 IF ( sbc_isf_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' ) 258 ! 259 ! initialisation 260 qisf(:,:) = 0._wp ; fwfisf (:,:) = 0._wp 261 risf_tsc(:,:,:) = 0._wp ; fwfisf_b(:,:) = 0._wp 262 ! 263 ! define isf tbl tickness, top and bottom indice 264 SELECT CASE ( nn_isf ) 265 CASE ( 1 ) 266 rhisf_tbl(:,:) = rn_hisf_tbl 267 misfkt(:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv 268 269 CASE ( 2 , 3 ) 270 ALLOCATE( sf_rnfisf(1), STAT=ierror ) 271 ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) ) 272 CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 273 274 ! read effective lenght (BG03) 275 IF (nn_isf == 2) THEN 276 CALL iom_open( sn_Leff_isf%clname, inum ) 277 cvarLeff = TRIM(sn_Leff_isf%clvar) 278 CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1) 279 CALL iom_close(inum) 280 ! 281 risfLeff = risfLeff*1000.0_wp !: convertion in m 282 END IF 283 284 ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth) 285 CALL iom_open( sn_depmax_isf%clname, inum ) 286 cvarhisf = TRIM(sn_depmax_isf%clvar) 287 CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base 288 CALL iom_close(inum) 289 ! 290 CALL iom_open( sn_depmin_isf%clname, inum ) 291 cvarzisf = TRIM(sn_depmin_isf%clvar) 292 CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base 293 CALL iom_close(inum) 294 ! 295 rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:) !: tickness isf boundary layer 296 297 !! compute first level of the top boundary layer 298 DO ji = 1, jpi 299 DO jj = 1, jpj 300 ik = 2 301 DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw_n(ji,jj,ik) < rzisf_tbl(ji,jj) ) ; ik = ik + 1 ; END DO 302 misfkt(ji,jj) = ik-1 303 END DO 304 END DO 305 306 CASE ( 4 ) 307 ! as in nn_isf == 1 308 rhisf_tbl(:,:) = rn_hisf_tbl 309 misfkt(:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv 310 311 ! load variable used in fldread (use for temporal interpolation of isf fwf forcing) 312 ALLOCATE( sf_fwfisf(1), STAT=ierror ) 313 ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) ) 314 CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 315 316 END SELECT 317 318 rhisf_tbl_0(:,:) = rhisf_tbl(:,:) 319 320 ! compute bottom level of isf tbl and thickness of tbl below the ice shelf 321 DO jj = 1,jpj 322 DO ji = 1,jpi 323 ikt = misfkt(ji,jj) 324 ikb = misfkt(ji,jj) 325 ! thickness of boundary layer at least the top level thickness 326 rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt)) 327 328 ! determine the deepest level influenced by the boundary layer 329 DO jk = ikt+1, mbkt(ji,jj) 330 IF ( (SUM(e3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 331 END DO 332 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 333 misfkb(ji,jj) = ikb ! last wet level of the tbl 334 r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 335 336 zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 337 ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer 338 END DO 339 END DO 340 341 END SUBROUTINE sbc_isf_init 342 343 SUBROUTINE sbc_isf_bg03(kt) 344 !!--------------------------------------------------------------------- 345 !! *** ROUTINE sbc_isf_bg03 *** 346 !! 347 !! ** Purpose : add net heat and fresh water flux from ice shelf melting 348 !! into the adjacent ocean 349 !! 350 !! ** Method : See reference 351 !! 352 !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean 353 !! interaction for climate models", Ocean Modelling 5(2003) 157-170. 354 !! (hereafter BG) 340 355 !! History : 341 !! !06-02 (C. Wang) Original code356 !! 06-02 (C. Wang) Original code 342 357 !!---------------------------------------------------------------------- 343 358 INTEGER, INTENT ( in ) :: kt 344 359 ! 345 INTEGER :: ji, jj, jk, jish !temporary integer 346 INTEGER :: ijkmin 347 INTEGER :: ii, ij, ik 348 INTEGER :: inum 349 350 REAL(wp) :: zt_sum ! sum of the temperature between 200m and 600m 351 REAL(wp) :: zt_ave ! averaged temperature between 200m and 600m 352 REAL(wp) :: zt_frz ! freezing point temperature at depth z 353 REAL(wp) :: zpress ! pressure to compute the freezing point in depth 354 355 !!---------------------------------------------------------------------- 356 IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03') 357 ! 358 359 ! This test is false only in the very first time step of a run (JMM ???- Initialy build to skip 1rst year of run ) 360 DO ji = 1, jpi 361 DO jj = 1, jpj 362 ik = misfkt(ji,jj) 363 !! Initialize arrays to 0 (each step) 364 zt_sum = 0.e0_wp 365 IF ( ik .GT. 1 ) THEN 366 ! 3. -----------the average temperature between 200m and 600m --------------------- 367 DO jk = misfkt(ji,jj),misfkb(ji,jj) 368 ! freezing point temperature at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!) 369 ! after verif with UNESCO, wrong sign in BG eq. 2 370 ! Calculate freezing temperature 371 zpress = grav*rau0*fsdept(ji,jj,ik)*1.e-04 372 CALL eos_fzp(tsb(ji,jj,ik,jp_sal), zt_frz, zpress) 373 zt_sum = zt_sum + (tsn(ji,jj,ik,jp_tem)-zt_frz) * fse3t(ji,jj,ik) * tmask(ji,jj,ik) ! sum temp 374 ENDDO 375 zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value 376 377 ! 4. ------------Net heat flux and fresh water flux due to the ice shelf 378 ! For those corresponding to zonal boundary 379 qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave & 380 & / (e1t(ji,jj) * e2t(ji,jj)) * tmask(ji,jj,ik) 360 INTEGER :: ji, jj, jk ! dummy loop index 361 INTEGER :: ik ! current level 362 REAL(wp) :: zt_sum ! sum of the temperature between 200m and 600m 363 REAL(wp) :: zt_ave ! averaged temperature between 200m and 600m 364 REAL(wp) :: zt_frz ! freezing point temperature at depth z 365 REAL(wp) :: zpress ! pressure to compute the freezing point in depth 366 !!---------------------------------------------------------------------- 367 368 IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03') 369 ! 370 DO ji = 1, jpi 371 DO jj = 1, jpj 372 ik = misfkt(ji,jj) 373 !! Initialize arrays to 0 (each step) 374 zt_sum = 0.e0_wp 375 IF ( ik > 1 ) THEN 376 ! 1. -----------the average temperature between 200m and 600m --------------------- 377 DO jk = misfkt(ji,jj),misfkb(ji,jj) 378 ! freezing point temperature at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!) 379 ! after verif with UNESCO, wrong sign in BG eq. 2 380 ! Calculate freezing temperature 381 CALL eos_fzp(stbl(ji,jj), zt_frz, zpress) 382 zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * e3t_n(ji,jj,jk) * tmask(ji,jj,jk) ! sum temp 383 END DO 384 zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value 385 ! 2. ------------Net heat flux and fresh water flux due to the ice shelf 386 ! For those corresponding to zonal boundary 387 qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave & 388 & * r1_e1e2t(ji,jj) * tmask(ji,jj,jk) 381 389 382 fwfisf(ji,jj) = qisf(ji,jj) /lfusisf !fresh water flux kg/(m2s)383 fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) )384 !add to salinity trend385 ELSE386 qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp387 END IF388 END DO389 END DO390 !391 IF( nn_timing == 1 ) CALL timing_stop('sbc_isf_bg03')390 fwfisf(ji,jj) = qisf(ji,jj) / rlfusisf !fresh water flux kg/(m2s) 391 fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) ) 392 !add to salinity trend 393 ELSE 394 qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp 395 END IF 396 END DO 397 END DO 398 ! 399 IF( nn_timing == 1 ) CALL timing_stop('sbc_isf_bg03') 392 400 ! 393 401 END SUBROUTINE sbc_isf_bg03 394 402 395 396 SUBROUTINE sbc_isf_cav( kt ) 403 SUBROUTINE sbc_isf_cav( kt ) 397 404 !!--------------------------------------------------------------------- 398 405 !! *** ROUTINE sbc_isf_cav *** … … 409 416 INTEGER, INTENT(in) :: kt ! ocean time step 410 417 ! 411 LOGICAL :: ln_isomip = .true. 412 REAL(wp), DIMENSION(:,:), POINTER :: zfrz,zpress,zti 413 REAL(wp), DIMENSION(:,:), POINTER :: zgammat2d, zgammas2d 414 !REAL(wp), DIMENSION(:,:), POINTER :: zqisf, zfwfisf 418 INTEGER :: ji, jj ! dummy loop indices 419 INTEGER :: nit 415 420 REAL(wp) :: zlamb1, zlamb2, zlamb3 416 421 REAL(wp) :: zeps1,zeps2,zeps3,zeps4,zeps6,zeps7 417 422 REAL(wp) :: zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac 418 REAL(wp) :: zfwflx, zhtflx, zhtflx_b 419 REAL(wp) :: zgammat, zgammas 420 REAL(wp) :: zeps = -1.e-20_wp !== Local constant initialization ==! 421 INTEGER :: ji, jj ! dummy loop indices 422 INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers 423 INTEGER :: ierror ! return error code 424 LOGICAL :: lit=.TRUE. 425 INTEGER :: nit 423 REAL(wp) :: zeps = 1.e-20_wp 424 REAL(wp) :: zerr 425 REAL(wp), DIMENSION(:,:), POINTER :: zfrz 426 REAL(wp), DIMENSION(:,:), POINTER :: zgammat, zgammas 427 REAL(wp), DIMENSION(:,:), POINTER :: zfwflx, zhtflx, zhtflx_b 428 LOGICAL :: lit 426 429 !!--------------------------------------------------------------------- 427 ! 428 ! coeficient for linearisation of tfreez429 zlamb1 =-0.0575430 zlamb2 =0.0901431 zlamb3 =-7.61e-04430 ! coeficient for linearisation of potential tfreez 431 ! Crude approximation for pressure (but commonly used) 432 zlamb1 =-0.0573_wp 433 zlamb2 = 0.0832_wp 434 zlamb3 =-7.53e-08_wp * grav * rau0 432 435 IF( nn_timing == 1 ) CALL timing_start('sbc_isf_cav') 433 436 ! 434 CALL wrk_alloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d ) 435 436 zcfac=0.0_wp 437 IF (ln_conserve) zcfac=1.0_wp 438 zpress(:,:)=0.0_wp 439 zgammat2d(:,:)=0.0_wp 440 zgammas2d(:,:)=0.0_wp 441 ! 442 ! 443 DO jj = 1, jpj 444 DO ji = 1, jpi 445 ! Crude approximation for pressure (but commonly used) 446 ! 1e-04 to convert from Pa to dBar 447 zpress(ji,jj)=grav*rau0*fsdepw(ji,jj,mikt(ji,jj))*1.e-04 448 ! 449 END DO 437 CALL wrk_alloc( jpi,jpj, zfrz , zgammat, zgammas ) 438 CALL wrk_alloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b ) 439 440 ! initialisation 441 zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0 442 zhtflx (:,:) = 0.0_wp ; zhtflx_b(:,:) = 0.0_wp 443 zfwflx (:,:) = 0.0_wp 444 445 ! compute ice shelf melting 446 nit = 1 ; lit = .TRUE. 447 DO WHILE ( lit ) ! maybe just a constant number of iteration as in blk_core is fine 448 SELECT CASE ( nn_isfblk ) 449 CASE ( 1 ) ! ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006) 450 ! Calculate freezing temperature 451 CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) ) 452 453 ! compute gammat every where (2d) 454 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 455 456 ! compute upward heat flux zhtflx and upward water flux zwflx 457 DO jj = 1, jpj 458 DO ji = 1, jpi 459 zhtflx(ji,jj) = zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-zfrz(ji,jj)) 460 zfwflx(ji,jj) = - zhtflx(ji,jj)/rlfusisf 461 END DO 462 END DO 463 464 ! Compute heat flux and upward fresh water flux 465 qisf (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 466 fwfisf(:,:) = zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 467 468 CASE ( 2 ) ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015) 469 ! compute gammat every where (2d) 470 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 471 472 ! compute upward heat flux zhtflx and upward water flux zwflx 473 ! Resolution of a 2d equation from equation 21, 22 and 23 to find Sb (Asay-Davis et al., 2015) 474 DO jj = 1, jpj 475 DO ji = 1, jpi 476 ! compute coeficient to solve the 2nd order equation 477 zeps1 = rcp*rau0*zgammat(ji,jj) 478 zeps2 = rlfusisf*rau0*zgammas(ji,jj) 479 zeps3 = rhoisf*rcpi*rkappa/MAX(risfdep(ji,jj),zeps) 480 zeps4 = zlamb2+zlamb3*risfdep(ji,jj) 481 zeps6 = zeps4-ttbl(ji,jj) 482 zeps7 = zeps4-tsurf 483 zaqe = zlamb1 * (zeps1 + zeps3) 484 zaqer = 0.5_wp/MIN(zaqe,-zeps) 485 zbqe = zeps1*zeps6+zeps3*zeps7-zeps2 486 zcqe = zeps2*stbl(ji,jj) 487 zdis = zbqe*zbqe-4.0_wp*zaqe*zcqe 488 489 ! Presumably zdis can never be negative because gammas is very small compared to gammat 490 ! compute s freeze 491 zsfrz=(-zbqe-SQRT(zdis))*zaqer 492 IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer 493 494 ! compute t freeze (eq. 22) 495 zfrz(ji,jj)=zeps4+zlamb1*zsfrz 496 497 ! zfwflx is upward water flux 498 ! zhtflx is upward heat flux (out of ocean) 499 ! compute the upward water and heat flux (eq. 28 and eq. 29) 500 zfwflx(ji,jj) = rau0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps) 501 zhtflx(ji,jj) = zgammat(ji,jj) * rau0 * rcp * (ttbl(ji,jj) - zfrz(ji,jj) ) 502 END DO 503 END DO 504 505 ! compute heat and water flux 506 qisf (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 507 fwfisf(:,:) = zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 508 509 END SELECT 510 511 ! define if we need to iterate (nn_gammablk 0/1 do not need iteration) 512 IF ( nn_gammablk < 2 ) THEN ; lit = .FALSE. 513 ELSE 514 ! check total number of iteration 515 IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 516 ELSE ; nit = nit + 1 517 END IF 518 519 ! compute error between 2 iterations 520 ! if needed save gammat and compute zhtflx_b for next iteration 521 zerr = MAXVAL(ABS(zhtflx-zhtflx_b)) 522 IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE. 523 ELSE ; zhtflx_b(:,:) = zhtflx(:,:) 524 END IF 525 END IF 450 526 END DO 451 452 ! Calculate in-situ temperature (ref to surface) 453 zti(:,:)=tinsitu( ttbl, stbl, zpress ) 454 ! Calculate freezing temperature 455 CALL eos_fzp( sss_m(:,:), zfrz(:,:), zpress ) 456 457 458 zhtflx=0._wp ; zfwflx=0._wp 459 IF (nn_isfblk == 1) THEN 460 DO jj = 1, jpj 461 DO ji = 1, jpi 462 IF (mikt(ji,jj) > 1 ) THEN 463 nit = 1; lit = .TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp 464 DO WHILE ( lit ) 465 ! compute gamma 466 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit) 467 ! zhtflx is upward heat flux (out of ocean) 468 zhtflx = zgammat*rcp*rau0*(zti(ji,jj)-zfrz(ji,jj)) 469 ! zwflx is upward water flux 470 zfwflx = - zhtflx/lfusisf 471 ! test convergence and compute gammat 472 IF ( (zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE. 473 474 nit = nit + 1 475 IF (nit .GE. 100) CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 476 477 ! save gammat and compute zhtflx_b 478 zgammat2d(ji,jj)=zgammat 479 zhtflx_b = zhtflx 480 END DO 481 482 qisf(ji,jj) = - zhtflx 483 ! For genuine ISOMIP protocol this should probably be something like 484 fwfisf(ji,jj) = zfwflx * ( soce / MAX(stbl(ji,jj),zeps)) 485 ELSE 486 fwfisf(ji,jj) = 0._wp 487 qisf(ji,jj) = 0._wp 488 END IF 489 ! 490 END DO 491 END DO 492 493 ELSE IF (nn_isfblk == 2 ) THEN 494 495 ! More complicated 3 equation thermodynamics as in MITgcm 496 DO jj = 2, jpj 497 DO ji = 2, jpi 498 IF (mikt(ji,jj) > 1 ) THEN 499 nit=1; lit=.TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp; zhtflx=0._wp 500 DO WHILE ( lit ) 501 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit) 502 503 zeps1=rcp*rau0*zgammat 504 zeps2=lfusisf*rau0*zgammas 505 zeps3=rhoisf*rcpi*kappa/risfdep(ji,jj) 506 zeps4=zlamb2+zlamb3*risfdep(ji,jj) 507 zeps6=zeps4-zti(ji,jj) 508 zeps7=zeps4-tsurf 509 zaqe=zlamb1 * (zeps1 + zeps3) 510 zaqer=0.5/zaqe 511 zbqe=zeps1*zeps6+zeps3*zeps7-zeps2 512 zcqe=zeps2*stbl(ji,jj) 513 zdis=zbqe*zbqe-4.0*zaqe*zcqe 514 ! Presumably zdis can never be negative because gammas is very small compared to gammat 515 zsfrz=(-zbqe-SQRT(zdis))*zaqer 516 IF (zsfrz .lt. 0.0) zsfrz=(-zbqe+SQRT(zdis))*zaqer 517 zfrz(ji,jj)=zeps4+zlamb1*zsfrz 518 519 ! zfwflx is upward water flux 520 zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / zsfrz ) 521 ! zhtflx is upward heat flux (out of ocean) 522 ! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value 523 zhtflx = ( zgammat*rau0 - zcfac*zfwflx ) * rcp * (zti(ji,jj) - zfrz(ji,jj) ) 524 ! zwflx is upward water flux 525 ! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz 526 zfwflx = ( zgammas*rau0 - zcfac*zfwflx ) * (zsfrz - stbl(ji,jj)) / stbl(ji,jj) 527 ! test convergence and compute gammat 528 IF (( zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE. 529 530 nit = nit + 1 531 IF (nit .GE. 51) THEN 532 WRITE(numout,*) "sbcisf : too many iteration ... ", & 533 & zhtflx, zhtflx_b, zgammat, zgammas, nn_gammablk, ji, jj, mikt(ji,jj), narea 534 CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 535 END IF 536 ! save gammat and compute zhtflx_b 537 zgammat2d(ji,jj)=zgammat 538 zgammas2d(ji,jj)=zgammas 539 zhtflx_b = zhtflx 540 541 END DO 542 ! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value 543 qisf(ji,jj) = - zhtflx 544 ! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz 545 fwfisf(ji,jj) = zfwflx 546 ELSE 547 fwfisf(ji,jj) = 0._wp 548 qisf(ji,jj) = 0._wp 549 ENDIF 550 ! 551 END DO 552 END DO 553 ENDIF 554 ! lbclnk 555 CALL lbc_lnk(zgammas2d(:,:),'T',1.) 556 CALL lbc_lnk(zgammat2d(:,:),'T',1.) 557 ! output 558 CALL iom_put('isfgammat', zgammat2d) 559 CALL iom_put('isfgammas', zgammas2d) 560 ! 561 CALL wrk_dealloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d ) 527 ! 528 CALL iom_put('isfgammat', zgammat) 529 CALL iom_put('isfgammas', zgammas) 530 ! 531 CALL wrk_dealloc( jpi,jpj, zfrz , zgammat, zgammas ) 532 CALL wrk_dealloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b ) 562 533 ! 563 534 IF( nn_timing == 1 ) CALL timing_stop('sbc_isf_cav') … … 565 536 END SUBROUTINE sbc_isf_cav 566 537 567 568 SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf, ji, jj, lit ) 538 SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf ) 569 539 !!---------------------------------------------------------------------- 570 540 !! ** Purpose : compute the coefficient echange for heat flux … … 575 545 !! Jenkins et al., 2010, JPO, p2298-2312 576 546 !!--------------------------------------------------------------------- 577 REAL(wp), INTENT(inout) :: gt, gs, zqhisf, zqwisf578 INTEGER , INTENT(in) :: ji,jj579 LOGICAL , INTENT(inout) :: lit580 581 INTEGER :: ikt! loop index582 REAL(wp) :: zut, zvt,zustar ! U, V at T point and friction velocity547 REAL(wp), DIMENSION(:,:), INTENT(out) :: pgt, pgs 548 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pqhisf, pqwisf 549 ! 550 INTEGER :: ikt 551 INTEGER :: ji, jj ! loop index 552 REAL(wp), DIMENSION(:,:), POINTER :: zustar ! U, V at T point and friction velocity 583 553 REAL(wp) :: zdku, zdkv ! U, V shear 584 554 REAL(wp) :: zPr, zSc, zRc ! Prandtl, Scmidth and Richardson number … … 590 560 REAL(wp) :: zcoef ! temporary coef 591 561 REAL(wp) :: zdep 592 REAL(wp), PARAMETER :: zxsiN = 0.052 ! dimensionless constant 593 REAL(wp), PARAMETER :: epsln = 1.0e-20 ! a small positive number 594 REAL(wp), PARAMETER :: znu = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1) 595 REAL(wp) :: rcs = 1.0e-3_wp ! conversion: mm/s ==> m/s 562 REAL(wp) :: zeps = 1.0e-20_wp 563 REAL(wp), PARAMETER :: zxsiN = 0.052_wp ! dimensionless constant 564 REAL(wp), PARAMETER :: znu = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1) 596 565 REAL(wp), DIMENSION(2) :: zts, zab 597 566 !!--------------------------------------------------------------------- 598 ! 599 IF( nn_gammablk == 0 ) THEN 600 !! gamma is constant (specified in namelist) 601 gt = rn_gammat0 602 gs = rn_gammas0 603 lit = .FALSE. 604 ELSE IF ( nn_gammablk == 1 ) THEN 605 !! gamma is assume to be proportional to u* 606 !! WARNING in case of Losh 2008 tbl parametrization, 607 !! you have to used the mean value of u in the boundary layer) 608 !! not yet coded 609 !! Jenkins et al., 2010, JPO, p2298-2312 610 ikt = mikt(ji,jj) 611 !! Compute U and V at T points 612 ! zut = 0.5 * ( utbl(ji-1,jj ) + utbl(ji,jj) ) 613 ! zvt = 0.5 * ( vtbl(ji ,jj-1) + vtbl(ji,jj) ) 614 zut = utbl(ji,jj) 615 zvt = vtbl(ji,jj) 616 617 !! compute ustar 618 zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) ) 619 !! Compute mean value over the TBL 620 621 !! Compute gammats 622 gt = zustar * rn_gammat0 623 gs = zustar * rn_gammas0 624 lit = .FALSE. 625 ELSE IF ( nn_gammablk == 2 ) THEN 626 !! gamma depends of stability of boundary layer 627 !! WARNING in case of Losh 2008 tbl parametrization, 628 !! you have to used the mean value of u in the boundary layer) 629 !! not yet coded 630 !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 631 !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO) 567 CALL wrk_alloc( jpi,jpj, zustar ) 568 ! 569 SELECT CASE ( nn_gammablk ) 570 CASE ( 0 ) ! gamma is constant (specified in namelist) 571 !! ISOMIP formulation (Hunter et al, 2006) 572 pgt(:,:) = rn_gammat0 573 pgs(:,:) = rn_gammas0 574 575 CASE ( 1 ) ! gamma is assume to be proportional to u* 576 !! Jenkins et al., 2010, JPO, p2298-2312 577 !! Adopted by Asay-Davis et al. (2015) 578 579 !! compute ustar (eq. 24) 580 zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) ) 581 582 !! Compute gammats 583 pgt(:,:) = zustar(:,:) * rn_gammat0 584 pgs(:,:) = zustar(:,:) * rn_gammas0 585 586 CASE ( 2 ) ! gamma depends of stability of boundary layer 587 !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 588 !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO) 589 !! compute ustar 590 zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) ) 591 592 !! compute Pr and Sc number (can be improved) 593 zPr = 13.8_wp 594 zSc = 2432.0_wp 595 596 !! compute gamma mole 597 zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp 598 zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp 599 600 !! compute gamma 601 DO ji=2,jpi 602 DO jj=2,jpj 632 603 ikt = mikt(ji,jj) 633 604 634 !! Compute U and V at T points 635 zut = 0.5 * ( utbl(ji-1,jj ) + utbl(ji,jj) ) 636 zvt = 0.5 * ( vtbl(ji ,jj-1) + vtbl(ji,jj) ) 637 638 !! compute ustar 639 zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) ) 640 IF (zustar == 0._wp) THEN ! only for kt = 1 I think 641 gt = rn_gammat0 642 gs = rn_gammas0 605 IF (zustar(ji,jj) == 0._wp) THEN ! only for kt = 1 I think 606 pgt = rn_gammat0 607 pgs = rn_gammas0 643 608 ELSE 644 !! compute Rc number (as done in zdfric.F90) 645 zcoef = 0.5 / fse3w(ji,jj,ikt) 646 ! ! shear of horizontal velocity 647 zdku = zcoef * ( un(ji-1,jj ,ikt ) + un(ji,jj,ikt ) & 648 & -un(ji-1,jj ,ikt+1) - un(ji,jj,ikt+1) ) 649 zdkv = zcoef * ( vn(ji ,jj-1,ikt ) + vn(ji,jj,ikt ) & 650 & -vn(ji ,jj-1,ikt+1) - vn(ji,jj,ikt+1) ) 651 ! ! richardson number (minimum value set to zero) 652 zRc = rn2(ji,jj,ikt+1) / ( zdku*zdku + zdkv*zdkv + 1.e-20 ) 653 654 !! compute Pr and Sc number (can be improved) 655 zPr = 13.8 656 zSc = 2432.0 657 658 !! compute gamma mole 659 zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0 660 zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0 661 662 !! compute bouyancy 663 zts(jp_tem) = ttbl(ji,jj) 664 zts(jp_sal) = stbl(ji,jj) 665 zdep = fsdepw(ji,jj,ikt) 666 ! 667 CALL eos_rab( zts, zdep, zab ) 609 !! compute Rc number (as done in zdfric.F90) 610 zcoef = 0.5_wp / e3w_n(ji,jj,ikt) 611 ! ! shear of horizontal velocity 612 zdku = zcoef * ( un(ji-1,jj ,ikt ) + un(ji,jj,ikt ) & 613 & -un(ji-1,jj ,ikt+1) - un(ji,jj,ikt+1) ) 614 zdkv = zcoef * ( vn(ji ,jj-1,ikt ) + vn(ji,jj,ikt ) & 615 & -vn(ji ,jj-1,ikt+1) - vn(ji,jj,ikt+1) ) 616 ! ! richardson number (minimum value set to zero) 617 zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps ) 618 619 !! compute bouyancy 620 zts(jp_tem) = ttbl(ji,jj) 621 zts(jp_sal) = stbl(ji,jj) 622 zdep = gdepw_n(ji,jj,ikt) 668 623 ! 669 !! compute length scale 670 zbuofdep = grav * ( zab(jp_tem) * zqhisf - zab(jp_sal) * zqwisf ) !!!!!!!!!!!!!!!!!!!!!!!!!!!! 671 672 !! compute Monin Obukov Length 673 ! Maximum boundary layer depth 674 zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) -0.001 675 ! Compute Monin obukhov length scale at the surface and Ekman depth: 676 zmob = zustar ** 3 / (vkarmn * (zbuofdep + epsln)) 677 zmols = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt) 678 679 !! compute eta* (stability parameter) 680 zetastar = 1 / ( SQRT(1 + MAX(zxsiN * zustar / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0))) 681 682 !! compute the sublayer thickness 683 zhnu = 5 * znu / zustar 684 !! compute gamma turb 685 zgturb = 1/vkarmn * LOG(zustar * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) & 686 & + 1 / ( 2 * zxsiN * zetastar ) - 1/vkarmn 687 688 !! compute gammats 689 gt = zustar / (zgturb + zgmolet) 690 gs = zustar / (zgturb + zgmoles) 624 CALL eos_rab( zts, zdep, zab ) 625 ! 626 !! compute length scale 627 zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) ) !!!!!!!!!!!!!!!!!!!!!!!!!!!! 628 629 !! compute Monin Obukov Length 630 ! Maximum boundary layer depth 631 zhmax = gdept_n(ji,jj,mbkt(ji,jj)) - gdepw_n(ji,jj,mikt(ji,jj)) -0.001_wp 632 ! Compute Monin obukhov length scale at the surface and Ekman depth: 633 zmob = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps)) 634 zmols = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt) 635 636 !! compute eta* (stability parameter) 637 zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0_wp))) 638 639 !! compute the sublayer thickness 640 zhnu = 5 * znu / zustar(ji,jj) 641 642 !! compute gamma turb 643 zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) & 644 & + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn 645 646 !! compute gammats 647 pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet) 648 pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles) 691 649 END IF 692 END IF 693 ! 694 END SUBROUTINE 695 696 697 SUBROUTINE sbc_isf_tbl( varin, varout, cptin ) 650 END DO 651 END DO 652 CALL lbc_lnk(pgt(:,:),'T',1.) 653 CALL lbc_lnk(pgs(:,:),'T',1.) 654 END SELECT 655 CALL wrk_dealloc( jpi,jpj, zustar ) 656 ! 657 END SUBROUTINE sbc_isf_gammats 658 659 SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin ) 698 660 !!---------------------------------------------------------------------- 699 661 !! *** SUBROUTINE sbc_isf_tbl *** 700 662 !! 701 !! ** Purpose : compute mean T/S/U/V in the boundary layer 702 !! 703 !!---------------------------------------------------------------------- 704 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: varin 705 REAL(wp), DIMENSION(:,:) , INTENT(out):: varout 663 !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point 664 !! 665 !!---------------------------------------------------------------------- 666 REAL(wp), DIMENSION(:,:,:), INTENT( in ) :: pvarin 667 REAL(wp), DIMENSION(:,:) , INTENT( out ) :: pvarout 668 CHARACTER(len=1), INTENT( in ) :: cd_ptin ! point of variable in/out 669 ! 670 REAL(wp) :: ze3, zhk 671 REAL(wp), DIMENSION(:,:), POINTER :: zhisf_tbl ! thickness of the tbl 672 673 INTEGER :: ji, jj, jk ! loop index 674 INTEGER :: ikt, ikb ! top and bottom index of the tbl 675 !!---------------------------------------------------------------------- 676 ! allocation 677 CALL wrk_alloc( jpi,jpj, zhisf_tbl) 706 678 707 CHARACTER(len=1), INTENT(in) :: cptin ! point of variable in/out 708 709 REAL(wp) :: ze3, zhk 710 REAL(wp), DIMENSION(:,:), POINTER :: zikt 711 712 INTEGER :: ji,jj,jk 713 INTEGER :: ikt,ikb 714 INTEGER, DIMENSION(:,:), POINTER :: mkt, mkb 715 716 CALL wrk_alloc( jpi,jpj, mkt, mkb ) 717 CALL wrk_alloc( jpi,jpj, zikt ) 718 719 ! get first and last level of tbl 720 mkt(:,:) = misfkt(:,:) 721 mkb(:,:) = misfkb(:,:) 722 723 varout(:,:)=0._wp 724 DO jj = 2,jpj 725 DO ji = 2,jpi 726 IF (ssmask(ji,jj) == 1) THEN 727 ikt = mkt(ji,jj) 728 ikb = mkb(ji,jj) 679 ! initialisation 680 pvarout(:,:)=0._wp 681 682 SELECT CASE ( cd_ptin ) 683 CASE ( 'U' ) ! compute U in the top boundary layer at T- point 684 DO jj = 1,jpj 685 DO ji = 1,jpi 686 ikt = miku(ji,jj) ; ikb = miku(ji,jj) 687 ! thickness of boundary layer at least the top level thickness 688 zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3u_n(ji,jj,ikt)) 689 690 ! determine the deepest level influenced by the boundary layer 691 DO jk = ikt+1, mbku(ji,jj) 692 IF ( (SUM(e3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk 693 END DO 694 zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 729 695 730 696 ! level fully include in the ice shelf boundary layer 731 697 DO jk = ikt, ikb - 1 732 ze3 = fse3t_n(ji,jj,jk) 733 IF (cptin == 'T' ) varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 734 IF (cptin == 'U' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji-1,jj,jk)) & 735 & * r1_hisf_tbl(ji,jj) * ze3 736 IF (cptin == 'V' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji,jj-1,jk)) & 737 & * r1_hisf_tbl(ji,jj) * ze3 698 ze3 = e3u_n(ji,jj,jk) 699 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 738 700 END DO 739 701 740 702 ! level partially include in ice shelf boundary layer 741 zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) 742 IF (cptin == 'T') & 743 & varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk) 744 IF (cptin == 'U') & 745 & varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk) 746 IF (cptin == 'V') & 747 & varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk) 748 END IF 749 END DO 750 END DO 751 752 CALL wrk_dealloc( jpi,jpj, mkt, mkb ) 753 CALL wrk_dealloc( jpi,jpj, zikt ) 754 755 IF (cptin == 'T') CALL lbc_lnk(varout,'T',1.) 756 IF (cptin == 'U' .OR. cptin == 'V') CALL lbc_lnk(varout,'T',-1.) 703 zhk = SUM( e3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 704 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 705 END DO 706 END DO 707 DO jj = 2,jpj 708 DO ji = 2,jpi 709 pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj)) 710 END DO 711 END DO 712 CALL lbc_lnk(pvarout,'T',-1.) 713 714 CASE ( 'V' ) ! compute V in the top boundary layer at T- point 715 DO jj = 1,jpj 716 DO ji = 1,jpi 717 ikt = mikv(ji,jj) ; ikb = mikv(ji,jj) 718 ! thickness of boundary layer at least the top level thickness 719 zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3v_n(ji,jj,ikt)) 720 721 ! determine the deepest level influenced by the boundary layer 722 DO jk = ikt+1, mbkv(ji,jj) 723 IF ( (SUM(e3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk 724 END DO 725 zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 726 727 ! level fully include in the ice shelf boundary layer 728 DO jk = ikt, ikb - 1 729 ze3 = e3v_n(ji,jj,jk) 730 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 731 END DO 732 733 ! level partially include in ice shelf boundary layer 734 zhk = SUM( e3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 735 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 736 END DO 737 END DO 738 DO jj = 2,jpj 739 DO ji = 2,jpi 740 pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1)) 741 END DO 742 END DO 743 CALL lbc_lnk(pvarout,'T',-1.) 744 745 CASE ( 'T' ) ! compute T in the top boundary layer at T- point 746 DO jj = 1,jpj 747 DO ji = 1,jpi 748 ikt = misfkt(ji,jj) 749 ikb = misfkb(ji,jj) 750 751 ! level fully include in the ice shelf boundary layer 752 DO jk = ikt, ikb - 1 753 ze3 = e3t_n(ji,jj,jk) 754 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 755 END DO 756 757 ! level partially include in ice shelf boundary layer 758 zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) 759 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 760 END DO 761 END DO 762 END SELECT 763 764 ! mask mean tbl value 765 pvarout(:,:) = pvarout(:,:) * ssmask(:,:) 766 767 ! deallocation 768 CALL wrk_dealloc( jpi,jpj, zhisf_tbl ) 757 769 ! 758 770 END SUBROUTINE sbc_isf_tbl … … 771 783 !! ** Action : phdivn decreased by the runoff inflow 772 784 !!---------------------------------------------------------------------- 773 REAL(wp), DIMENSION(:,:,:), INTENT( inout) :: phdivn ! horizontal divergence774 ! !785 REAL(wp), DIMENSION(:,:,:), INTENT( inout ) :: phdivn ! horizontal divergence 786 ! 775 787 INTEGER :: ji, jj, jk ! dummy loop indices 776 788 INTEGER :: ikt, ikb 777 INTEGER :: nk_isf 778 REAL(wp) :: zhk, z1_hisf_tbl, zhisf_tbl 779 REAL(wp) :: zfact ! local scalar 789 REAL(wp) :: zhk 790 REAL(wp) :: zfact ! local scalar 780 791 !!---------------------------------------------------------------------- 781 792 ! 782 793 zfact = 0.5_wp 783 794 ! 784 IF (lk_vvl) THEN ! need to re compute level distribution of isf fresh water795 IF(.NOT.ln_linssh ) THEN ! need to re compute level distribution of isf fresh water 785 796 DO jj = 1,jpj 786 797 DO ji = 1,jpi … … 788 799 ikb = misfkt(ji,jj) 789 800 ! thickness of boundary layer at least the top level thickness 790 rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t(ji,jj,ikt))801 rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt)) 791 802 792 803 ! determine the deepest level influenced by the boundary layer 793 ! test on tmask useless ?????794 804 DO jk = ikt, mbkt(ji,jj) 795 IF ( (SUM( fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk796 END DO 797 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM( fse3t(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.805 IF ( (SUM(e3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 806 END DO 807 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 798 808 misfkb(ji,jj) = ikb ! last wet level of the tbl 799 809 r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 800 810 801 zhk = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 802 ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer 803 END DO 804 END DO 805 END IF ! vvl case 806 ! 811 zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 812 ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer 813 END DO 814 END DO 815 END IF 816 ! 817 !== ice shelf melting distributed over several levels ==! 807 818 DO jj = 1,jpj 808 819 DO ji = 1,jpi … … 812 823 DO jk = ikt, ikb - 1 813 824 phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) & 814 & 825 & * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact 815 826 END DO 816 827 ! level partially include in ice shelf boundary layer 817 828 phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) & 818 & + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 819 !== ice shelf melting mass distributed over several levels ==! 829 & + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 820 830 END DO 821 831 END DO 822 832 ! 823 833 END SUBROUTINE sbc_isf_div 824 825 826 FUNCTION tinsitu( ptem, psal, ppress ) RESULT( pti ) 827 !!---------------------------------------------------------------------- 828 !! *** ROUTINE eos_init *** 829 !! 830 !! ** Purpose : Compute the in-situ temperature [Celcius] 831 !! 832 !! ** Method : 833 !! 834 !! Reference : Bryden,h.,1973,deep-sea res.,20,401-408 835 !!---------------------------------------------------------------------- 836 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ptem ! potential temperature [Celcius] 837 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu] 838 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ppress ! pressure [dBar] 839 REAL(wp), DIMENSION(:,:), POINTER :: pti ! in-situ temperature [Celcius] 840 ! REAL(wp) :: fsatg 841 ! REAL(wp) :: pfps, pfpt, pfphp 842 REAL(wp) :: zt, zs, zp, zh, zq, zxk 843 INTEGER :: ji, jj ! dummy loop indices 844 ! 845 CALL wrk_alloc( jpi,jpj, pti ) 846 ! 847 DO jj=1,jpj 848 DO ji=1,jpi 849 zh = ppress(ji,jj) 850 ! Theta1 851 zt = ptem(ji,jj) 852 zs = psal(ji,jj) 853 zp = 0.0 854 zxk= zh * fsatg( zs, zt, zp ) 855 zt = zt + 0.5 * zxk 856 zq = zxk 857 ! Theta2 858 zp = zp + 0.5 * zh 859 zxk= zh*fsatg( zs, zt, zp ) 860 zt = zt + 0.29289322 * ( zxk - zq ) 861 zq = 0.58578644 * zxk + 0.121320344 * zq 862 ! Theta3 863 zxk= zh * fsatg( zs, zt, zp ) 864 zt = zt + 1.707106781 * ( zxk - zq ) 865 zq = 3.414213562 * zxk - 4.121320344 * zq 866 ! Theta4 867 zp = zp + 0.5 * zh 868 zxk= zh * fsatg( zs, zt, zp ) 869 pti(ji,jj) = zt + ( zxk - 2.0 * zq ) / 6.0 870 END DO 871 END DO 872 ! 873 CALL wrk_dealloc( jpi,jpj, pti ) 874 ! 875 END FUNCTION tinsitu 876 877 878 FUNCTION fsatg( pfps, pfpt, pfphp ) 879 !!---------------------------------------------------------------------- 880 !! *** FUNCTION fsatg *** 881 !! 882 !! ** Purpose : Compute the Adiabatic laspse rate [Celcius].[decibar]^-1 883 !! 884 !! ** Reference : Bryden,h.,1973,deep-sea res.,20,401-408 885 !! 886 !! ** units : pressure pfphp decibars 887 !! temperature pfpt deg celsius (ipts-68) 888 !! salinity pfps (ipss-78) 889 !! adiabatic fsatg deg. c/decibar 890 !!---------------------------------------------------------------------- 891 REAL(wp) :: pfps, pfpt, pfphp 892 REAL(wp) :: fsatg 893 ! 894 fsatg = (((-2.1687e-16*pfpt+1.8676e-14)*pfpt-4.6206e-13)*pfphp & 895 & +((2.7759e-12*pfpt-1.1351e-10)*(pfps-35.)+((-5.4481e-14*pfpt & 896 & +8.733e-12)*pfpt-6.7795e-10)*pfpt+1.8741e-8))*pfphp & 897 & +(-4.2393e-8*pfpt+1.8932e-6)*(pfps-35.) & 898 & +((6.6228e-10*pfpt-6.836e-8)*pfpt+8.5258e-6)*pfpt+3.5803e-5 899 ! 900 END FUNCTION fsatg 901 !!====================================================================== 834 !!====================================================================== 902 835 END MODULE sbcisf -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
r5836 r6140 17 17 18 18 !!---------------------------------------------------------------------- 19 !! sbc_init 20 !! sbc 19 !! sbc_init : read namsbc namelist 20 !! sbc : surface ocean momentum, heat and freshwater boundary conditions 21 21 !!---------------------------------------------------------------------- 22 USE oce ! ocean dynamics and tracers 23 USE dom_oce ! ocean space and time domain 24 USE phycst ! physical constants 25 USE sbc_oce ! Surface boundary condition: ocean fields 26 USE trc_oce ! shared ocean-passive tracers variables 27 USE sbc_ice ! Surface boundary condition: ice fields 28 USE sbcdcy ! surface boundary condition: diurnal cycle 29 USE sbcssm ! surface boundary condition: sea-surface mean variables 30 USE sbcana ! surface boundary condition: analytical formulation 31 USE sbcflx ! surface boundary condition: flux formulation 32 USE sbcblk_clio ! surface boundary condition: bulk formulation : CLIO 33 USE sbcblk_core ! surface boundary condition: bulk formulation : CORE 34 USE sbcblk_mfs ! surface boundary condition: bulk formulation : MFS 35 USE sbcice_if ! surface boundary condition: ice-if sea-ice model 36 USE sbcice_lim ! surface boundary condition: LIM 3.0 sea-ice model 37 USE sbcice_lim_2 ! surface boundary condition: LIM 2.0 sea-ice model 38 USE sbcice_cice ! surface boundary condition: CICE sea-ice model 39 USE sbccpl ! surface boundary condition: coupled florulation 40 USE cpl_oasis3 ! OASIS routines for coupling 41 USE sbcssr ! surface boundary condition: sea surface restoring 42 USE sbcrnf ! surface boundary condition: runoffs 43 USE sbcisf ! surface boundary condition: ice shelf 44 USE sbcfwb ! surface boundary condition: freshwater budget 45 USE closea ! closed sea 46 USE icbstp ! Icebergs! 47 48 USE prtctl ! Print control (prt_ctl routine) 49 USE iom ! IOM library 50 USE in_out_manager ! I/O manager 51 USE lib_mpp ! MPP library 52 USE timing ! Timing 53 USE sbcwave ! Wave module 54 USE bdy_par ! Require lk_bdy 22 USE oce ! ocean dynamics and tracers 23 USE dom_oce ! ocean space and time domain 24 USE phycst ! physical constants 25 USE sbc_oce ! Surface boundary condition: ocean fields 26 USE trc_oce ! shared ocean-passive tracers variables 27 USE sbc_ice ! Surface boundary condition: ice fields 28 USE sbcdcy ! surface boundary condition: diurnal cycle 29 USE sbcssm ! surface boundary condition: sea-surface mean variables 30 USE sbcana ! surface boundary condition: analytical formulation 31 USE sbcflx ! surface boundary condition: flux formulation 32 USE sbcblk_clio ! surface boundary condition: bulk formulation : CLIO 33 USE sbcblk_core ! surface boundary condition: bulk formulation : CORE 34 USE sbcblk_mfs ! surface boundary condition: bulk formulation : MFS 35 USE sbcice_if ! surface boundary condition: ice-if sea-ice model 36 USE sbcice_lim ! surface boundary condition: LIM 3.0 sea-ice model 37 USE sbcice_lim_2 ! surface boundary condition: LIM 2.0 sea-ice model 38 USE sbcice_cice ! surface boundary condition: CICE sea-ice model 39 USE sbccpl ! surface boundary condition: coupled florulation 40 USE cpl_oasis3 ! OASIS routines for coupling 41 USE sbcssr ! surface boundary condition: sea surface restoring 42 USE sbcrnf ! surface boundary condition: runoffs 43 USE sbcisf ! surface boundary condition: ice shelf 44 USE sbcfwb ! surface boundary condition: freshwater budget 45 USE closea ! closed sea 46 USE icbstp ! Icebergs 47 USE traqsr ! active tracers: light penetration 48 USE sbcwave ! Wave module 49 USE bdy_par ! Require lk_bdy 50 ! 51 USE prtctl ! Print control (prt_ctl routine) 52 USE iom ! IOM library 53 USE in_out_manager ! I/O manager 54 USE lib_mpp ! MPP library 55 USE timing ! Timing 56 57 USE diurnal_bulk, ONLY: & 58 & ln_diurnal_only 55 59 56 60 IMPLICIT NONE … … 62 66 INTEGER :: nsbc ! type of surface boundary condition (deduced from namsbc informations) 63 67 64 !! * Substitutions65 # include "domzgr_substitute.h90"66 68 !!---------------------------------------------------------------------- 67 69 !! NEMO/OPA 4.0 , NEMO-consortium (2011) … … 85 87 INTEGER :: icpt ! local integer 86 88 !! 87 NAMELIST/namsbc/ nn_fsbc , ln_ana , ln_flx, ln_blk_clio, ln_blk_core, ln_mixcpl, & 88 & ln_blk_mfs, ln_apr_dyn, nn_ice, nn_ice_embd, ln_dm2dc , ln_rnf , & 89 & ln_ssr , nn_isf , nn_fwb, ln_cdgw , ln_wave , ln_sdw , & 90 & nn_lsm , nn_limflx , nn_components, ln_cpl 89 NAMELIST/namsbc/ nn_fsbc , ln_ana , ln_flx, ln_blk_clio, ln_blk_core, ln_blk_mfs, & 90 & ln_cpl , ln_mixcpl, nn_components , nn_limflx , & 91 & ln_traqsr, ln_dm2dc , & 92 & nn_ice , nn_ice_embd, & 93 & ln_rnf , ln_ssr , ln_isf , nn_fwb , ln_apr_dyn, & 94 & ln_wave , & 95 & nn_lsm 91 96 INTEGER :: ios 92 97 INTEGER :: ierr, ierr0, ierr1, ierr2, ierr3, jpm 93 98 LOGICAL :: ll_purecpl 94 99 !!---------------------------------------------------------------------- 95 100 ! 96 101 IF(lwp) THEN 97 102 WRITE(numout,*) … … 99 104 WRITE(numout,*) '~~~~~~~~ ' 100 105 ENDIF 101 106 ! 102 107 REWIND( numnam_ref ) ! Namelist namsbc in reference namelist : Surface boundary 103 108 READ ( numnam_ref, namsbc, IOSTAT = ios, ERR = 901) 104 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc in reference namelist', lwp )105 109 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc in reference namelist', lwp ) 110 ! 106 111 REWIND( numnam_cfg ) ! Namelist namsbc in configuration namelist : Parameters of the run 107 112 READ ( numnam_cfg, namsbc, IOSTAT = ios, ERR = 902 ) 108 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc in configuration namelist', lwp )109 IF(lwm) WRITE 110 113 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc in configuration namelist', lwp ) 114 IF(lwm) WRITE( numond, namsbc ) 115 ! 111 116 ! ! overwrite namelist parameter using CPP key information 112 117 IF( Agrif_Root() ) THEN ! AGRIF zoom … … 119 124 nn_ice = 0 120 125 ENDIF 121 126 ! 122 127 IF(lwp) THEN ! Control print 123 128 WRITE(numout,*) ' Namelist namsbc (partly overwritten with CPP key setting)' 124 129 WRITE(numout,*) ' frequency update of sbc (and ice) nn_fsbc = ', nn_fsbc 125 WRITE(numout,*) ' Type of sbc : ' 126 WRITE(numout,*) ' analytical formulation ln_ana = ', ln_ana 127 WRITE(numout,*) ' flux formulation ln_flx = ', ln_flx 128 WRITE(numout,*) ' CLIO bulk formulation ln_blk_clio = ', ln_blk_clio 129 WRITE(numout,*) ' CORE bulk formulation ln_blk_core = ', ln_blk_core 130 WRITE(numout,*) ' MFS bulk formulation ln_blk_mfs = ', ln_blk_mfs 131 WRITE(numout,*) ' ocean-atmosphere coupled formulation ln_cpl = ', ln_cpl 132 WRITE(numout,*) ' forced-coupled mixed formulation ln_mixcpl = ', ln_mixcpl 133 WRITE(numout,*) ' OASIS coupling (with atm or sas) lk_oasis = ', lk_oasis 134 WRITE(numout,*) ' components of your executable nn_components = ', nn_components 135 WRITE(numout,*) ' Multicategory heat flux formulation (LIM3) nn_limflx = ', nn_limflx 130 WRITE(numout,*) ' Type of air-sea fluxes : ' 131 WRITE(numout,*) ' analytical formulation ln_ana = ', ln_ana 132 WRITE(numout,*) ' flux formulation ln_flx = ', ln_flx 133 WRITE(numout,*) ' CLIO bulk formulation ln_blk_clio = ', ln_blk_clio 134 WRITE(numout,*) ' CORE bulk formulation ln_blk_core = ', ln_blk_core 135 WRITE(numout,*) ' MFS bulk formulation ln_blk_mfs = ', ln_blk_mfs 136 WRITE(numout,*) ' Type of coupling (Ocean/Ice/Atmosphere) : ' 137 WRITE(numout,*) ' ocean-atmosphere coupled formulation ln_cpl = ', ln_cpl 138 WRITE(numout,*) ' forced-coupled mixed formulation ln_mixcpl = ', ln_mixcpl 139 WRITE(numout,*) ' OASIS coupling (with atm or sas) lk_oasis = ', lk_oasis 140 WRITE(numout,*) ' components of your executable nn_components = ', nn_components 141 WRITE(numout,*) ' Multicategory heat flux formulation (LIM3) nn_limflx = ', nn_limflx 142 WRITE(numout,*) ' Sea-ice : ' 143 WRITE(numout,*) ' ice management in the sbc (=0/1/2/3) nn_ice = ', nn_ice 144 WRITE(numout,*) ' ice-ocean embedded/levitating (=0/1/2) nn_ice_embd = ', nn_ice_embd 136 145 WRITE(numout,*) ' Misc. options of sbc : ' 137 WRITE(numout,*) ' Patm gradient added in ocean & ice Eqs. ln_apr_dyn = ', ln_apr_dyn 138 WRITE(numout,*) ' ice management in the sbc (=0/1/2/3) nn_ice = ', nn_ice 139 WRITE(numout,*) ' ice-ocean embedded/levitating (=0/1/2) nn_ice_embd = ', nn_ice_embd 140 WRITE(numout,*) ' daily mean to diurnal cycle qsr ln_dm2dc = ', ln_dm2dc 141 WRITE(numout,*) ' runoff / runoff mouths ln_rnf = ', ln_rnf 142 WRITE(numout,*) ' iceshelf formulation nn_isf = ', nn_isf 143 WRITE(numout,*) ' Sea Surface Restoring on SST and/or SSS ln_ssr = ', ln_ssr 144 WRITE(numout,*) ' FreshWater Budget control (=0/1/2) nn_fwb = ', nn_fwb 145 WRITE(numout,*) ' closed sea (=0/1) (set in namdom) nn_closea = ', nn_closea 146 WRITE(numout,*) ' n. of iterations if land-sea-mask applied nn_lsm = ', nn_lsm 147 ENDIF 148 149 ! LIM3 Multi-category heat flux formulation 150 SELECT CASE ( nn_limflx) 151 CASE ( -1 ) 152 IF(lwp) WRITE(numout,*) ' Use of per-category fluxes (nn_limflx = -1) ' 153 CASE ( 0 ) 154 IF(lwp) WRITE(numout,*) ' Average per-category fluxes (nn_limflx = 0) ' 155 CASE ( 1 ) 156 IF(lwp) WRITE(numout,*) ' Average then redistribute per-category fluxes (nn_limflx = 1) ' 157 CASE ( 2 ) 158 IF(lwp) WRITE(numout,*) ' Redistribute a single flux over categories (nn_limflx = 2) ' 159 END SELECT 160 ! 161 IF ( nn_components /= jp_iam_nemo .AND. .NOT. lk_oasis ) & 162 & CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but key_oasis3 disabled' ) 163 IF ( nn_components == jp_iam_opa .AND. ln_cpl ) & 164 & CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_cpl = T in OPA' ) 165 IF ( nn_components == jp_iam_opa .AND. ln_mixcpl ) & 166 & CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl = T in OPA' ) 167 IF ( ln_cpl .AND. .NOT. lk_oasis ) & 168 & CALL ctl_stop( 'STOP', 'sbc_init : OASIS-coupled atmosphere model, but key_oasis3 disabled' ) 146 WRITE(numout,*) ' Light penetration in temperature Eq. ln_traqsr = ', ln_traqsr 147 WRITE(numout,*) ' daily mean to diurnal cycle qsr ln_dm2dc = ', ln_dm2dc 148 WRITE(numout,*) ' Sea Surface Restoring on SST and/or SSS ln_ssr = ', ln_ssr 149 WRITE(numout,*) ' FreshWater Budget control (=0/1/2) nn_fwb = ', nn_fwb 150 WRITE(numout,*) ' Patm gradient added in ocean & ice Eqs. ln_apr_dyn = ', ln_apr_dyn 151 WRITE(numout,*) ' runoff / runoff mouths ln_rnf = ', ln_rnf 152 WRITE(numout,*) ' iceshelf formulation ln_isf = ', ln_isf 153 WRITE(numout,*) ' closed sea (=0/1) (set in namdom) nn_closea = ', nn_closea 154 WRITE(numout,*) ' nb of iterations if land-sea-mask applied nn_lsm = ', nn_lsm 155 WRITE(numout,*) ' surface wave ln_wave = ', ln_wave 156 ENDIF 157 ! 158 IF(lwp) THEN 159 WRITE(numout,*) 160 SELECT CASE ( nn_limflx ) ! LIM3 Multi-category heat flux formulation 161 CASE ( -1 ) ; WRITE(numout,*) ' LIM3: use per-category fluxes (nn_limflx = -1) ' 162 CASE ( 0 ) ; WRITE(numout,*) ' LIM3: use average per-category fluxes (nn_limflx = 0) ' 163 CASE ( 1 ) ; WRITE(numout,*) ' LIM3: use average then redistribute per-category fluxes (nn_limflx = 1) ' 164 CASE ( 2 ) ; WRITE(numout,*) ' LIM3: Redistribute a single flux over categories (nn_limflx = 2) ' 165 END SELECT 166 ENDIF 167 ! 168 IF( nn_components /= jp_iam_nemo .AND. .NOT. lk_oasis ) & 169 & CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but key_oasis3 disabled' ) 170 IF( nn_components == jp_iam_opa .AND. ln_cpl ) & 171 & CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but ln_cpl = T in OPA' ) 172 IF( nn_components == jp_iam_opa .AND. ln_mixcpl ) & 173 & CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl = T in OPA' ) 174 IF( ln_cpl .AND. .NOT. lk_oasis ) & 175 & CALL ctl_stop( 'sbc_init : OASIS-coupled atmosphere model, but key_oasis3 disabled' ) 169 176 IF( ln_mixcpl .AND. .NOT. lk_oasis ) & 170 177 & CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) requires the cpp key key_oasis3' ) … … 178 185 179 186 ! ! Checks: 180 IF( nn_isf .EQ. 0) THEN ! variable initialisation if no ice shelf181 IF( sbc_isf_alloc() /= 0 ) CALL ctl_stop( ' sbc_init : unable to allocate sbc_isf arrays' )187 IF( .NOT. ln_isf ) THEN ! variable initialisation if no ice shelf 188 IF( sbc_isf_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_isf arrays' ) 182 189 fwfisf (:,:) = 0.0_wp ; fwfisf_b (:,:) = 0.0_wp 183 190 risf_tsc(:,:,:) = 0.0_wp ; risf_tsc_b(:,:,:) = 0.0_wp 184 rdivisf = 0.0_wp185 191 END IF 186 IF( nn_ice == 0 .AND. nn_components /= jp_iam_opa ) fr_i(:,:) = 0. e0! no ice in the domain, ice fraction is always zero187 188 sfx(:,:) = 0. 0_wp! the salt flux due to freezing/melting will be computed (i.e. will be non-zero)192 IF( nn_ice == 0 .AND. nn_components /= jp_iam_opa ) fr_i(:,:) = 0._wp ! no ice in the domain, ice fraction is always zero 193 194 sfx(:,:) = 0._wp ! the salt flux due to freezing/melting will be computed (i.e. will be non-zero) 189 195 ! only if sea-ice is present 190 196 191 fmmflx(:,:) = 0. 0_wp! freezing-melting array initialisation197 fmmflx(:,:) = 0._wp ! freezing-melting array initialisation 192 198 193 taum(:,:) = 0. 0_wp! Initialise taum for use in gls in case of reduced restart199 taum(:,:) = 0._wp ! Initialise taum for use in gls in case of reduced restart 194 200 195 201 ! ! restartability … … 214 220 & CALL ctl_stop( 'diurnal cycle into qsr field from daily values requires a flux or core-bulk formulation' ) 215 221 216 IF ( ln_wave ) THEN217 !Activated wave module but neither drag nor stokes drift activated218 IF ( .NOT.(ln_cdgw .OR. ln_sdw) ) THEN219 CALL ctl_warn( 'Ask for wave coupling but nor drag coefficient (ln_cdgw=F) neither stokes drift activated (ln_sdw=F)' )220 !drag coefficient read from wave model definable only with mfs bulk formulae and core221 ELSEIF (ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) ) THEN222 CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core')223 ENDIF224 ELSE225 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) ')228 ENDIF229 222 ! ! Choice of the Surface Boudary Condition (set nsbc) 230 223 ll_purecpl = ln_cpl .AND. .NOT. ln_mixcpl … … 245 238 IF(lwp) THEN 246 239 WRITE(numout,*) 247 IF( nsbc == jp_gyre ) WRITE(numout,*) ' GYRE analytical formulation' 248 IF( nsbc == jp_ana ) WRITE(numout,*) ' analytical formulation' 249 IF( nsbc == jp_flx ) WRITE(numout,*) ' flux formulation' 250 IF( nsbc == jp_clio ) WRITE(numout,*) ' CLIO bulk formulation' 251 IF( nsbc == jp_core ) WRITE(numout,*) ' CORE bulk formulation' 252 IF( nsbc == jp_purecpl ) WRITE(numout,*) ' pure coupled formulation' 253 IF( nsbc == jp_mfs ) WRITE(numout,*) ' MFS Bulk formulation' 254 IF( nsbc == jp_none ) WRITE(numout,*) ' OPA coupled to SAS via oasis' 255 IF( ln_mixcpl ) WRITE(numout,*) ' + forced-coupled mixed formulation' 240 SELECT CASE( nsbc ) 241 CASE( jp_gyre ) ; WRITE(numout,*) ' GYRE analytical formulation' 242 CASE( jp_ana ) ; WRITE(numout,*) ' analytical formulation' 243 CASE( jp_flx ) ; WRITE(numout,*) ' flux formulation' 244 CASE( jp_clio ) ; WRITE(numout,*) ' CLIO bulk formulation' 245 CASE( jp_core ) ; WRITE(numout,*) ' CORE bulk formulation' 246 CASE( jp_purecpl ) ; WRITE(numout,*) ' pure coupled formulation' 247 CASE( jp_mfs ) ; WRITE(numout,*) ' MFS Bulk formulation' 248 CASE( jp_none ) ; WRITE(numout,*) ' OPA coupled to SAS via oasis' 249 IF( ln_mixcpl ) WRITE(numout,*) ' + forced-coupled mixed formulation' 250 END SELECT 256 251 IF( nn_components/= jp_iam_nemo ) & 257 & WRITE(numout,*) '+ OASIS coupled SAS'252 & WRITE(numout,*) ' + OASIS coupled SAS' 258 253 ENDIF 259 254 ! 260 255 IF( lk_oasis ) CALL sbc_cpl_init (nn_ice) ! OASIS initialisation. must be done before: (1) first time step 261 256 ! ! (2) the use of nn_fsbc 262 263 ! nn_fsbc initialization if OPA-SAS coupling via OASIS 264 ! sas model time step has to be declared in OASIS (mandatory) -> nn_fsbc has to be modified accordingly 265 IF ( nn_components /= jp_iam_nemo ) THEN 266 IF ( nn_components == jp_iam_opa ) nn_fsbc = cpl_freq('O_SFLX') / NINT(rdt) 267 IF ( nn_components == jp_iam_sas ) nn_fsbc = cpl_freq('I_SFLX') / NINT(rdt) 257 ! nn_fsbc initialization if OPA-SAS coupling via OASIS 258 ! sas model time step has to be declared in OASIS (mandatory) -> nn_fsbc has to be modified accordingly 259 IF( nn_components /= jp_iam_nemo ) THEN 260 IF( nn_components == jp_iam_opa ) nn_fsbc = cpl_freq('O_SFLX') / NINT(rdt) 261 IF( nn_components == jp_iam_sas ) nn_fsbc = cpl_freq('I_SFLX') / NINT(rdt) 268 262 ! 269 263 IF(lwp)THEN … … 273 267 ENDIF 274 268 ENDIF 275 269 ! 276 270 IF( MOD( nitend - nit000 + 1, nn_fsbc) /= 0 .OR. & 277 271 MOD( nstock , nn_fsbc) /= 0 ) THEN … … 286 280 IF( ln_dm2dc .AND. ( ( NINT(rday) / ( nn_fsbc * NINT(rdt) ) ) < 8 ) ) & 287 281 & CALL ctl_warn( 'diurnal cycle for qsr: the sampling of the diurnal cycle is too small...' ) 288 289 290 ! 291 IF( ln_ssr 292 ! 293 294 ! 295 IF( nn_ice == 3 296 297 IF( nn_ice == 4 298 282 ! 283 CALL sbc_ssm_init ! Sea-surface mean fields initialisation 284 ! 285 IF( ln_ssr ) CALL sbc_ssr_init ! Sea-Surface Restoring initialisation 286 ! 287 CALL sbc_rnf_init ! Runof initialisation 288 ! 289 IF( nn_ice == 3 ) CALL sbc_lim_init ! LIM3 initialisation 290 ! 291 IF( nn_ice == 4 ) CALL cice_sbc_init( nsbc ) ! CICE initialisation 292 ! 299 293 END SUBROUTINE sbc_init 300 294 … … 327 321 vtau_b(:,:) = vtau(:,:) ! (except at nit000 where before fields 328 322 qns_b (:,:) = qns (:,:) ! are set at the end of the routine) 329 ! The 3D heat content due to qsr forcing is treated in traqsr 330 ! qsr_b (:,:) = qsr (:,:) 331 emp_b(:,:) = emp(:,:) 332 sfx_b(:,:) = sfx(:,:) 323 emp_b (:,:) = emp (:,:) 324 sfx_b (:,:) = sfx (:,:) 333 325 ENDIF 334 326 ! ! ---------------------------------------- ! … … 336 328 ! ! ---------------------------------------- ! 337 329 ! 338 IF( nn_components /= jp_iam_sas ) CALL sbc_ssm ( kt )! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)330 IF( nn_components /= jp_iam_sas ) CALL sbc_ssm ( kt ) ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m) 339 331 ! ! averaged over nf_sbc time-step 340 341 IF (ln_wave) CALL sbc_wave( kt ) 332 IF( ln_wave ) CALL sbc_wave( kt ) ! surface waves 333 334 342 335 !== sbc formulation ==! 343 336 … … 357 350 CASE( jp_mfs ) ; CALL sbc_blk_mfs ( kt ) ! bulk formulation : MFS for the ocean 358 351 CASE( jp_none ) 359 IF( nn_components == jp_iam_opa ) &360 352 IF( nn_components == jp_iam_opa ) & 353 & CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice ) ! OPA-SAS coupling: OPA receiving fields from SAS 361 354 END SELECT 362 355 363 356 IF( ln_mixcpl ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice ) ! forced-coupled mixed formulation after forcing 364 357 365 358 ! 366 359 ! !== Misc. Options ==! 367 360 ! 368 361 SELECT CASE( nn_ice ) ! Update heat and freshwater fluxes over sea-ice areas 369 362 CASE( 1 ) ; CALL sbc_ice_if ( kt ) ! Ice-cover climatology ("Ice-if" model) … … 375 368 IF( ln_icebergs ) CALL icb_stp( kt ) ! compute icebergs 376 369 377 IF( nn_isf /= 0 ) CALL sbc_isf( kt )! compute iceshelves370 IF( ln_isf ) CALL sbc_isf( kt ) ! compute iceshelves 378 371 379 372 IF( ln_rnf ) CALL sbc_rnf( kt ) ! add runoffs to fresh water fluxes … … 383 376 IF( nn_fwb /= 0 ) CALL sbc_fwb( kt, nn_fwb, nn_fsbc ) ! control the freshwater budget 384 377 385 IF( nn_closea == 1 ) CALL sbc_clo( kt ) ! treatment of closed sea in the model domain 386 ! ! (update freshwater fluxes) 378 ! treatment of closed sea in the model domain 379 ! (update freshwater fluxes) 380 ! Should not be ran if ln_diurnal_only 381 IF( .NOT.(ln_diurnal_only) .AND. (nn_closea == 1) ) CALL sbc_clo( kt ) 382 387 383 !RBbug do not understand why see ticket 667 388 384 !clem: it looks like it is necessary for the north fold (in certain circumstances). Don't know why. … … 430 426 CALL iom_rstput( kt, nitrst, numrow, 'sfx_b' , sfx ) 431 427 ENDIF 432 433 428 ! ! ---------------------------------------- ! 434 429 ! ! Outputs and control print ! … … 452 447 ! 453 448 IF(ln_ctl) THEN ! print mean trends (used for debugging) 454 CALL prt_ctl(tab2d_1=fr_i , clinfo1=' fr_i 455 CALL prt_ctl(tab2d_1=(emp-rnf + fwfisf), clinfo1=' emp-rnf 456 CALL prt_ctl(tab2d_1=(sfx-rnf + fwfisf), clinfo1=' sfx-rnf 449 CALL prt_ctl(tab2d_1=fr_i , clinfo1=' fr_i - : ', mask1=tmask, ovlap=1 ) 450 CALL prt_ctl(tab2d_1=(emp-rnf + fwfisf), clinfo1=' emp-rnf - : ', mask1=tmask, ovlap=1 ) 451 CALL prt_ctl(tab2d_1=(sfx-rnf + fwfisf), clinfo1=' sfx-rnf - : ', mask1=tmask, ovlap=1 ) 457 452 CALL prt_ctl(tab2d_1=qns , clinfo1=' qns - : ', mask1=tmask, ovlap=1 ) 458 453 CALL prt_ctl(tab2d_1=qsr , clinfo1=' qsr - : ', mask1=tmask, ovlap=1 ) -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
r5836 r6140 12 12 13 13 !!---------------------------------------------------------------------- 14 !! sbc_rnf : monthly runoffs read in a NetCDF file15 !! sbc_rnf_init : runoffs initialisation16 !! rnf_mouth : set river mouth mask14 !! sbc_rnf : monthly runoffs read in a NetCDF file 15 !! sbc_rnf_init : runoffs initialisation 16 !! rnf_mouth : set river mouth mask 17 17 !!---------------------------------------------------------------------- 18 USE dom_oce ! ocean space and time domain 19 USE phycst ! physical constants 20 USE sbc_oce ! surface boundary condition variables 21 USE sbcisf ! PM we could remove it I think 22 USE closea ! closed seas 23 USE fldread ! read input field at current time step 24 USE in_out_manager ! I/O manager 25 USE iom ! I/O module 26 USE lib_mpp ! MPP library 27 USE eosbn2 28 USE wrk_nemo ! Memory allocation 18 USE dom_oce ! ocean space and time domain 19 USE phycst ! physical constants 20 USE sbc_oce ! surface boundary condition variables 21 USE sbcisf ! PM we could remove it I think 22 USE closea ! closed seas 23 USE eosbn2 ! Equation Of State 24 ! 25 USE in_out_manager ! I/O manager 26 USE fldread ! read input field at current time step 27 USE iom ! I/O module 28 USE lib_mpp ! MPP library 29 USE wrk_nemo ! Memory allocation 29 30 30 31 IMPLICIT NONE 31 32 PRIVATE 32 33 33 PUBLIC sbc_rnf ! routinecalled in sbcmod module34 PUBLIC sbc_rnf_div ! routinecalled in divhor module35 PUBLIC sbc_rnf_alloc ! routinecalled in sbcmod module36 PUBLIC sbc_rnf_init ! routinecalled in sbcmod module34 PUBLIC sbc_rnf ! called in sbcmod module 35 PUBLIC sbc_rnf_div ! called in divhor module 36 PUBLIC sbc_rnf_alloc ! called in sbcmod module 37 PUBLIC sbc_rnf_init ! called in sbcmod module 37 38 38 ! 39 CHARACTER(len=100) :: cn_dir !: Root directory for location of rnf files39 ! !!* namsbc_rnf namelist * 40 CHARACTER(len=100) :: cn_dir !: Root directory for location of rnf files 40 41 LOGICAL :: ln_rnf_depth !: depth river runoffs attribute specified in a file 41 LOGICAL :: ln_rnf_depth_ini !: depth river runoffs computed at the initialisation42 REAL(wp) :: rn_rnf_max !: maximum value of the runoff climatologie ( ln_rnf_depth_ini = .true)43 REAL(wp) :: rn_dep_max !: depth over which runoffs is spread ( ln_rnf_depth_ini = .true)44 INTEGER :: nn_rnf_depth_file !: create (=1) a runoff depth file or not (=0)45 LOGICAL :: ln_rnf_tem !: temperature river runoffs attribute specified in a file46 LOGICAL , PUBLIC :: ln_rnf_sal !: salinity river runoffs attribute specified in a file47 TYPE(FLD_N) , PUBLIC :: sn_rnf !: information about the runoff file to be read48 TYPE(FLD_N) :: sn_cnf !: information about the runoff mouth file to be read49 TYPE(FLD_N) :: sn_s_rnf !: information about the salinities of runoff file to be read50 TYPE(FLD_N) :: sn_t_rnf !: information about the temperatures of runoff file to be read51 TYPE(FLD_N) :: sn_dep_rnf !: information about the depth which river inflow affects52 LOGICAL , PUBLIC :: ln_rnf_mouth !: specific treatment in mouths vicinity53 REAL(wp) :: rn_hrnf !: runoffs, depth over which enhanced vertical mixing is used54 REAL(wp) , PUBLIC :: rn_avt_rnf !: runoffs, value of the additional vertical mixing coef. [m2/s]55 REAL(wp) :: rn_rfact!: multiplicative factor for runoff56 57 LOGICAL , PUBLIC :: l_rnfcpl = .false. !runoffs recieved from oasis58 59 INTEGER , PUBLIC :: nkrnf = 0 !: nb of levels over which Kz is increased at river mouths42 LOGICAL :: ln_rnf_depth_ini !: depth river runoffs computed at the initialisation 43 REAL(wp) :: rn_rnf_max !: maximum value of the runoff climatologie (ln_rnf_depth_ini =T) 44 REAL(wp) :: rn_dep_max !: depth over which runoffs is spread (ln_rnf_depth_ini =T) 45 INTEGER :: nn_rnf_depth_file !: create (=1) a runoff depth file or not (=0) 46 LOGICAL :: ln_rnf_tem !: temperature river runoffs attribute specified in a file 47 LOGICAL , PUBLIC :: ln_rnf_sal !: salinity river runoffs attribute specified in a file 48 TYPE(FLD_N) , PUBLIC :: sn_rnf !: information about the runoff file to be read 49 TYPE(FLD_N) :: sn_cnf !: information about the runoff mouth file to be read 50 TYPE(FLD_N) :: sn_s_rnf !: information about the salinities of runoff file to be read 51 TYPE(FLD_N) :: sn_t_rnf !: information about the temperatures of runoff file to be read 52 TYPE(FLD_N) :: sn_dep_rnf !: information about the depth which river inflow affects 53 LOGICAL , PUBLIC :: ln_rnf_mouth !: specific treatment in mouths vicinity 54 REAL(wp) :: rn_hrnf !: runoffs, depth over which enhanced vertical mixing is used 55 REAL(wp) , PUBLIC :: rn_avt_rnf !: runoffs, value of the additional vertical mixing coef. [m2/s] 56 REAL(wp) , PUBLIC :: rn_rfact !: multiplicative factor for runoff 57 58 LOGICAL , PUBLIC :: l_rnfcpl = .false. !: runoffs recieved from oasis 59 INTEGER , PUBLIC :: nkrnf = 0 !: nb of levels over which Kz is increased at river mouths 60 60 61 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rnfmsk !: river mouth mask (hori.) 61 62 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: rnfmsk_z !: river mouth mask (vert.) … … 68 69 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_t_rnf ! structure: river runoff temperature (file information, fields read) 69 70 70 !! * Substitutions71 # include "domzgr_substitute.h90"72 71 !!---------------------------------------------------------------------- 73 72 !! NEMO/OPA 3.3 , NEMO Consortium (2010) … … 127 126 IF( ln_rnf_sal ) CALL fld_read ( kt, nn_fsbc, sf_s_rnf ) ! idem for runoffs salinity if required 128 127 ! 129 ! Runoff reduction only associated to the ORCA2_LIM configuration130 ! when reading the NetCDF file runoff_1m_nomask.nc131 IF( cp_cfg == 'orca' .AND. jp_cfg == 2 .AND. .NOT. l_rnfcpl ) THEN132 WHERE( 40._wp < gphit(:,:) .AND. gphit(:,:) < 65._wp )133 sf_rnf(1)%fnow(:,:,1) = 0.85 * sf_rnf(1)%fnow(:,:,1)134 END WHERE135 ENDIF136 !137 128 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 138 129 ! … … 142 133 IF( ln_rnf_tem ) THEN ! use runoffs temperature data 143 134 rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0 135 CALL eos_fzp( sss_m(:,:), ztfrz(:,:) ) 144 136 WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999._wp ) ! if missing data value use SST as runoffs temperature 145 137 rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0 146 138 END WHERE 147 139 WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp ) ! where fwf comes from melting of ice shelves or iceberg 148 ztfrz(:,:) = -1.9 !tfreez( sss_m(:,:) ) !PM to be discuss (trouble if sensitivity study) 149 rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * lfusisf * r1_rau0_rcp 140 rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * rlfusisf * r1_rau0_rcp 150 141 END WHERE 151 142 ELSE ! use SST as runoffs temperature … … 211 202 ! 212 203 IF( ln_rnf_depth .OR. ln_rnf_depth_ini ) THEN !== runoff distributed over several levels ==! 213 IF( lk_vvl ) THEN ! variable volume case 204 IF( ln_linssh ) THEN !* constant volume case : just apply the runoff input flow 205 DO jj = 1, jpj 206 DO ji = 1, jpi 207 DO jk = 1, nk_rnf(ji,jj) 208 phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rau0 / h_rnf(ji,jj) 209 END DO 210 END DO 211 END DO 212 ELSE !* variable volume case 214 213 DO jj = 1, jpj ! update the depth over which runoffs are distributed 215 214 DO ji = 1, jpi 216 215 h_rnf(ji,jj) = 0._wp 217 216 DO jk = 1, nk_rnf(ji,jj) ! recalculates h_rnf to be the depth in metres 218 h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk) ! to the bottom of the relevant grid box217 h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk) ! to the bottom of the relevant grid box 219 218 END DO 220 219 ! ! apply the runoff input flow … … 224 223 END DO 225 224 END DO 226 ELSE ! constant volume case : just apply the runoff input flow227 DO jj = 1, jpj228 DO ji = 1, jpi229 DO jk = 1, nk_rnf(ji,jj)230 phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rau0 / h_rnf(ji,jj)231 END DO232 END DO233 END DO234 225 ENDIF 235 226 ELSE !== runoff put only at the surface ==! 236 IF( lk_vvl ) THEN ! variable volume case 237 h_rnf(:,:) = fse3t(:,:,1) ! recalculate h_rnf to be depth of top box 238 ENDIF 239 phdivn(:,:,1) = phdivn(:,:,1) - ( rnf(:,:) + rnf_b(:,:) ) * zfact * r1_rau0 / fse3t(:,:,1) 227 h_rnf (:,:) = e3t_n (:,:,1) ! update h_rnf to be depth of top box 228 phdivn(:,:,1) = phdivn(:,:,1) - ( rnf(:,:) + rnf_b(:,:) ) * zfact * r1_rau0 / e3t_n(:,:,1) 240 229 ENDIF 241 230 ! … … 377 366 h_rnf(ji,jj) = 0._wp 378 367 DO jk = 1, nk_rnf(ji,jj) 379 h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk)368 h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk) 380 369 END DO 381 370 END DO … … 435 424 h_rnf(ji,jj) = 0._wp 436 425 DO jk = 1, nk_rnf(ji,jj) 437 h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk)426 h_rnf(ji,jj) = h_rnf(ji,jj) + e3t_n(ji,jj,jk) 438 427 END DO 439 428 END DO … … 448 437 ELSE ! runoffs applied at the surface 449 438 nk_rnf(:,:) = 1 450 h_rnf (:,:) = fse3t(:,:,1)439 h_rnf (:,:) = e3t_n(:,:,1) 451 440 ENDIF 452 441 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90
r5407 r6140 6 6 !! History : 9.0 ! 2006-07 (G. Madec) Original code 7 7 !! 3.3 ! 2010-10 (C. Bricaud, G. Madec) add the Patm forcing for sea-ice 8 !!---------------------------------------------------------------------- 9 10 !!---------------------------------------------------------------------- 11 !! sbc_ssm : calculate sea surface mean currents, temperature, 12 !! and salinity over nn_fsbc time-step 13 !!---------------------------------------------------------------------- 14 USE oce ! ocean dynamics and tracers 15 USE dom_oce ! ocean space and time domain 16 USE sbc_oce ! surface boundary condition: ocean fields 17 USE sbcapr ! surface boundary condition: atmospheric pressure 18 USE eosbn2 ! equation of state and related derivatives 8 !! 3.7 ! 2015-11 (G. Madec) non linear free surface by default: e3t_m always computed 9 !!---------------------------------------------------------------------- 10 11 !!---------------------------------------------------------------------- 12 !! sbc_ssm : calculate sea surface mean currents, temperature, 13 !! and salinity over nn_fsbc time-step 14 !!---------------------------------------------------------------------- 15 USE oce ! ocean dynamics and tracers 16 USE dom_oce ! ocean space and time domain 17 USE sbc_oce ! surface boundary condition: ocean fields 18 USE sbcapr ! surface boundary condition: atmospheric pressure 19 USE eosbn2 ! equation of state and related derivatives 19 20 ! 20 USE in_out_manager 21 USE prtctl 22 USE iom 21 USE in_out_manager ! I/O manager 22 USE prtctl ! Print control 23 USE iom ! IOM library 23 24 24 25 IMPLICIT NONE 25 26 PRIVATE 26 27 27 PUBLIC sbc_ssm ! routine called by step.F90 28 PUBLIC sbc_ssm_init ! routine called by sbcmod.F90 29 30 LOGICAL, SAVE :: l_ssm_mean = .FALSE. ! keep track of whether means have been read 31 ! from restart file 28 PUBLIC sbc_ssm ! routine called by step.F90 29 PUBLIC sbc_ssm_init ! routine called by sbcmod.F90 30 31 LOGICAL, SAVE :: l_ssm_mean = .FALSE. ! keep track of whether means have been read from restart file 32 32 33 !! * Substitutions34 # include "domzgr_substitute.h90"35 33 !!---------------------------------------------------------------------- 36 34 !! NEMO/OPA 3.3 , NEMO Consortium (2010) … … 59 57 REAL(wp), DIMENSION(jpi,jpj,jpts) :: zts 60 58 !!--------------------------------------------------------------------- 61 59 ! 62 60 ! !* surface T-, U-, V- ocean level variables (T, S, depth, velocity) 63 61 DO jj = 1, jpj … … 81 79 ENDIF 82 80 ! 83 IF( lk_vvl ) e3t_m(:,:) = fse3t_n(:,:,1)81 e3t_m(:,:) = e3t_n(:,:,1) 84 82 ! 85 83 frq_m(:,:) = fraqsr_1lev(:,:) … … 103 101 ENDIF 104 102 ! 105 IF( lk_vvl ) e3t_m(:,:) = zcoef * fse3t_n(:,:,1)103 e3t_m(:,:) = zcoef * e3t_n(:,:,1) 106 104 ! 107 105 frq_m(:,:) = zcoef * fraqsr_1lev(:,:) … … 109 107 ELSEIF( MOD( kt - 2 , nn_fsbc ) == 0 ) THEN ! Initialisation: New mean computation ! 110 108 ! ! ---------------------------------------- ! 111 ssu_m(:,:) = 0. e0! reset to zero ocean mean sbc fields112 ssv_m(:,:) = 0. e0113 sst_m(:,:) = 0. e0114 sss_m(:,:) = 0. e0115 ssh_m(:,:) = 0. e0116 IF( lk_vvl ) e3t_m(:,:) = 0.e0117 frq_m(:,:) = 0. e0109 ssu_m(:,:) = 0._wp ! reset to zero ocean mean sbc fields 110 ssv_m(:,:) = 0._wp 111 sst_m(:,:) = 0._wp 112 sss_m(:,:) = 0._wp 113 ssh_m(:,:) = 0._wp 114 e3t_m(:,:) = 0._wp 115 frq_m(:,:) = 0._wp 118 116 ENDIF 119 117 ! ! ---------------------------------------- ! … … 131 129 ENDIF 132 130 ! 133 IF( lk_vvl ) e3t_m(:,:) = fse3t_m(:,:) + fse3t_n(:,:,1)134 ! 135 frq_m(:,:) = 131 e3t_m(:,:) = e3t_m(:,:) + e3t_n(:,:,1) 132 ! 133 frq_m(:,:) = frq_m(:,:) + fraqsr_1lev(:,:) 136 134 137 135 ! ! ---------------------------------------- ! … … 139 137 ! ! ---------------------------------------- ! 140 138 zcoef = 1. / REAL( nn_fsbc, wp ) 141 sst_m(:,:) = sst_m(:,:) * zcoef 142 sss_m(:,:) = sss_m(:,:) * zcoef 143 ssu_m(:,:) = ssu_m(:,:) * zcoef 144 ssv_m(:,:) = ssv_m(:,:) * zcoef 145 ssh_m(:,:) = ssh_m(:,:) * zcoef 146 IF( lk_vvl ) e3t_m(:,:) = fse3t_m(:,:) * zcoef! mean vertical scale factor [m]147 frq_m(:,:) = frq_m(:,:) * zcoef ! mean fraction of solar net radiation absorbed in the 1st T level [-]139 sst_m(:,:) = sst_m(:,:) * zcoef ! mean SST [Celcius] 140 sss_m(:,:) = sss_m(:,:) * zcoef ! mean SSS [psu] 141 ssu_m(:,:) = ssu_m(:,:) * zcoef ! mean suface current [m/s] 142 ssv_m(:,:) = ssv_m(:,:) * zcoef ! 143 ssh_m(:,:) = ssh_m(:,:) * zcoef ! mean SSH [m] 144 e3t_m(:,:) = e3t_m(:,:) * zcoef ! mean vertical scale factor [m] 145 frq_m(:,:) = frq_m(:,:) * zcoef ! mean fraction of solar net radiation absorbed in the 1st T level [-] 148 146 ! 149 147 ENDIF … … 162 160 CALL iom_rstput( kt, nitrst, numrow, 'sss_m' , sss_m ) 163 161 CALL iom_rstput( kt, nitrst, numrow, 'ssh_m' , ssh_m ) 164 IF( lk_vvl )CALL iom_rstput( kt, nitrst, numrow, 'e3t_m' , e3t_m )162 CALL iom_rstput( kt, nitrst, numrow, 'e3t_m' , e3t_m ) 165 163 CALL iom_rstput( kt, nitrst, numrow, 'frq_m' , frq_m ) 166 164 ! … … 175 173 CALL iom_put( 'sss_m', sss_m ) 176 174 CALL iom_put( 'ssh_m', ssh_m ) 177 IF( lk_vvl )CALL iom_put( 'e3t_m', e3t_m )175 CALL iom_put( 'e3t_m', e3t_m ) 178 176 CALL iom_put( 'frq_m', frq_m ) 179 177 ENDIF 180 178 ! 181 179 END SUBROUTINE sbc_ssm 180 182 181 183 182 SUBROUTINE sbc_ssm_init … … 189 188 !! ** Action : - read parameters 190 189 !!---------------------------------------------------------------------- 191 REAL(wp) :: zcoef, zf_sbc 190 REAL(wp) :: zcoef, zf_sbc ! local scalar 192 191 !!---------------------------------------------------------------------- 193 192 ! 194 193 IF( nn_fsbc == 1 ) THEN 195 194 ! … … 206 205 IF( ln_rstart .AND. iom_varid( numror, 'nn_fsbc', ldstop = .FALSE. ) > 0 ) THEN 207 206 l_ssm_mean = .TRUE. 208 CALL iom_get( numror , 'nn_fsbc', zf_sbc ) ! sbc frequency of previous run209 CALL iom_get( numror, jpdom_autoglo, 'ssu_m' , ssu_m ) ! sea surface mean velocity (T-point)210 CALL iom_get( numror, jpdom_autoglo, 'ssv_m' , ssv_m ) ! " " velocity (V-point)211 CALL iom_get( numror, jpdom_autoglo, 'sst_m' , sst_m ) ! " " temperature (T-point)212 CALL iom_get( numror, jpdom_autoglo, 'sss_m' , sss_m ) ! " " salinity (T-point)213 CALL iom_get( numror, jpdom_autoglo, 'ssh_m' , ssh_m ) ! " " height (T-point)214 IF( lk_vvl ) CALL iom_get( numror, jpdom_autoglo, 'e3t_m', e3t_m)207 CALL iom_get( numror , 'nn_fsbc', zf_sbc ) ! sbc frequency of previous run 208 CALL iom_get( numror, jpdom_autoglo, 'ssu_m' , ssu_m ) ! sea surface mean velocity (U-point) 209 CALL iom_get( numror, jpdom_autoglo, 'ssv_m' , ssv_m ) ! " " velocity (V-point) 210 CALL iom_get( numror, jpdom_autoglo, 'sst_m' , sst_m ) ! " " temperature (T-point) 211 CALL iom_get( numror, jpdom_autoglo, 'sss_m' , sss_m ) ! " " salinity (T-point) 212 CALL iom_get( numror, jpdom_autoglo, 'ssh_m' , ssh_m ) ! " " height (T-point) 213 CALL iom_get( numror, jpdom_autoglo, 'e3t_m' , e3t_m ) ! 1st level thickness (T-point) 215 214 ! fraction of solar net radiation absorbed in 1st T level 216 215 IF( iom_varid( numror, 'frq_m', ldstop = .FALSE. ) > 0 ) THEN … … 229 228 sss_m(:,:) = zcoef * sss_m(:,:) 230 229 ssh_m(:,:) = zcoef * ssh_m(:,:) 231 IF( lk_vvl ) e3t_m(:,:) = zcoef * fse3t_m(:,:)230 e3t_m(:,:) = zcoef * e3t_m(:,:) 232 231 frq_m(:,:) = zcoef * frq_m(:,:) 233 232 ELSE … … 245 244 ELSE ; sst_m(:,:) = tsn(:,:,1,jp_tem) 246 245 ENDIF 247 sss_m(:,:) = tsn (:,:,1,jp_sal)248 ssh_m(:,:) = sshn (:,:)249 IF( lk_vvl ) e3t_m(:,:) = fse3t_n(:,:,1)246 sss_m(:,:) = tsn (:,:,1,jp_sal) 247 ssh_m(:,:) = sshn (:,:) 248 e3t_m(:,:) = e3t_n(:,:,1) 250 249 frq_m(:,:) = 1._wp 251 250 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90
r5836 r6140 19 19 ! 20 20 USE fldread ! read input fields 21 USE in_out_manager ! I/O manager 21 22 USE iom ! I/O manager 22 USE in_out_manager ! I/O manager23 23 USE lib_mpp ! distribued memory computing library 24 24 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 47 47 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_sss ! structure of input SSS (file informations, fields read) 48 48 49 !! * Substitutions50 # include "domzgr_substitute.h90"51 49 !!---------------------------------------------------------------------- 52 50 !! NEMO/OPA 4.0 , NEMO Consortium (2011) -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbctide.F90
r5930 r6140 6 6 !! History : 9.0 ! 2007 (O. Le Galloudec) Original code 7 7 !!---------------------------------------------------------------------- 8 USE oce 9 USE dom_oce 10 USE phycst 11 USE daymod 12 USE tideini 8 USE oce ! ocean dynamics and tracers variables 9 USE dom_oce ! ocean space and time domain 10 USE phycst ! physical constant 11 USE daymod ! calandar 12 USE tideini ! 13 13 ! 14 USE i om15 USE i n_out_manager ! I/O units16 USE ioipsl 17 USE lbclnk 14 USE in_out_manager ! I/O units 15 USE iom ! xIOs server 16 USE ioipsl ! NetCDF IPSL library 17 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 18 18 19 19 IMPLICIT NONE … … 46 46 INTEGER, INTENT( in ) :: kt ! ocean time-step 47 47 INTEGER :: jk ! dummy loop index 48 INTEGER :: nsec_day_orig ! Temporary variable 48 49 !!---------------------------------------------------------------------- 49 50 IF( nsec_day == NINT(0.5_wp * rdt tra(1))) THEN ! start a new day50 51 IF( nsec_day == NINT(0.5_wp * rdt) .OR. kt == nit000 ) THEN ! start a new day 51 52 ! 52 53 IF( kt == nit000 ) THEN … … 59 60 pot_astro(:,:) = 0._wp 60 61 ! 62 ! If the run does not start from midnight then need to initialise tides 63 ! at the start of the current day (only occurs when kt==nit000) 64 ! Temporarily set nsec_day to beginning of day. 65 nsec_day_orig = nsec_day 66 IF ( nsec_day /= NINT(0.5_wp * rdt) ) THEN 67 kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt 68 nsec_day = NINT(0.5_wp * rdt) 69 ELSE 70 kt_tide = kt 71 ENDIF 61 72 CALL tide_harmo( omega_tide, v0tide, utide, ftide, ntide, nb_harmo ) 62 73 ! 63 kt_tide = kt64 74 ! 65 75 IF(lwp) THEN … … 74 84 IF( ln_tide_pot ) CALL tide_init_potential 75 85 ! 86 ! Reset nsec_day 87 nsec_day = nsec_day_orig 76 88 ENDIF 77 89 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
r5860 r6140 40 40 41 41 !! * Substitutions 42 # include "domzgr_substitute.h90"43 42 # include "vectopt_loop_substitute.h90" 44 43 !!---------------------------------------------------------------------- … … 72 71 REAL(wp), DIMENSION(:,:,:), POINTER :: zusd_t, zvsd_t, ze3hdiv ! 3D workspace 73 72 !! 74 NAMELIST/namsbc_wave/ sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn 73 NAMELIST/namsbc_wave/ sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn, ln_cdgw , ln_sdw 75 74 !!--------------------------------------------------------------------- 76 75 ! … … 81 80 READ ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901) 82 81 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp ) 83 82 ! 84 83 REWIND( numnam_cfg ) ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model 85 84 READ ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 ) … … 87 86 IF(lwm) WRITE ( numond, namsbc_wave ) 88 87 ! 89 IF ( ln_cdgw ) THEN 88 IF(lwp) THEN ! Control print 89 WRITE(numout,*) ' Namelist namsbc_wave : surface wave setting' 90 WRITE(numout,*) ' wave drag coefficient ln_cdgw = ', ln_cdgw 91 WRITE(numout,*) ' wave stokes drift ln_sdw = ', ln_sdw 92 ENDIF 93 ! 94 IF( .NOT.( ln_cdgw .OR. ln_sdw ) ) & 95 & CALL ctl_warn( 'ln_sbcwave=T but nor drag coefficient (ln_cdgw=F) neither stokes drift activated (ln_sdw=F)' ) 96 IF( ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) ) & 97 & CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core') 98 ! 99 IF( ln_cdgw ) THEN 90 100 ALLOCATE( sf_cd(1), STAT=ierror ) !* allocate and fill sf_wave with sn_cdg 91 101 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) … … 97 107 cdn_wave(:,:) = 0.0 98 108 ENDIF 99 IF 109 IF( ln_sdw ) THEN 100 110 slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd; slf_i(jp_wn) = sn_wn 101 111 ALLOCATE( sf_sd(3), STAT=ierror ) !* allocate and fill sf_wave with sn_cdg … … 115 125 ENDIF 116 126 ! 117 IF ( ln_cdgw ) THEN!== Neutral drag coefficient ==!127 IF( ln_cdgw ) THEN !== Neutral drag coefficient ==! 118 128 CALL fld_read( kt, nn_fsbc, sf_cd ) ! read from external forcing 119 129 cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 120 130 ENDIF 121 131 ! 122 IF ( ln_sdw ) THEN !== Computation of the 3d Stokes Drift ==! 132 IF( ln_sdw ) THEN !== Computation of the 3d Stokes Drift ==! 133 ! 134 CALL wrk_alloc( jpi,jpj,jpk, zusd_t, zvsd_t, ze3hdiv ) 123 135 ! 124 136 CALL fld_read( kt, nn_fsbc, sf_sd ) !* read drag coefficient from external forcing 125 137 ! 126 ! 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) ) 138 DO jk = 1, jpkm1 !* distribute it on the vertical 139 zusd_t(:,:,jk) = sf_sd(jp_usd)%fnow(:,:,1) * EXP( -2._wp * sf_sd(jp_wn)%fnow(:,:,1) * gdept_n(:,:,jk) ) 140 zvsd_t(:,:,jk) = sf_sd(jp_vsd)%fnow(:,:,1) * EXP( -2._wp * sf_sd(jp_wn)%fnow(:,:,1) * gdept_n(:,:,jk) ) 132 141 END DO 133 ! !* interpolate the stokes drift from t-point to u- and v-points 134 DO jk = 1, jpkm1 142 DO jk = 1, jpkm1 !* interpolate the stokes drift from t-point to u- and v-points 135 143 DO jj = 1, jpjm1 136 144 DO ji = 1, jpim1 … … 146 154 DO jj = 2, jpjm1 147 155 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)156 ze3hdiv(ji,jj,jk) = ( e2u(ji ,jj) * e3u_n(ji ,jj,jk) * usd3d(ji ,jj,jk) & 157 & - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * usd3d(ji-1,jj,jk) & 158 & + e1v(ji,jj ) * e3v_n(ji,jj ,jk) * vsd3d(ji,jj ,jk) & 159 & - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vsd3d(ji,jj-1,jk) ) * r1_e1e2t(ji,jj) 152 160 END DO 153 161 END DO -
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90
r5930 r6140 6 6 !! History : 1.0 ! 2007 (O. Le Galloudec) Original code 7 7 !!---------------------------------------------------------------------- 8 USE oce 9 USE dom_oce 10 USE phycst 11 USE daymod 12 USE tide_mod 8 USE oce ! ocean dynamics and tracers variables 9 USE dom_oce ! ocean space and time domain 10 USE phycst ! physical constant 11 USE daymod ! calandar 12 USE tide_mod ! 13 13 ! 14 USE i om15 USE i n_out_manager ! I/O units16 USE ioipsl 17 USE lbclnk 14 USE in_out_manager ! I/O units 15 USE iom ! xIOs server 16 USE ioipsl ! NetCDF IPSL library 17 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 18 18 19 19 IMPLICIT NONE … … 27 27 LOGICAL , PUBLIC :: ln_tide_pot !: 28 28 LOGICAL , PUBLIC :: ln_tide_ramp !: 29 INTEGER , PUBLIC :: nb_harmo 30 INTEGER , PUBLIC :: kt_tide 31 REAL(wp), PUBLIC :: rdttideramp 29 INTEGER , PUBLIC :: nb_harmo !: 30 INTEGER , PUBLIC :: kt_tide !: 31 REAL(wp), PUBLIC :: rdttideramp !: 32 32 33 33 INTEGER , PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntide !: … … 40 40 CONTAINS 41 41 42 SUBROUTINE tide_init ( kt ) 43 !!---------------------------------------------------------------------- 44 !! *** ROUTINE tide_init *** 45 !!---------------------------------------------------------------------- 46 !! * Local declarations 47 INTEGER :: ji, jk 48 INTEGER, INTENT( in ) :: kt ! ocean time-step 49 CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: clname 50 INTEGER :: ios ! Local integer output status for namelist read 51 ! 52 NAMELIST/nam_tide/ln_tide_pot, ln_tide_ramp, rdttideramp, clname 53 !!---------------------------------------------------------------------- 54 55 IF ( kt == nit000 ) THEN 56 ! 57 IF(lwp) THEN 58 WRITE(numout,*) 59 WRITE(numout,*) 'tide_init : Initialization of the tidal components' 60 WRITE(numout,*) '~~~~~~~~~ ' 61 ENDIF 62 ! 63 CALL tide_init_Wave 64 ! 65 ! Read Namelist nam_tide 66 REWIND( numnam_ref ) ! Namelist nam_tide in reference namelist : Tides 67 READ ( numnam_ref, nam_tide, IOSTAT = ios, ERR = 901) 68 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_tide in reference namelist', lwp ) 69 70 REWIND( numnam_cfg ) ! Namelist nam_tide in configuration namelist : Tides 71 READ ( numnam_cfg, nam_tide, IOSTAT = ios, ERR = 902 ) 72 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_tide in configuration namelist', lwp ) 73 IF(lwm) WRITE ( numond, nam_tide ) 74 ! 75 nb_harmo=0 76 DO jk = 1, jpmax_harmo 77 DO ji = 1,jpmax_harmo 78 IF( TRIM(clname(jk)) == Wave(ji)%cname_tide ) nb_harmo = nb_harmo + 1 79 END DO 80 END DO 81 ! 82 ! Ensure that tidal components have been set in namelist_cfg 83 IF( nb_harmo .EQ. 0 ) CALL ctl_stop( 'tide_init : No tidal components set in nam_tide' ) 84 ! 85 IF(lwp) THEN 86 WRITE(numout,*) ' Namelist nam_tide' 87 WRITE(numout,*) ' Apply astronomical potential : ln_tide_pot =', ln_tide_pot 88 WRITE(numout,*) ' nb_harmo = ', nb_harmo 89 WRITE(numout,*) ' ln_tide_ramp = ', ln_tide_ramp 90 WRITE(numout,*) ' rdttideramp = ', rdttideramp 91 ENDIF 92 IF( ln_tide_ramp.AND.((nitend-nit000+1)*rdt/rday < rdttideramp) ) & 93 & CALL ctl_stop('rdttideramp must be lower than run duration') 94 IF( ln_tide_ramp.AND.(rdttideramp<0.) ) & 95 & CALL ctl_stop('rdttideramp must be positive') 96 ! 97 ! 98 ALLOCATE( ntide(nb_harmo) ) 99 DO jk = 1, nb_harmo 100 DO ji = 1, jpmax_harmo 101 IF( TRIM(clname(jk)) .eq. Wave(ji)%cname_tide ) THEN 102 ntide(jk) = ji 103 EXIT 104 END IF 105 END DO 106 END DO 107 ! 108 ALLOCATE( omega_tide(nb_harmo), v0tide (nb_harmo), & 109 & utide (nb_harmo), ftide (nb_harmo) ) 110 kt_tide = kt 111 ! 42 SUBROUTINE tide_init 43 !!---------------------------------------------------------------------- 44 !! *** ROUTINE tide_init *** 45 !!---------------------------------------------------------------------- 46 INTEGER :: ji, jk 47 CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: clname 48 INTEGER :: ios ! Local integer output status for namelist read 49 ! 50 NAMELIST/nam_tide/ln_tide_pot, ln_tide_ramp, rdttideramp, clname 51 !!---------------------------------------------------------------------- 52 ! 53 IF(lwp) THEN 54 WRITE(numout,*) 55 WRITE(numout,*) 'tide_init : Initialization of the tidal components' 56 WRITE(numout,*) '~~~~~~~~~ ' 112 57 ENDIF 58 ! 59 CALL tide_init_Wave 60 ! 61 ! Read Namelist nam_tide 62 REWIND( numnam_ref ) ! Namelist nam_tide in reference namelist : Tides 63 READ ( numnam_ref, nam_tide, IOSTAT = ios, ERR = 901) 64 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_tide in reference namelist', lwp ) 65 ! 66 REWIND( numnam_cfg ) ! Namelist nam_tide in configuration namelist : Tides 67 READ ( numnam_cfg, nam_tide, IOSTAT = ios, ERR = 902 ) 68 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_tide in configuration namelist', lwp ) 69 IF(lwm) WRITE ( numond, nam_tide ) 70 ! 71 nb_harmo=0 72 DO jk = 1, jpmax_harmo 73 DO ji = 1,jpmax_harmo 74 IF( TRIM(clname(jk)) == Wave(ji)%cname_tide ) nb_harmo = nb_harmo + 1 75 END DO 76 END DO 77 ! 78 ! Ensure that tidal components have been set in namelist_cfg 79 IF( nb_harmo == 0 ) CALL ctl_stop( 'tide_init : No tidal components set in nam_tide' ) 80 ! 81 IF(lwp) THEN 82 WRITE(numout,*) ' Namelist nam_tide' 83 WRITE(numout,*) ' Apply astronomical potential : ln_tide_pot =', ln_tide_pot 84 WRITE(numout,*) ' nb_harmo = ', nb_harmo 85 WRITE(numout,*) ' ln_tide_ramp = ', ln_tide_ramp 86 WRITE(numout,*) ' rdttideramp = ', rdttideramp 87 ENDIF 88 IF( ln_tide_ramp.AND.((nitend-nit000+1)*rdt/rday < rdttideramp) ) & 89 & CALL ctl_stop('rdttideramp must be lower than run duration') 90 IF( ln_tide_ramp.AND.(rdttideramp<0.) ) & 91 & CALL ctl_stop('rdttideramp must be positive') 92 ! 93 ALLOCATE( ntide(nb_harmo) ) 94 DO jk = 1, nb_harmo 95 DO ji = 1, jpmax_harmo 96 IF( TRIM(clname(jk)) == Wave(ji)%cname_tide ) THEN 97 ntide(jk) = ji 98 EXIT 99 ENDIF 100 END DO 101 END DO 102 ! 103 ALLOCATE( omega_tide(nb_harmo), v0tide (nb_harmo), & 104 & utide (nb_harmo), ftide (nb_harmo) ) 105 kt_tide = nit000 113 106 ! 114 107 END SUBROUTINE tide_init
Note: See TracChangeset
for help on using the changeset viewer.