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

Last change on this file since 1009 was 740, checked in by tnavarro, 12 years ago

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