Index: trunk/LMDZ.MARS/libf/phymars/surfini.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/surfini.F (revision 2501) +++ trunk/LMDZ.MARS/libf/phymars/surfini.F (revision 2502) @@ -4,5 +4,6 @@ use netcdf use tracer_mod, only: nqmx, noms - use geometry_mod, only: longitude, latitude ! in radians + use geometry_mod, only: longitude, latitude, ! in radians + & cell_area ! for watercaptag diagnosis use surfdat_h, only: watercaptag, frost_albedo_threshold, & albedo_h2o_ice, inert_h2o_ice, albedodat, @@ -36,10 +37,12 @@ LOGICAL, DIMENSION(100000) :: longwatercaptag -! There are 3 different modes for ice distribution: +! There are 4 different modes for ice distribution: ! icelocationmode = 1 ---> based on data from surface.nc ! icelocationmode = 2 ---> directly predefined for GCM resolutions 32x24 or 64x48 ! icelocationmode = 3 ---> based on logical relations for latitude and longitude +! icelocationmode = 4 ---> predefined 64x48 but usable with every +! resolution, and easily adaptable for dynamico ! For visualisation : > /u/tnalmd/bin/watercaps gcm_txt_output_file - INTEGER,SAVE :: icelocationmode = 2 + INTEGER,SAVE :: icelocationmode = 4 @@ -430,9 +433,18 @@ end do ! of (klon_glo) + ELSE IF (icelocationmode .eq. 4) THEN + + print*,'icelocationmode = 4' + print*,'Surfini: ice caps defined using manual 64x48 settings' + print*,'(although, it should work with any resolution)' + call locate_watercaptag(klon_glo,lati_glo, + & long_glo,watercaptag_glo) + +! print*,'watercaptag_glo(:), ',watercaptag_glo(:) ELSE print*, 'In surfini.F, icelocationmode is ', icelocationmode - print*, 'It should be 1, 2 or 3.' + print*, 'It should be 1, 2, 3 or 4 (default is 4)' call abort_physic("surfini","wrong icelocationmode",1) @@ -448,4 +460,10 @@ print*,'surfini: ice water cap', lati_glo(ig)*180./pi, & long_glo(ig)*180./pi, ig +! write(1,*),ig, lati_glo(ig)*180./pi, +! & cell_area(ig) +! write(2,*), lati_glo(ig)*180./pi, +! & long_glo(ig)*180./pi,cell_area(ig) +! write(3,*), ig, lati_glo(ig)*180./pi, +! & long_glo(ig)*180./pi,cell_area(ig) endif enddo @@ -464,4 +482,122 @@ ! end of #else of #ifndef MESOSCALE #endif - - END +! END SUBROUTINE surfini(ngrid,piceco2,qsurf) + END !SUBROUTINE surfini(ngrid,piceco2,qsurf) + + SUBROUTINE locate_watercaptag(klon_glo,lati_glo, + & long_glo,watercaptag_glo) + + USE comcstfi_h, ONLY: pi + + integer, intent(in) :: klon_glo + real, intent(in) :: lati_glo(klon_glo) + real, intent(in) :: long_glo(klon_glo) + logical, intent(out) :: watercaptag_glo(klon_glo) + integer :: ig,i +! real, dimension(klon_glo,120) :: wcap + real, dimension(120,2) :: latedge + real, dimension(120,2) :: lonedge + +! In icelocationmode=2 there are 120 manually predefined grid points where +! watercaptag is true (for the 64x48 resolution). The grid cells corners +! coordinates in latitude and longitude are written below. With this +! routine, we check if the grid cell center is in between any of those +! points. If so, watercaptag = true. + + + + + latedge(:,1)=(/ + & 88.125, 84.375, 84.375, 84.375, 84.375, 84.375,84.375, 84.375, + & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, + & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, + & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, + & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, + & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, + & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, + & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, + & 84.375, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, + & 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, + & 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, + & 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, + & 80.625, 80.625, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, + & 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 73.125, + & 73.125, 73.125, 73.125, 73.125, 73.125, 73.125, 73.125, 73.125/) + + + latedge(:,2)=(/ + & 90. , 88.125, 88.125, 88.125, 88.125, 88.125,88.125, 88.125, + & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, + & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, + & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, + & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, + & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, + & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, + & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, + & 88.125, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, + & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, + & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, + & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, + & 84.375, 84.375, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, + & 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 76.875, + & 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875/) + + + lonedge(:,1)=(/ + &-180. , -180. , -177.1875, -171.5625,-165.9375, -160.3125, + &-154.6875, -149.0625, -143.4375, -137.8125, -132.1875,-126.5625, + &-120.9375, -115.3125, -109.6875, -104.0625, -98.4375, -92.8125, + & -87.1875, -81.5625, -75.9375, -70.3125, -64.6875, -59.0625, + & -53.4375, -47.8125, -42.1875, -36.5625, -30.9375, -25.3125, + & -19.6875, -14.0625, -8.4375, -2.8125, 2.8125, 8.4375, + & 14.0625, 19.6875, 25.3125, 30.9375, 36.5625, 42.1875, + & 47.8125, 53.4375, 59.0625, 64.6875, 70.3125, 75.9375, + & 81.5625, 87.1875, 92.8125, 98.4375, 104.0625, 109.6875, + & 115.3125, 120.9375, 126.5625, 132.1875, 137.8125, 143.4375, + & 149.0625, 154.6875, 160.3125, 165.9375, 171.5625,-132.1875, + &-126.5625, -120.9375, -115.3125, -109.6875, -104.0625, -98.4375, + & -92.8125, -87.1875, -81.5625, -75.9375, -30.9375, -25.3125, + & -19.6875, -14.0625, -8.4375, -2.8125, 2.8125, 8.4375, + & 14.0625, 19.6875, 25.3125, 30.9375, 36.5625, 42.1875, + & 47.8125, 53.4375, 59.0625, 64.6875, 70.3125, 75.9375, + & 81.5625, 87.1875, -149.0625, -137.8125, -126.5625,-115.3125, + & -8.4375, 2.8125, 14.0625, 115.3125, 126.5625, 137.8125, + & 149.0625, 160.3125, 171.5625, -180. , -132.1875,-109.6875, + & 98.4375, 109.6875, 132.1875, 143.4375, 154.6875,165.9375/) + + lonedge(:,2)=(/ + & 180. , -180. , -171.5625, -165.9375,-160.3125, -154.6875, + &-149.0625,-143.4375, -137.8125, -132.1875, -126.5625, -120.9375, + &-115.3125,-109.6875, -104.0625, -98.4375, -92.8125, -87.1875, + & -81.5625, -75.9375, -70.3125, -64.6875, -59.0625, -53.4375, + & -47.8125, -42.1875, -36.5625, -30.9375, -25.3125, -19.6875, + & -14.0625, -8.4375, -2.8125, 2.8125, 8.4375, 14.0625, + & 19.6875, 25.3125, 30.9375, 36.5625, 42.1875, 47.8125, + & 53.4375, 59.0625, 64.6875, 70.3125, 75.9375, 81.5625, + & 87.1875, 92.8125, 98.4375, 104.0625, 109.6875, 115.3125, + & 120.9375, 126.5625, 132.1875, 137.8125, 143.4375, 149.0625, + & 154.6875, 160.3125, 165.9375, 171.5625, 177.1875, -126.5625, + &-120.9375,-115.3125, -109.6875, -104.0625, -98.4375, -92.8125, + & -87.1875, -81.5625, -75.9375, -70.3125, -25.3125, -19.6875, + & -14.0625, -8.4375, -2.8125, 2.8125, 8.4375, 14.0625, + & 19.6875, 25.3125, 30.9375, 36.5625, 42.1875, 47.8125, + & 53.4375, 59.0625, 64.6875, 70.3125, 75.9375, 81.5625, + & 87.1875, 92.8125, -143.4375, -132.1875, -120.9375, -109.6875, + & -2.8125, 8.4375, 19.6875, 120.9375, 132.1875, 143.4375, + & 154.6875, 165.9375, 177.1875, -177.1875, -126.5625, -104.0625, + & 104.0625, 115.3125, 137.8125, 149.0625, 160.3125,171.5625/) + + + watercaptag_glo(:) = .false. + DO ig=1, klon_glo + DO i=1, 120 + if ((long_glo(ig)*180./pi.ge.lonedge(i,1)) + & .and.(long_glo(ig)*180./pi.le.lonedge(i,2)) + & .and.(lati_glo(ig)*180./pi.ge.latedge(i,1)) + & .and.(lati_glo(ig)*180./pi.le.latedge(i,2))) then + watercaptag_glo(ig) = .true. + endif + ENDDO !i=1, 120 + ENDDO ! ig=1, klon_glo + + END SUBROUTINE locate_watercaptag