C $Header: /u/gcmpack/MITgcm/pkg/profiles/profiles_init_fixed.F,v 1.39 2017/04/03 23:16:38 ou.wang Exp $ C $Name: $ #include "PROFILES_OPTIONS.h" #include "AD_CONFIG.h" C *==========================================================* C | subroutine profiles_init_fixed C | o initialization for netcdf profiles data C | started: Gael Forget 15-March-2006 C | extended: Gael Forget 14-June-2007 C *==========================================================* SUBROUTINE PROFILES_INIT_FIXED( myThid ) implicit none C ==================== Global Variables =========================== #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "DYNVARS.h" #ifdef ALLOW_CAL #include "cal.h" #endif #ifdef ALLOW_PROFILES # include "PROFILES_SIZE.h" # include "profiles.h" # include "netcdf.inc" #endif C ==================== Routine Arguments ========================== integer myThid C ==================== Routine Variables ========================== #ifdef ALLOW_PROFILES integer i,j,k,l,m,bi,bj,iG,jG,num_file,ProfNo_tile integer stopProfiles integer fid, dimid, varid1, varid1a, varid1b integer varid2,varid3 _RL tmpyymmdd(1000),tmphhmmss(1000),diffsecs _RL yymmddMin,yymmddMax _RL hhmmssMin,hhmmssMax integer tmpdate(4),tmpdiff(4),profIsInRunTime _RL tmp_lon, tmp_lon2(1000), tmp_lat2(1000), lon_cur, lat_cur _RL lon_1, lon_2, lat_1, lat_2 _RL lon_tmp1, lon_tmp2 _RL lat_fac, lon_fac integer prof_i, prof_j integer vec_start(2), vec_count(2), profno_div1000, kk character*(80) profilesfile, fnamedatanc character*(80) fnameequinc, adfnameequinc integer IL, JL, KL, err logical exst integer varid_intp1, varid_intp2, varid_intp11 , varid_intp22 integer varid_intp3, varid_intp4, varid_intp5, q, iINTERP _RL tmp_i(1000,NUM_INTERP_POINTS) _RL tmp_j(1000,NUM_INTERP_POINTS) _RL tmp_weights(1000,NUM_INTERP_POINTS),tmp_sum_weights _RL tmp_xC11(1000),tmp_yC11(1000) _RL tmp_xCNINJ(1000),tmp_yCNINJ(1000) integer stopGenericGrid Real*8 xy_buffer_r8(0:sNx+1,0:sNy+1) integer vec_start2(2), vec_count2(2) integer hh, ProfNo_hh #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST integer varid4 _RL tmp_avgbin(1000) #endif c == external functions == integer ILNBLNK EXTERNAL integer MDS_RECLEN EXTERNAL character*(max_len_mbuf) msgbuf c-- == end of interface == write(msgbuf,'(a)') ' ' call PRINT_MESSAGE( msgbuf, standardmessageunit, & SQUEEZE_RIGHT , mythid) write(msgbuf,'(a)') &'// =======================================================' call PRINT_MESSAGE( msgbuf, standardmessageunit, & SQUEEZE_RIGHT , mythid) write(msgbuf,'(a)') &'// insitu profiles model sampling >>> START <<<' call PRINT_MESSAGE( msgbuf, standardmessageunit, & SQUEEZE_RIGHT , mythid) write(msgbuf,'(a)') &'// =======================================================' call PRINT_MESSAGE( msgbuf, standardmessageunit, & SQUEEZE_RIGHT , mythid) write(msgbuf,'(a)') ' ' call PRINT_MESSAGE( msgbuf, standardmessageunit, & SQUEEZE_RIGHT , mythid) stopProfiles=0 stopGenericGrid=0 IF ( (.NOT.profilesDoGenGrid).AND. & (.NOT.usingSphericalPolarGrid .OR. rotateGrid) ) THEN WRITE(msgBuf,'(2A)') 'PROFILES_INIT_FIXED: ', & 'profilesDoGenGrid=.true. is required' CALL PRINT_ERROR( msgBuf , myThid ) WRITE(msgBuf,'(2A)') 'PROFILES_INIT_FIXED: ', & 'unless usingSphericalGrid=.TRUE. and rotateGrid=.FALSE.' CALL PRINT_ERROR( msgBuf , myThid ) CALL ALL_PROC_DIE( myThid ) STOP 'ABNORMAL END: S/R PROFILES_INIT_FIXED' ENDIF write(msgbuf,'(a)') ' ' call PRINT_MESSAGE( msgbuf, standardmessageunit, & SQUEEZE_RIGHT , mythid) write(msgbuf,'(a)') 'general packages parameters :' JL = ILNBLNK( profilesDir ) if (JL.NE.0) then write(msgbuf,'(a,a)') ' profilesDir ',profilesDir(1:JL) else write(msgbuf,'(a,a)') ' profilesDir ','./' endif call PRINT_MESSAGE( & msgbuf, standardmessageunit, SQUEEZE_RIGHT , mythid) write(msgbuf,'(a,l5)') ' profilesDoGenGrid ',profilesDoGenGrid call PRINT_MESSAGE( & msgbuf, standardmessageunit, SQUEEZE_RIGHT , mythid) write(msgbuf,'(a,l5)') ' profilesDoNcOutput ',profilesDoNcOutput call PRINT_MESSAGE( & msgbuf, standardmessageunit, SQUEEZE_RIGHT , mythid) write(msgbuf,'(a)') ' ' call PRINT_MESSAGE( msgbuf, standardmessageunit, & SQUEEZE_RIGHT , mythid) _BEGIN_MASTER( mythid ) DO bj=1,nSy DO bi=1,nSx profiles_curfile_buff(bi,bj)=0 yymmddMin=modelstartdate(1) yymmddMax=modelenddate(1) hhmmssMin=modelstartdate(2) hhmmssMax=modelenddate(2) do m=1,NLEVELMAX do l=1,1000 do k=1,NVARMAX profiles_data_buff(m,l,k,bi,bj)=0. _d 0 profiles_weight_buff(m,l,k,bi,bj)=0. _d 0 enddo enddo enddo do num_file=1,NFILESPROFMAX ProfNo_hh=0 profilesfile=' ' IL = ILNBLNK( profilesfiles(num_file) ) if (IL.NE.0) then write(profilesfile(1:80),'(1a)') & profilesfiles(num_file)(1:IL) write(msgbuf,'(a)') ' ' call PRINT_MESSAGE( msgbuf, standardmessageunit, & SQUEEZE_RIGHT , mythid) write(msgbuf,'(a,i3,a,a)') & 'profiles file #',num_file,' is ', profilesfile(1:IL) call PRINT_MESSAGE( & msgbuf, standardmessageunit, SQUEEZE_RIGHT , mythid) endif IL = ILNBLNK( profilesfile ) if (IL.NE.0) then C=========================================================== c open data files and read information C=========================================================== write(fnamedatanc(1:80),'(2a)') profilesfile(1:IL),'.nc' err = NF_OPEN(fnamedatanc, 0, fiddata(num_file,bi,bj)) c1) read the number of profiles : fid=fiddata(num_file,bi,bj) err = NF_INQ_DIMID(fid,'iPROF', dimid ) err = NF_INQ_DIMLEN(fid, dimid, ProfNo(num_file,bi,bj) ) err = NF_INQ_DIMID(fid,'iDEPTH', dimid ) if (err.NE.NF_NOERR) then err = NF_INQ_DIMID(fid,'Z', dimid ) endif err = NF_INQ_DIMLEN(fid, dimid, ProfDepthNo(num_file,bi,bj) ) err = NF_INQ_DIMID(fid,'iINTERP', dimid ) if (err.EQ.NF_NOERR) then err = NF_INQ_DIMLEN(fid, dimid, iINTERP ) else iINTERP=NUM_INTERP_POINTS endif write(msgbuf,'(a,i4,a,i4)') & ' current tile is bi,bj =', & bi,',',bj call PRINT_MESSAGE( & msgbuf, standardmessageunit, SQUEEZE_RIGHT , mythid) write(msgbuf,'(a,i9)') & ' # of depth levels in file =', & ProfDepthNo(num_file,bi,bj) call PRINT_MESSAGE( & msgbuf, standardmessageunit, SQUEEZE_RIGHT , mythid) write(msgbuf,'(a,i9)') & ' # of profiles in file =', & ProfNo(num_file,bi,bj) call PRINT_MESSAGE( & msgbuf, standardmessageunit, SQUEEZE_RIGHT , mythid) c2) read the dates and positions : err = NF_INQ_VARID(fid,'prof_depth', varid1a ) if (err.NE.NF_NOERR) then c if no prof_depth is found, then try old variable name: err = NF_INQ_VARID(fid,'depth', varid1a ) endif if (err.NE.NF_NOERR) then c if neither is found, then stop IL = ILNBLNK( profilesfile ) WRITE(msgBuf,'(3A)') & 'PROFILES_INIT_FIXED: file ', profilesfile(1:IL), & '.nc is not in the pkg/profiles format (no prof_depth etc.)' CALL PRINT_ERROR( msgBuf, myThid) stopProfiles=1 endif do k=1,ProfDepthNo(num_file,bi,bj) err = NF_GET_VAR1_DOUBLE(fid,varid1a,k, & prof_depth(num_file,k,bi,bj)) enddo err = NF_INQ_VARID(fid,'prof_YYYYMMDD', varid1a ) err = NF_INQ_VARID(fid,'prof_HHMMSS', varid1b ) err = NF_INQ_VARID(fid,'prof_lon', varid2 ) err = NF_INQ_VARID(fid,'prof_lat', varid3 ) #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST err = NF_INQ_VARID(fid,'prof_bin_id_a', varid4 ) #endif if (err.NE.NF_NOERR) then IL = ILNBLNK( profilesfile ) WRITE(msgBuf,'(3A)') & 'PROFILES_INIT_FIXED: file ', profilesfile(1:IL), & '.nc is not in the pkg/profiles format (no prof_YYYYMMDD etc.)' CALL PRINT_ERROR( msgBuf, myThid) stopProfiles=1 endif if (profilesDoGenGrid) then c3) read interpolattion information (grid points, coeffs, etc.) err = NF_INQ_VARID(fid,'prof_interp_XC11',varid_intp1) err = NF_INQ_VARID(fid,'prof_interp_YC11',varid_intp2) err = NF_INQ_VARID(fid,'prof_interp_XCNINJ',varid_intp11) err = NF_INQ_VARID(fid,'prof_interp_YCNINJ',varid_intp22) err = NF_INQ_VARID(fid,'prof_interp_weights',varid_intp3) err = NF_INQ_VARID(fid,'prof_interp_i',varid_intp4) err = NF_INQ_VARID(fid,'prof_interp_j',varid_intp5) if (err.NE.NF_NOERR) then IL = ILNBLNK( profilesfile ) WRITE(msgBuf,'(3A)') & 'PROFILES_INIT_FIXED: file ', profilesfile(1:IL), & '.nc is missing interpolation information (profilesDoGenGrid)' CALL PRINT_ERROR( msgBuf, myThid) stopGenericGrid=2 endif endif c4) default values do k=1,NOBSGLOB prof_time(num_file,k,bi,bj)=-999. _d 0 prof_lon(num_file,k,bi,bj)=-999. _d 0 prof_lat(num_file,k,bi,bj)=-999. _d 0 prof_ind_glob(num_file,k,bi,bj)=0 #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST prof_ind_avgbin(num_file,k,bi,bj)=-999 #endif do q = 1,NUM_INTERP_POINTS prof_interp_i(num_file,k,q,bi,bj) = 1 prof_interp_j(num_file,k,q,bi,bj) = 1 prof_interp_weights(num_file,k,q,bi,bj) = 0. _d 0 enddo prof_interp_xC11(num_file,k,bi,bj)=-999. _d 0 prof_interp_yC11(num_file,k,bi,bj)=-999. _d 0 prof_interp_xCNINJ(num_file,k,bi,bj)=-999. _d 0 prof_interp_yCNINJ(num_file,k,bi,bj)=-999. _d 0 enddo c5) main loop: look for profiles in this tile ProfNo_tile=0 profno_div1000=max(0,int(ProfNo(num_file,bi,bj)/1000)) do kk=1,profno_div1000+1 if (min(ProfNo(num_file,bi,bj), 1000*kk).GE. & 1+1000*(kk-1)) then c5.1) read a chunk vec_start(1)=1 vec_start(2)=1+1000*(kk-1) vec_count(1)=1 vec_count(2)=min(1000,ProfNo(num_file,bi,bj)-1000*(kk-1)) if ( (vec_count(2).LE.0).OR.(vec_count(2).GT.1000).OR. & (vec_start(2).LE.0).OR. & (vec_count(2)+vec_start(2)-1.GT.ProfNo(num_file,bi,bj)) ) & then IL = ILNBLNK( profilesfile ) WRITE(msgBuf,'(3A)') & 'PROFILES_INIT_FIXED: file ', profilesfile(1:IL), & '.nc was not read properly (case 1).' CALL PRINT_ERROR( msgBuf, myThid) stopProfiles=1 endif err = NF_GET_VARA_DOUBLE(fid,varid1a,vec_start(2), & vec_count(2), tmpyymmdd) err = NF_GET_VARA_DOUBLE(fid,varid1b,vec_start(2), & vec_count(2), tmphhmmss) err = NF_GET_VARA_DOUBLE(fid,varid2,vec_start(2), & vec_count(2), tmp_lon2) err = NF_GET_VARA_DOUBLE(fid,varid3,vec_start(2), & vec_count(2), tmp_lat2) #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST err = NF_GET_VARA_DOUBLE(fid,varid4,vec_start(2), & vec_count(2), tmp_avgbin) #endif if (err.NE.NF_NOERR) then WRITE(msgBuf,'(3A)') & 'PROFILES_INIT_FIXED: file ', profilesfile(1:IL), & '.nc was not read properly (case 2).' CALL PRINT_ERROR( msgBuf, myThid) stopProfiles=1 endif c if profilesDoGenGrid then also read in the interpolation coeffs and indices if (profilesDoGenGrid) then err = NF_GET_VARA_DOUBLE(fid,varid_intp1,vec_start(2), & vec_count(2), tmp_xC11) err = NF_GET_VARA_DOUBLE(fid,varid_intp2,vec_start(2), & vec_count(2), tmp_yC11) err = NF_GET_VARA_DOUBLE(fid,varid_intp11,vec_start(2), & vec_count(2), tmp_xCNINJ) err = NF_GET_VARA_DOUBLE(fid,varid_intp22,vec_start(2), & vec_count(2), tmp_yCNINJ) do q=1,iINTERP vec_start2(1)=q vec_start2(2)=1+1000*(kk-1) vec_count2(1)=1 vec_count2(2)=min(1000,ProfNo(num_file,bi,bj)-1000*(kk-1)) err = NF_GET_VARA_DOUBLE(fid,varid_intp3,vec_start2, & vec_count2, tmp_weights(1,q)) err = NF_GET_VARA_DOUBLE(fid,varid_intp4,vec_start2, & vec_count2, tmp_i(1,q)) err = NF_GET_VARA_DOUBLE(fid,varid_intp5,vec_start2, & vec_count2, tmp_j(1,q)) enddo endif c5.2) loop through this chunk do k=1,min(1000,ProfNo(num_file,bi,bj)-1000*(kk-1)) if ( stopProfiles .EQ. 0) then profIsInRunTime=1 IF (( ( tmpyymmdd(k).GT.yymmddMin ).OR.(( tmpyymmdd(k).EQ. & yymmddMin ).AND.( tmphhmmss(k).GT.hhmmssMin ))).AND. & ( ( tmpyymmdd(k).LT.yymmddMax ).OR.(( tmpyymmdd(k).EQ. & yymmddMax ).AND.( tmphhmmss(k).LT.hhmmssMax ))) ) THEN hh = int(tmphhmmss(k))/10000 IF ( hh.LT.hoursPerDay ) THEN profIsInRunTime=1 call CAL_FULLDATE( int(tmpyymmdd(k)),int(tmphhmmss(k)), & tmpdate,mythid ) call CAL_TIMEPASSED( modelstartdate,tmpdate,tmpdiff,mythid ) call CAL_TOSECONDS (tmpdiff,diffsecs,mythid) diffsecs=diffsecs+nIter0*deltaTclock ELSE c if tmphhmmss is out of range then disregard profile profIsInRunTime=0 diffsecs=-deltaTclock ProfNo_hh=ProfNo_hh+1 ENDIF ELSE profIsInRunTime=0 diffsecs=-deltaTclock ENDIF c ============================================================================== c 5.2a) determine whether profiles is in current tile domain (lat-lon grid case) if ((.NOT.profilesDoGenGrid).AND.(profIsInRunTime.EQ.1)) then if (xC(sNx+1,1,bi,bj).LT.xC(1,1,bi,bj)) then tmp_lon=xC(sNx+1,1,bi,bj)+360. _d 0 else tmp_lon=xC(sNx+1,1,bi,bj) endif if ((xC(1,1,bi,bj).LE.tmp_lon2(k)).AND. & (tmp_lon.GT.tmp_lon2(k)).AND. & (yC(1,1,bi,bj).LE.tmp_lat2(k)).AND. & (yC(1,sNy+1,bi,bj).GT.tmp_lat2(k)) ) then lon_cur=tmp_lon2(k) lat_cur=tmp_lat2(k) elseif ((xC(sNx+1,1,bi,bj).LT.xC(1,1,bi,bj)).AND. & (xC(1,1,bi,bj).LE.tmp_lon2(k)+360. _d 0).AND. & (tmp_lon.GT.tmp_lon2(k)+360. _d 0).AND. & (yC(1,1,bi,bj).LE.tmp_lat2(k)).AND. & (yC(1,sNy+1,bi,bj).GT.tmp_lat2(k)) & ) then lon_cur=tmp_lon2(k)+360. _d 0 lat_cur=tmp_lat2(k) else profIsInRunTime=0 endif c now determine value of i,j to the south-ouest of data point prof_i=-10 prof_j=-10 lon_1=-10 lon_2=-10 lat_1=-10 lat_2=-10 if (profIsInRunTime.EQ.1) then DO j=1,sNy+1 DO i=1,sNx+1 c value of j, south of the data point: if ((yC(i,j,bi,bj).LE.lat_cur).AND. & (yC(i,j+1,bi,bj).GT.lat_cur)) then prof_j=j lat_1=yC(i,j,bi,bj) lat_2=yC(i,j+1,bi,bj) endif c value of i, west of the data point: if (xC(i+1,j,bi,bj).LT.xC(1,j,bi,bj)) then lon_tmp2=xC(i+1,j,bi,bj)+360 else lon_tmp2=xC(i+1,j,bi,bj) endif if (xC(i,j,bi,bj).LT.xC(1,j,bi,bj)) then lon_tmp1=xC(i,j,bi,bj)+360 else lon_tmp1=xC(i,j,bi,bj) endif if ((lon_tmp1.LE.lon_cur).AND.(lon_tmp2.GT.lon_cur)) then prof_i=i lon_1=lon_tmp1 lon_2=lon_tmp2 endif ENDDO ENDDO endif if ((prof_i.EQ.-10).OR.(prof_j.EQ.-10)) profIsInRunTime=0 if (profIsInRunTime.EQ.1) then c if yes then store prof_time and longitude and latitude: ProfNo_tile=ProfNo_tile+1 prof_time(num_file,ProfNo_tile,bi,bj)=diffsecs prof_lon(num_file,ProfNo_tile,bi,bj)=lon_cur prof_lat(num_file,ProfNo_tile,bi,bj)=lat_cur prof_ind_glob(num_file,ProfNo_tile,bi,bj)=k+1000*(kk-1) #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST prof_ind_avgbin(num_file,ProfNo_tile,bi,bj)=tmp_avgbin(k) #endif c then store interpolation coeffs and indices lon_fac=(lon_cur-lon_1)/(lon_2-lon_1) lat_fac=(lat_cur-lat_1)/(lat_2-lat_1) prof_interp_weights(num_file,ProfNo_tile,1,bi,bj)= & (1-lon_fac)*(1-lat_fac) prof_interp_i(num_file,ProfNo_tile,1,bi,bj)=prof_i prof_interp_j(num_file,ProfNo_tile,1,bi,bj)=prof_j prof_interp_weights(num_file,ProfNo_tile,2,bi,bj)= & lon_fac*(1-lat_fac) prof_interp_i(num_file,ProfNo_tile,2,bi,bj)=prof_i+1 prof_interp_j(num_file,ProfNo_tile,2,bi,bj)=prof_j prof_interp_weights(num_file,ProfNo_tile,3,bi,bj)= & (1-lon_fac)*lat_fac prof_interp_i(num_file,ProfNo_tile,3,bi,bj)=prof_i prof_interp_j(num_file,ProfNo_tile,3,bi,bj)=prof_j+1 prof_interp_weights(num_file,ProfNo_tile,4,bi,bj)= & lon_fac*lat_fac prof_interp_i(num_file,ProfNo_tile,4,bi,bj)=prof_i+1 prof_interp_j(num_file,ProfNo_tile,4,bi,bj)=prof_j+1 endif c ============================================================================== c 5.2a) determine whether profiles is in current tile domain (generic grid case) elseif (profIsInRunTime.EQ.1) then if (stopGenericGrid.EQ.0) then if ( ( abs( tmp_xC11(k) - xC(1,1,bi,bj) ).LT.0.0001 _d 0 ) .AND. & ( abs( tmp_yC11(k) - yC(1,1,bi,bj) ).LT.0.0001 _d 0) .AND. & ( abs( tmp_xCNINJ(k) - xC(sNx,sNy,bi,bj) ).LT.0.0001 _d 0 ) .AND. & ( abs( tmp_yCNINJ(k) - yC(sNx,sNy,bi,bj) ).LT.0.0001 _d 0 ) & .AND.(profIsInRunTime.EQ.1)) then c if yes then store prof_time and interpolation coeffs and indices: ProfNo_tile=ProfNo_tile+1 prof_time(num_file,ProfNo_tile,bi,bj)=diffsecs #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST prof_ind_avgbin(num_file,ProfNo_tile,bi,bj)=tmp_avgbin(k) #endif prof_interp_xC11(num_file,ProfNo_tile,bi,bj)=tmp_xC11(k) prof_interp_yC11(num_file,ProfNo_tile,bi,bj)=tmp_yC11(k) prof_interp_xCNINJ(num_file,ProfNo_tile,bi,bj)=tmp_xCNINJ(k) prof_interp_yCNINJ(num_file,ProfNo_tile,bi,bj)=tmp_yCNINJ(k) tmp_sum_weights=0. _d 0 do q = 1,iINTERP prof_interp_weights(num_file,ProfNo_tile,q,bi,bj) & =tmp_weights(k,q) prof_interp_i(num_file,ProfNo_tile,q,bi,bj) & =tmp_i(k,q) prof_interp_j(num_file,ProfNo_tile,q,bi,bj) & =tmp_j(k,q) tmp_sum_weights=tmp_sum_weights+tmp_weights(k,q) c more test of the inputs: is the offline-computed c interpolation information consistent (self and with grid) if ( (tmp_i(k,q).LT.0).OR.(tmp_j(k,q).LT.0) & .OR.(tmp_i(k,q).GT.sNx+1).OR.(tmp_j(k,q).GT.sNy+1) ) then WRITE(msgBuf,'(4A)') & 'PROFILES_INIT_FIXED: file ', profilesfile(1:IL), & '.nc includes inconsistent interpolation ', & 'points (profilesDoGenGrid; out of tile)' CALL PRINT_ERROR( msgBuf, myThid) stopGenericGrid=1 endif #ifdef ALLOW_PROFILES_EXCLUDE_CORNERS if ( tmp_weights(k,q) .NE. 0. _d 0) then if ( ((tmp_i(k,q).EQ.0).AND.(tmp_j(k,q).EQ.0)) & .OR.((tmp_i(k,q).EQ.sNx+1).AND.(tmp_j(k,q).EQ.sNy+1)) & .OR.((tmp_i(k,q).EQ.0).AND.(tmp_j(k,q).EQ.sNy+1)) & .OR.((tmp_i(k,q).EQ.sNx+1).AND.(tmp_j(k,q).EQ.0)) ) then WRITE(msgBuf,'(4A)') & 'PROFILES_INIT_FIXED: file ', profilesfile(1:IL), & '.nc includes inconsistent interpolation ', & 'points (profilesDoGenGrid; using overlap corners)' CALL PRINT_ERROR( msgBuf, myThid) stopGenericGrid=1 endif endif #endif /* ALLOW_PROFILES_EXCLUDE_CORNERS */ if ( (tmp_weights(k,q).LT.0. _d 0).OR. & (tmp_weights(k,q).GT.1. _d 0) ) then WRITE(msgBuf,'(4A)') & 'PROFILES_INIT_FIXED: file ', profilesfile(1:IL), & '.nc includes inconsistent interpolation ', & 'weights (profilesDoGenGrid; sum oustide 0-1)' CALL PRINT_ERROR( msgBuf, myThid) stopGenericGrid=1 endif enddo if ( abs(tmp_sum_weights -1. _d 0 ) .GT. 0.0001 _d 0) then WRITE(msgBuf,'(4A)') & 'PROFILES_INIT_FIXED: file ', profilesfile(1:IL), & '.nc includes inconsistent interpolation ', & 'weights (profilesDoGenGrid; dont add up to 1)' CALL PRINT_ERROR( msgBuf, myThid) stopGenericGrid=1 endif prof_ind_glob(num_file,ProfNo_tile,bi,bj)=k+1000*(kk-1) endif endif endif !if (.NOT.profilesDoGenGrid) then c ============================================================================== c check that maximum size was not reached: if (ProfNo_tile.GE.NOBSGLOB) then WRITE(msgBuf,'(3A)') & 'PROFILES_INIT_FIXED: file ', profilesfile(1:IL), & '.nc was not read properly (increase NOBSGLOB).' CALL PRINT_ERROR( msgBuf, myThid) stopProfiles=1 endif endif !if ( stopProfiles .EQ. 0) then enddo !do k=1,min(1000,ProfNo(num_file,bi,bj)-1000*(kk-1)) endif !if (min(ProfNo(num_file,bi,bj), 1000... enddo !do kk=1,profno_div1000+1 ProfNo(num_file,bi,bj)=ProfNo_tile write(msgbuf,'(a,i9)') & ' # of profiles with erroneous HHMMSS values =', & ProfNo_hh call PRINT_MESSAGE( & msgbuf, standardmessageunit, SQUEEZE_RIGHT , mythid) write(msgbuf,'(a,i9)') & ' # of profiles within tile and time period =', & ProfNo(num_file,bi,bj) call PRINT_MESSAGE( & msgbuf, standardmessageunit, SQUEEZE_RIGHT , mythid) c6) available variablesin the data set do k=1,NVARMAX prof_num_var_cur(num_file,k,bi,bj)=0 enddo prof_num_var_tot(num_file,bi,bj)=0 do k=1,NVARMAX JL = ILNBLNK( prof_names(num_file,k) ) err = NF_INQ_VARID(fid,prof_names(num_file,k)(1:JL), varid1 ) if (err.EQ.NF_NOERR) then vec_quantities(num_file,k,bi,bj)=.TRUE. prof_num_var_tot(num_file,bi,bj)= & prof_num_var_tot(num_file,bi,bj)+1 prof_num_var_cur(num_file,k,bi,bj)= & prof_num_var_tot(num_file,bi,bj) else vec_quantities(num_file,k,bi,bj)=.FALSE. endif enddo do k=1,NVARMAX if (vec_quantities(num_file,k,bi,bj)) then KL = ILNBLNK( prof_names(num_file,k) ) JL = ILNBLNK( prof_namesmod(num_file,k) ) if (prof_namesmod(num_file,k).EQ.'pTracer') then write(msgbuf,'(a,I3,5a,I3)') ' variable #',k,' is ' , & prof_names(num_file,k)(1:KL),' and ', & prof_namesmod(num_file,k)(1:JL),' #', & prof_itracer(num_file,k) else write(msgbuf,'(a,I3,4a)') ' variable #',k, & ' is ' , & prof_names(num_file,k)(1:KL),' and ', & prof_namesmod(num_file,k)(1:JL) endif call PRINT_MESSAGE(msgbuf, & standardmessageunit, SQUEEZE_RIGHT , mythid) endif enddo C=========================================================== c create files for model counterparts to observations C=========================================================== if (ProfNo(num_file,bi,bj).GT.0) then iG=bi+(myXGlobalLo-1)/sNx ! Kludge until unstructered tiles jG=bj+(myYGlobalLo-1)/sNy ! Kludge until unstructered tiles JL = ILNBLNK( profilesDir ) if (profilesDoNcOutput) then write(fnameequinc(1:80),'(3a,i3.3,a,i3.3,a)') & profilesDir(1:JL),profilesfile(1:IL),'.',iG,'.',jG,'.equi.nc' write(adfnameequinc(1:80),'(4a,i3.3,a,i3.3,a)') & profilesDir(1:JL),'ad', & profilesfile(1:IL),'.',iG,'.',jG,'.equi.nc' inquire( file=fnameequinc, exist=exst ) if (.NOT.exst) then call PROFILES_INIT_NCFILE(num_file, & fiddata(num_file,bi,bj),fnameequinc, & fidforward(num_file,bi,bj),ProfNo(num_file,bi,bj), & ProfDepthNo(num_file,bi,bj), & bi,bj,myThid) else err = NF_OPEN(fnameequinc,NF_WRITE,fidforward(num_file,bi,bj)) endif #ifdef ALLOW_ADJOINT_RUN inquire( file=adfnameequinc, exist=exst ) if (.NOT.exst) then call PROFILES_INIT_NCFILE(num_file,fiddata(num_file,bi,bj), & adfnameequinc, fidadjoint(num_file,bi,bj), & ProfNo(num_file,bi,bj),ProfDepthNo(num_file,bi,bj), & bi,bj, myThid) else err = NF_OPEN(adfnameequinc,NF_WRITE,fidadjoint(num_file,bi,bj)) endif #endif else write(fnameequinc(1:80),'(3a,i3.3,a,i3.3,a)') & profilesDir(1:JL),profilesfile(1:IL),'.',iG,'.',jG,'.equi.data' write(adfnameequinc(1:80),'(4a,i3.3,a,i3.3,a)') & profilesDir(1:JL),'ad', & profilesfile(1:IL),'.',iG,'.',jG,'.equi.data' inquire( file=fnameequinc, exist=exst ) #ifdef PROFILES_USE_MDSFINDUNITS call MDSFINDUNIT( fidforward(num_file,bi,bj) , mythid ) #else call PROFILES_FINDUNIT( fidforward(num_file,bi,bj) , mythid ) #endif if (.NOT.exst) then call PROFILES_INIT_NCFILE(num_file,fiddata(num_file,bi,bj), & fnameequinc,fidforward(num_file,bi,bj), & ProfNo(num_file,bi,bj),ProfDepthNo(num_file,bi,bj), & bi,bj,myThid) else open( fidforward(num_file,bi,bj),file=fnameequinc, & form ='unformatted',status='unknown', access='direct', & recl= (ProfDepthNo(num_file,bi,bj)+1)*WORDLENGTH*2 ) endif #ifdef ALLOW_ADJOINT_RUN inquire( file=adfnameequinc, exist=exst ) #ifdef PROFILES_USE_MDSFINDUNITS call MDSFINDUNIT( fidadjoint(num_file,bi,bj) , mythid ) #else call PROFILES_FINDUNIT( fidadjoint(num_file,bi,bj) , mythid ) #endif if (.NOT.exst) then call PROFILES_INIT_NCFILE(num_file,fiddata(num_file,bi,bj), & adfnameequinc, fidadjoint(num_file,bi,bj), & ProfNo(num_file,bi,bj),ProfDepthNo(num_file,bi,bj), & bi,bj, myThid) else open( fidadjoint(num_file,bi,bj),file=adfnameequinc, & form ='unformatted',status='unknown', access='direct', & recl= (ProfDepthNo(num_file,bi,bj)+1)*WORDLENGTH*2 ) endif #endif endif endif C=========================================================== else ProfNo(num_file,bi,bj)=0 do k=1,NVARMAX prof_num_var_cur(num_file,k,bi,bj)=0 vec_quantities(num_file,k,bi,bj)=.FALSE. enddo prof_num_var_tot(num_file,bi,bj)=0 do k=1,NOBSGLOB prof_time(num_file,k,bi,bj)=-999. _d 0 prof_lon(num_file,k,bi,bj)=-999. _d 0 prof_lat(num_file,k,bi,bj)=-999. _d 0 prof_ind_glob(num_file,k,bi,bj)=0 #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST prof_ind_avgbin(num_file,k,bi,bj)=-999 #endif do q = 1,NUM_INTERP_POINTS prof_interp_i(num_file,k,q,bi,bj) = 1 prof_interp_j(num_file,k,q,bi,bj) = 1 prof_interp_weights(num_file,k,q,bi,bj) = 0. _d 0 enddo prof_interp_xC11(num_file,k,bi,bj)=-999. _d 0 prof_interp_yC11(num_file,k,bi,bj)=-999. _d 0 prof_interp_xCNINJ(num_file,k,bi,bj)=-999. _d 0 prof_interp_yCNINJ(num_file,k,bi,bj)=-999. _d 0 enddo endif !if (IL.NE.0) then enddo ! do num_file=1,NFILESPROFMAX #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST C Find the unique depth levels from all profile datasets C initialize prof_depth_comb if(bi.EQ.1.and.bj.EQ.1)then NLEVELCOMB = 0 NLEVELCOMBRL = NLEVELCOMB endif do m=1,NLEVELCOMBMAX prof_depth_comb(m,bi,bj)=-999. _d 0 enddo m = 1 do num_file=1,NFILESPROFMAX do K=1,ProfDepthNo(num_file,bi,bj) if(m.EQ.1) then prof_depth_comb(m,bi,bj) = prof_depth(num_file, k,bi,bj) m = m + 1 else C sort do l=1,NLEVELCOMBMAX-1 if(prof_depth_comb(l,bi,bj) .ne. -999. _d 0) then if(prof_depth(num_file, k,bi,bj).lt. & prof_depth_comb(l,bi,bj).and. & l.EQ.1) then prof_depth_comb(NLEVELCOMBMAX,bi,bj) = & prof_depth_comb(l,bi,bj) prof_depth_comb(l,bi,bj)= & prof_depth(num_file, k,bi,bj) do il = NLEVELCOMBMAX-1, l+2,-1 prof_depth_comb(il,bi,bj)= & prof_depth_comb(il-1,bi,bj) enddo prof_depth_comb(l+1,bi,bj)= & prof_depth_comb(NLEVELCOMBMAX,bi,bj) else if(prof_depth(num_file, k,bi,bj).gt. & prof_depth_comb(l,bi,bj).and. & prof_depth(num_file, k,bi,bj).lt. & prof_depth_comb(l+1,bi,bj)) then prof_depth_comb(NLEVELCOMBMAX,bi,bj) = & prof_depth_comb(l+1,bi,bj) prof_depth_comb(l+1,bi,bj)= & prof_depth(num_file, k,bi,bj) do il = NLEVELCOMBMAX-1, l+3,-1 prof_depth_comb(il,bi,bj)= & prof_depth_comb(il-1,bi,bj) enddo prof_depth_comb(l+2,bi,bj)= & prof_depth_comb(NLEVELCOMBMAX,bi,bj) else if ( prof_depth(num_file, k,bi,bj).gt. & prof_depth_comb(l,bi,bj).and. & prof_depth_comb(l+1,bi,bj).eq.-999. _d 0) then prof_depth_comb(l+1,bi,bj) = & prof_depth(num_file, k,bi,bj) endif endif enddo endif if(m.GE.NLEVELCOMBMAX-2)then WRITE(msgBuf,'(A)') & 'increase NLEVELCOMBMAX' CALL PRINT_ERROR( msgBuf, myThid) endif enddo ! do K=1,ProfDepthNo(num_file,bi,bj) enddo ! do num_file=1,NFILESPROFMAX prof_depth_comb(NLEVELCOMBMAX,bi,bj) = -999. _d 0 C diagnostics output do m=1,NLEVELCOMBMAX if(prof_depth_comb(m,bi,bj) .GE. 0. _d 0 & .AND. NLEVELCOMB.LT.m)then NLEVELCOMB = m if(m.GE.NLEVELCOMBMAX-2)then WRITE(msgBuf,'(A,2i6)') & 'increase NLEVELCOMBMAX: m,NLEVELCOMBMA ', & m, NLEVELCOMBMAX CALL PRINT_ERROR( msgBuf, myThid) endif endif enddo WRITE(msgBuf,'(A, i6,d20.5)') & 'NLEVELCOMB = ', NLEVELCOMB call PRINT_MESSAGE( msgbuf, standardmessageunit, & SQUEEZE_RIGHT , mythid) #endif C=========================================================== C error cases: C=========================================================== c1) you want to provide interpolation information if ( stopGenericGrid.EQ.2) then iG=bi+(myXGlobalLo-1)/sNx ! Kludge until unstructered tiles jG=bj+(myYGlobalLo-1)/sNy ! Kludge until unstructered tiles cgf XC grid call MDSFINDUNIT( fid , mythid ) write(fnameequinc(1:80),'(a,i3.3,a,i3.3,a,i4.4,a,i4.4,a)') & 'profilesXCincl1PointOverlap.',iG,'.',jG,'.',sNx,'.',sNy,'.data' k=MDS_RECLEN(64,(sNx+2)*(sNy+2),mythid) WRITE(standardMessageUnit,'(A,/,2A)') & 'PROFILES_INIT_FIXED: creating grid from profiles; file:', & fnameequinc open( fid, file= fnameequinc, form ='unformatted', & status='unknown',access='direct', recl= k) DO m=0,sNy+1 DO l=0,sNx+1 xy_buffer_r8(l,m)=xC(l,m,bi,bj) ENDDO ENDDO #ifdef _BYTESWAPIO call MDS_BYTESWAPR8((sNx+2)*(sNy+2),xy_buffer_r8) #endif write(fid,rec=1) xy_buffer_r8 close(fid) cgf YC grid call MDSFINDUNIT( fid , mythid ) write(fnameequinc(1:80),'(a,i3.3,a,i3.3,a,i4.4,a,i4.4,a)') & 'profilesYCincl1PointOverlap.',iG,'.',jG,'.',sNx,'.',sNy,'.data' k=MDS_RECLEN(64,(sNx+2)*(sNy+2),mythid) WRITE(standardMessageUnit,'(A,/,A)') & 'PROFILES_INIT_FIXED: creating grid from profiles; file:', & fnameequinc open( fid, file= fnameequinc, form ='unformatted', & status='unknown', access='direct', recl= k) DO m=0,sNy+1 DO l=0,sNx+1 xy_buffer_r8(l,m)=yC(l,m,bi,bj) ENDDO ENDDO #ifdef _BYTESWAPIO call MDS_BYTESWAPR8((sNx+2)*(sNy+2),xy_buffer_r8) #endif write(fid,rec=1) xy_buffer_r8 close(fid) WRITE(msgBuf,'(3A)') & 'PROFILES_INIT_FIXED : ', & 'when using profilesDoGenGrid ', & 'you have to provide interpolation coeffs etc. ' CALL PRINT_ERROR( msgBuf, myThid) WRITE(msgBuf,'(2A)') & 'and some of your nc files dont have them. ', & 'You could use profiles_prep_mygrid.m and/or' CALL PRINT_ERROR( msgBuf, myThid) WRITE(msgBuf,'(A)') & 'use the grid info in profiles*incl1PointOverlap*data' CALL PRINT_ERROR( msgBuf, myThid) stopProfiles=1 endif ENDDO ENDDO #ifdef ALLOW_PROFILES_SAMPLESPLIT_COST NLEVELCOMBRL = NLEVELCOMB _GLOBAL_MAX_RL( NLEVELCOMBRL, myThid ) NLEVELCOMB = NLEVELCOMBRL #endif _END_MASTER( mythid ) _BARRIER c2) stop after other kind of errors CALL GLOBAL_SUM_INT( stopProfiles , myThid ) if ( stopProfiles.GE.1) then CALL ALL_PROC_DIE( myThid ) STOP 'ABNORMAL END: S/R PROFILES_INIT_FIXED' endif CALL GLOBAL_SUM_INT( stopGenericGrid , myThid ) if ( stopGenericGrid.GE.1) then CALL ALL_PROC_DIE( myThid ) STOP 'ABNORMAL END: S/R PROFILES_INIT_FIXED' endif write(msgbuf,'(a)') ' ' call PRINT_MESSAGE( msgbuf, standardmessageunit, & SQUEEZE_RIGHT , mythid) write(msgbuf,'(a)') &'// =======================================================' call PRINT_MESSAGE( msgbuf, standardmessageunit, & SQUEEZE_RIGHT , mythid) write(msgbuf,'(a)') &'// insitu profiles model sampling >>> END <<<' call PRINT_MESSAGE( msgbuf, standardmessageunit, & SQUEEZE_RIGHT , mythid) write(msgbuf,'(a)') &'// =======================================================' call PRINT_MESSAGE( msgbuf, standardmessageunit, & SQUEEZE_RIGHT , mythid) write(msgbuf,'(a)') ' ' call PRINT_MESSAGE( msgbuf, standardmessageunit, & SQUEEZE_RIGHT , mythid) #endif RETURN END