source: trunk/LMDZ.MARS/libf/phymars/surfini.F @ 706

Last change on this file since 706 was 697, checked in by emillour, 12 years ago

Mars GCM:

  • Minor cleanup in surfini.F
  • Corrected the polar mesh surface area which was wrong in the physics (changes in phyetat0.F, calfis.F and newstart.F)

EM
File size: 17.7 KB
Line 
1      SUBROUTINE surfini(ngrid,piceco2,qsurf,psolaralb)
2   ! to use  'getin'
3      USE ioipsl_getincom
4      use netcdf
5      IMPLICIT NONE
6c=======================================================================
7c
8c   creation des calottes pour l'etat initial
9c
10c=======================================================================
11c-----------------------------------------------------------------------
12c   Declarations:
13c   -------------
14#include "dimensions.h"
15#include "dimphys.h"
16#include "surfdat.h"
17#include "callkeys.h"
18#include "tracer.h"
19#include "comgeomfi.h"
20#include "comcstfi.h"
21
22#include "datafile.h"
23
24      INTEGER ngrid,ig,icap,iq,alternate
25      REAL  piceco2(ngrid),psolaralb(ngrid,2)
26      REAL qsurf(ngrid,nqmx) !tracer on surface (kg/m2)
27      REAL icedryness ! ice dryness
28     
29      ! longwatercaptag is watercaptag. Trick for some compilers
30      LOGICAL, DIMENSION(100000) :: longwatercaptag
31
32      EXTERNAL ISMIN,ISMAX
33      INTEGER ISMIN,ISMAX
34     
35! There are 3 different modes for ice distribution:
36! icelocationmode = 1 ---> based on data from surface.nc
37! icelocationmode = 2 ---> directly predefined for GCM resolutions 32x24 or 64x48
38! icelocationmode = 3 ---> based on logical relations for latitude and longitude
39! For visualisation : > /u/tnalmd/bin/watercaps gcm_txt_output_file
40      INTEGER,SAVE :: icelocationmode = 2
41       
42       
43      !in case icelocationmode == 1
44      INTEGER i,j
45      INTEGER     imd,jmd
46      PARAMETER   (imd=360,jmd=180)
47      REAL        zdata(imd*jmd)
48      REAL        zelat,zelon
49! on a data(i,j) = zdata((j-1)*jmd + i) avec i longitude de 1 a imd (E vers O) et j lat de 1 a jmd (N vers S)
50
51      INTEGER nb_ice(ngrid,2)              ! number of counts | detected ice for GCM grid
52      INTEGER latice(jjm,2),lonice (iim,2) ! number of counts | detected ice along lat & lon axis
53
54      REAL step,count,ratiolat
55
56      INTEGER   ierr,nid,nvarid
57     
58      REAL,SAVE :: min_icevalue = 500.
59      character(len=50) :: string = 'thermal'
60     
61      character (len=100) :: zedatafile
62c
63c=======================================================================
64
65c     water ice outliers
66c     ------------------------------------------
67
68      IF ((water) .and. (caps)) THEN
69     
70c Perennial H20 north cap defined by watercaptag=true (allows surface to be
71c hollowed by sublimation in vdifc).
72
73c We might not want albedodat to be modified because it is used to write
74c restart files. Instead, albedo is directly modified when needed (i.e.
75c if we have watercaptag and no co2 ice), below and in albedocaps.F90
76
77c       "Dryness coefficient" controlling the evaporation and
78c        sublimation from the ground water ice (close to 1)
79c        HERE, the goal is to correct for the fact
80c        that the simulated permanent water ice polar caps
81c        is larger than the actual cap and the atmospheric
82c        opacity not always realistic.
83
84         alternate = 0
85         
86         if (ngridmx .ne. 1) then
87           watercaptag(:) = .false.
88           longwatercaptag(:) = .false.
89         endif
90         
91         write(*,*) "Ice dryness ?"
92         icedryness=1. ! default value
93         call getin("icedryness",icedryness)
94         write(*,*) " icedryness = ",icedryness
95         dryness (:) = icedryness
96         
97       
98#ifdef MESOSCALE
99
100      do ig=1,ngridmx
101
102         !write(*,*) "all qsurf to zero. dirty."
103         do iq=1,nqmx
104         qsurf(ig,iq)=0.  !! on jette les inputs GCM
105                          !! on regle juste watercaptag
106                          !! il faudrait garder les inputs GCM
107                          !! si elles sont consequentes
108         enddo
109         if ( ( lati(ig)*180./pi .gt. 70. ) .and.
110     .        ( albedodat(ig) .ge. 0.26   ) )  then
111                 write(*,*)"outlier ",ig
112                 watercaptag(ig)  = .true.
113                 dryness(ig)      = 1.
114                 albedodat(ig)    = albedo_h2o_ice  !! pour output
115         else
116                 watercaptag(ig)  = .false.
117                 dryness(ig)      = 1.
118         endif
119         
120      enddo
121#else
122
123
124
125      IF (ngridmx .eq. 1) THEN ! special case for 1d --> do nothing
126     
127         print*, 'ngridmx = 1, do no put ice caps in surfini.F'
128
129      ELSE IF (icelocationmode .eq. 1) THEN
130     
131         print*,'Surfini: ice caps defined from surface.nc'
132           
133! This method detects ice as gridded value above min_icevalue in the field "string" from surface.nc
134! Typically, it is for thermal inertia above 500 tiu.
135! Two conditions are verified:
136! 1. GCM ice caps are defined such as area is conserved for a given latitude
137! (the approximation is that all points within the GCM latitude resolution have the same area).
138! 2. caps are placed to fill the GCM points with the most detected ice first.
139     
140
141           
142        zedatafile = trim(datafile)
143 
144       
145        ierr=nf90_open(trim(zedatafile)//'/surface.nc',
146     &  NF90_NOWRITE,nid)
147     
148      IF (ierr.NE.nf90_noerr) THEN
149       write(*,*)'Error : cannot open file surface.nc '
150       write(*,*)'(in phymars/surfini.F)'
151       write(*,*)'It should be in :',trim(zedatafile),'/'
152       write(*,*)'1) You can set this path in the callphys.def file:'
153       write(*,*)'   datadir=/path/to/the/datafiles'
154       write(*,*)'2) If necessary, surface.nc (and other datafiles)'
155       write(*,*)'   can be obtained online on:'
156       write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
157       CALL ABORT
158      ENDIF
159     
160     
161      ierr=nf90_inq_varid(nid, string, nvarid)
162      if (ierr.ne.nf90_noerr) then
163        write(*,*) 'surfini error, cannot find ',trim(string)
164        write(*,*) ' in file ',trim(zedatafile),'/surface.nc'
165        write(*,*)trim(nf90_strerror(ierr))
166        stop
167      endif
168
169      ierr=nf90_get_var(nid, nvarid, zdata)
170
171      if (ierr.ne.nf90_noerr) then
172        write(*,*) 'surfini: error failed loading ',trim(string)
173        write(*,*)trim(nf90_strerror(ierr))
174        stop
175      endif
176 
177                     
178      ierr=nf90_close(nid)
179 
180
181      nb_ice(:,1) = 1 ! default: there is no ice
182      latice(:,1) = 1
183      lonice(:,1) = 1
184      nb_ice(:,2) = 0
185      latice(:,2) = 0
186      lonice(:,2) = 0
187      !print*,'jjm,iim',jjm,iim ! jjm =  nb lati , iim = nb longi
188
189      ! loop over the GCM grid - except for poles (ig=1 and ngridmx)
190      do ig=2,ngridmx-1
191     
192        ! loop over the surface file grid
193        do i=1,imd*jmd
194         zelat = 90. - 90./real(jmd) - (i-1)/imd      ! surface file latitude
195         zelon = mod(i-1,imd) - (180.-180./real(imd)) ! surface file longitude
196                 
197          if ((abs(lati(ig)*180./pi-zelat) .le. 90./real(jjm)) .and.
198     &        (abs(long(ig)*180./pi-zelon) .le. 180./real(iim))) then
199              ! count all points in that GCM grid point
200              nb_ice(ig,1) = nb_ice(ig,1) + 1
201              if (zdata(i) > min_icevalue)
202                 ! count all detected points in that GCM grid point
203     &           nb_ice(ig,2) = nb_ice(ig,2) + 1
204          endif
205         
206        enddo ! of do i=1,imd*jmd
207       
208        ! projection of nb_ice on GCM lat and lon axes
209        latice(1+(ig-2)/iim,:) =
210     &     latice(1+(ig-2)/iim,:) + nb_ice(ig,:)
211        lonice(1+mod(ig-2,iim),:) =
212     &     lonice(1+mod(ig-2,iim),:) + nb_ice(ig,:) ! lonice is USELESS ...
213
214      enddo ! of do ig=2,ngridmx-1
215     
216
217     
218      ! special case for poles
219      nb_ice(1,2)   = 1  ! ice prescribed on north pole
220      latice(1,:)   = nb_ice(1,:)
221      lonice(1,:)   = nb_ice(1,:)
222      latice(jjm,:) = nb_ice(ngridmx,:)
223      lonice(iim,:) = nb_ice(ngridmx,:)
224     
225     
226  !    print*, 'latice TOT', latice(:,1)
227  !    print*, 'latice FOUND', latice(:,2)
228  !    print*, 'lonice TOT', lonice(:,1)
229  !    print*, 'lonice FOUND', lonice(:,2)
230     
231  !    print*, 'lat ratio', int(real(latice(:,2))/real(latice(:,1))*iim)
232  !    print*, 'lon ratio', int(real(lonice(:,2))/real(lonice(:,1))*jjm)
233     
234  !    print*,''
235  !    print*,'sum lat', sum(latice(:,1)), sum(lonice(:,1))
236  !    print*,'sum lon', sum(latice(:,2)), sum(lonice(:,2))
237     
238   
239      ! loop over GCM latitudes. CONSIDER ONLY NORTHERN HEMISPHERE
240      do i=1,jjm/2
241      step  = 1. ! threshold to add ice cap
242      count = 0. ! number of ice GCM caps at this latitude
243      ! ratiolat is the ratio of area covered by ice within this GCM latitude range
244      ratiolat  = real(latice(i,2))/real(latice(i,1))
245      !print*,'i',i,(i-1)*iim+2,i*iim+1
246     
247        ! put ice caps while there is not enough ice,
248        ! as long as the threshold is above 20%
249        do while ( (count .le. ratiolat*iim ) .and. (step .ge. 0.2))
250          count = 0.
251          ! loop over GCM longitudes
252          do j=1,iim
253            ! if the detected ice ratio in the GCM grid point
254            ! is more than 'step', then add ice
255            if (real(nb_ice((i-1)*iim+1+j,2))
256     &        / real(nb_ice((i-1)*iim+1+j,1)) .ge. step) then
257                  watercaptag((i-1)*iim+1+j) = .true.
258                  count = count + 1
259            endif
260          enddo ! of do j=1,iim
261          !print*, 'step',step,count,ratiolat*iim
262          step = step - 0.01
263        enddo ! of do while
264      !print*, 'step',step,count,ratiolat*iim
265
266      enddo ! of do i=1,jjm/2
267           
268
269      ELSE IF (icelocationmode .eq. 2) THEN
270     
271        print*,'Surfini: predefined ice caps'
272     
273        if ((iim .eq. 32) .and. (jjm .eq. 24)) then ! 32x24
274           
275          print*,'water ice caps distribution for 32x24 resolution'
276          longwatercaptag(1:9)    = .true. ! central cap - core
277          longwatercaptag(26:33)  = .true. ! central cap
278          longwatercaptag(1:33)  = .true. ! central cap
279          longwatercaptag(56)  = .true. ! central cap
280          longwatercaptag(58)  = .true. ! central cap
281          longwatercaptag(60)  = .true. ! central cap
282          longwatercaptag(62)  = .true. ! central cap
283          longwatercaptag(64)  = .true. ! central cap
284!---------------------   OUTLIERS  ----------------------------
285
286        else if ((iim .eq. 64) .and. (jjm .eq. 48)) then ! 64x48
287
288          print*,'water ice caps distribution for 64x48 resolution'
289          longwatercaptag(1:65)   = .true. ! central cap - core
290          longwatercaptag(75:85)  = .true. ! central cap
291          longwatercaptag(93:114) = .true. ! central cap
292!---------------------   OUTLIERS  ----------------------------
293          if (.true.) then
294          longwatercaptag(136)    = .true. ! outlier, lat = 78.75
295          longwatercaptag(138)    = .true. ! outlier, lat = 78.75
296          longwatercaptag(140)    = .true. ! outlier, lat = 78.75
297          longwatercaptag(142)    = .true. ! outlier, lat = 78.75
298          longwatercaptag(161)    = .true. ! outlier, lat = 78.75
299          longwatercaptag(163)    = .true. ! outlier, lat = 78.75
300          longwatercaptag(165)    = .true. ! outlier, lat = 78.75
301          longwatercaptag(183)    = .true. ! outlier, lat = 78.75
302          longwatercaptag(185)    = .true. ! outlier, lat = 78.75
303          longwatercaptag(187)    = .true. ! outlier, lat = 78.75
304          longwatercaptag(189)    = .true. ! outlier, lat = 78.75
305          longwatercaptag(191)    = .true. ! outlier, lat = 78.75
306          longwatercaptag(193)    = .true. ! outlier, lat = 78.75
307          longwatercaptag(194)    = .true. ! outlier, lat = 75
308          longwatercaptag(203)    = .true. ! outlier, lat = 75
309          longwatercaptag(207)    = .true. ! outlier, lat = 75
310          longwatercaptag(242)    = .true. ! outlier, lat = 75
311          longwatercaptag(244)    = .true. ! outlier, lat = 75
312          longwatercaptag(246)    = .true. ! outlier, lat = 75
313          longwatercaptag(248)    = .true. ! outlier, lat = 75
314          longwatercaptag(250)    = .true. ! outlier, lat = 75
315          longwatercaptag(252)    = .true. ! outlier, lat = 75
316          longwatercaptag(254)    = .true. ! outlier, lat = 75
317          longwatercaptag(256:257)= .true. ! outlier, lat = 75
318          endif
319!--------------------------------------------------------------       
320
321
322           
323        else if (ngridmx .ne. 1) then
324       
325      print*,'No predefined ice location for this resolution :',iim,jjm
326      print*,'Please change icelocationmode in surfini.F'
327      print*,'Or add some new definitions ...'
328      call abort
329         
330        endif
331
332        do ig=1,ngridmx
333          if (longwatercaptag(ig)) watercaptag(ig) = .true.
334        enddo
335
336
337      ELSE IF (icelocationmode .eq. 3) THEN
338     
339        print*,'Surfini: ice caps defined by lat and lon values'
340
341         do ig=1,ngridmx
342         
343c-------- Towards olympia planitia water caps -----------
344c--------------------------------------------------------
345
346       if ( ( ( lati(ig)*180./pi .ge. 77.  ) .and. ! cap #2
347     .           ( lati(ig)*180./pi .le. 80.  ) .and.
348     .           ( long(ig)*180./pi .ge. 110. ) .and.
349     .           ( long(ig)*180./pi .le. 181. ) )
350     .         .or.
351
352     .         ( ( lati(ig)*180./pi .ge. 75.  ) .and. ! cap #4 (Korolev crater)
353     .           ( lati(ig)*180./pi .le. 76.  ) .and.
354     .           ( long(ig)*180./pi .ge. 150. ) .and.
355     .           ( long(ig)*180./pi .le. 168. ) )
356     .         .or.
357     .         ( ( lati(ig)*180./pi .ge. 77 ) .and. ! cap #5
358     .           ( lati(ig)*180./pi .le. 80.  ) .and.
359     .           ( long(ig)*180./pi .ge. -150.) .and.
360     .           ( long(ig)*180./pi .le. -110.) ) )
361     .         then
362             
363               if ((alternate .eq. 0)) then  ! 1/2 en 64x48 sinon trop large en lat
364              !    watercaptag(ig)=.true.
365                  alternate = 1
366               else
367                  alternate = 0
368               endif !end if alternate = 0
369               
370       endif
371
372c----------- Opposite olympia planitia water cap --------
373c--------------------------------------------------------
374
375        if ( ( ( lati(ig)*180./pi     .ge.  80 ) .and.
376     .         ( lati(ig)*180./pi     .le.  84 ) )
377     .         .and.
378     .       ( ( long(ig)*180./pi .lt. -95. ) .or.       !!! 32x24
379     .         ( long(ig)*180./pi .gt.  85. ) ) ) then   !!! 32x24
380!     .     ( ( ( long(ig)*180./pi .ge. -29. ) .and.       !!! 64x48
381!     .         ( long(ig)*180./pi .le.  90. ) ) .or.      !!! 64x48
382!     .       ( ( long(ig)*180./pi .ge. -77. ) .and.       !!! 64x48
383!     .         ( long(ig)*180./pi .le. -70. ) ) ) ) then  !!! 64x48
384        !   watercaptag(ig)=.true.
385        endif
386
387
388c -------------------- Central cap ----------------------
389c--------------------------------------------------------
390
391      if (abs(lati(ig)*180./pi).gt.80)
392     .       watercaptag(ig)=.true.
393           
394c--------------------------------------------------------
395c--------------------------------------------------------
396      end do ! of (ngridmx)
397
398
399       ELSE
400     
401         print*, 'In surfini.F, icelocationmode is ', icelocationmode
402         print*, 'It should be 1, 2 or 3.'
403         call abort
404
405       ENDIF ! of if (icelocation)
406       
407       
408        ! print caps locations - useful for plots too
409        print*,'latitude | longitude | ig'
410        do ig=1,ngridmx
411          dryness (ig) = icedryness
412
413          if (watercaptag(ig)) then
414             print*,'ice water cap', lati(ig)*180./pi,
415     .              long(ig)*180./pi, ig
416          endif
417        enddo
418       
419#endif     
420
421       ENDIF ! (caps & water)
422       
423
424c ===============================================================
425c      INITIAL ALBEDO
426c ===============================================================
427
428         write(*,*)"surfini: water frost thickness",
429     s     frost_albedo_threshold
430         write(*,*)"surfini: water ice albedo:", albedo_h2o_ice
431         write(*,*)"surfini: water ice TI:", inert_h2o_ice
432
433c        To start with : Initial albedo = observed dataset
434c        -------------------------------------------------
435         DO ig=1,ngrid
436              psolaralb(ig,1)=albedodat(ig)
437              psolaralb(ig,2)=albedodat(ig)
438         END DO
439         PRINT*,'minimum albedo sans water caps',
440     s     albedodat(ISMIN(ngrid,albedodat,1))
441         PRINT*,'maximum albedo sans water caps',
442     s     albedodat(ISMAX(ngrid,albedodat,1))
443
444c        initial albedo if permanent H2O ice is present
445c        ------------------------------------------------
446         IF ((water) .and. (caps)) THEN
447           DO ig=1,ngrid
448            IF (watercaptag(ig)) THEN
449              psolaralb(ig,1) = albedo_h2o_ice
450              psolaralb(ig,2) = albedo_h2o_ice
451            ENDIF
452           END DO
453           PRINT*,'minimum albedo avec water caps',
454     s     psolaralb(ISMIN(ngrid,psolaralb,1),1)
455           PRINT*,'maximum albedo avec water caps',
456     s     psolaralb(ISMAX(ngrid,psolaralb,1),1)
457
458         ENDIF
459
460c      changing initial albedo if CO2 ice is present
461c      -------------------------------------------
462
463       DO ig=1,ngrid
464         IF (piceco2(ig) .GT. 0.) THEN
465             IF(ig.GT.ngrid/2+1) THEN
466                icap=2
467             ELSE
468                icap=1
469             ENDIF
470             psolaralb(ig,1) = albedice(icap)
471             psolaralb(ig,2) = albedice(icap)   
472         END IF
473       END DO
474
475c      changing initial albedo if water ice frost is present
476c      -------------------------------------------
477       IF (water) THEN
478          do iq=1,nqmx
479c          if there is frost and surface albedo is set to albedo_h2o_ice
480           if(noms(iq).eq."h2o_ice") then
481             do ig=1,ngrid
482              if ((piceco2(ig) .eq. 0.).and.
483     &          (qsurf(ig,iq).gt.frost_albedo_threshold)) then
484                     psolaralb(ig,1) = albedo_h2o_ice
485                     psolaralb(ig,2) = albedo_h2o_ice
486c                     PRINT*,'surfini.F frost',
487c     &                  lati(ig)*180./pi, long(ig)*180./pi
488               endif
489              enddo
490           endif
491          end do
492          PRINT*,'minimum albedo avec givre et co2',
493     s     psolaralb(ISMIN(ngrid,psolaralb,1),1)
494          PRINT*,'maximum albedo avec givre et co2',
495     s     psolaralb(ISMAX(ngrid,psolaralb,1),1)
496       END IF
497         
498      RETURN
499      END
Note: See TracBrowser for help on using the repository browser.